## Abstract

We present a mathematical model of biofilm response to antibiotics, controlled by a quorum sensing system that confers increased resistance. The model is a highly nonlinear system of partial differential equations that we investigate in computer simulations. Our results suggest that an adaptive, quorum sensing-controlled, mechanism to switch between modes of fast growth with little protection and protective modes of slow growth may confer benefits to biofilm populations. It enhances the formation of micro-niches in the inner regions of the biofilm in which bacteria are not easily reached by antibiotics. Whereas quorum sensing inhibitors can delay the onset of increased resistance, their advantage is lost after up-regulation. This emphasizes the importance of timing for treatment of biofilms with antibiotics.

This is a preview of subscription content, access via your institution.

## References

Abel zur Wiesch P, Clarelli F, Cohen T (2017) Using chemical reaction kinetics to predict optimal antibiotic treatment strategies. PLoS Comput Biol 13(1):e1005321. https://doi.org/10.1371/journal.pcbi.1005321

Anderl JN, Franklin MJ, Stewart PS (2000) Role of antibiotic penetration limitation in

*Klebsiella pneumoniae*biofilm resistance to ampicillin and ciprofloxacin. Antimicrob Agents Chemother 4(7):1818–1824Balaban NQ, Merrin J, Chait R, Kowalik L, Leibler S (2004) Bacterial persistence as a phenotypic switch. Science 305:1622–1625

Brackman G, Cos P, Maes L, Nelis HJ, Coenye T (2011) Quorum sensing inhibitors increase the susceptibility of bacterial biofilms to antibiotics in vitro and in vivo. Antimicrob Agents Chemother 55:2655–2661

Brackman G, Coenye T (2015) Quorum sensing inhibitors as anti-biofilm agents. Curr Pharm Des 21(1):5–11

Brown MR et al (1988) Resistance of bacterial biofilms to antibiotics: a growth-rate related effect? J Antimicrob Chemother 22:777–780

Bruchmann J, Kirchen S, Schwartz T (2013) Sub-inhibitory concentrations of antibiotics and wastewater influencing biofilm formation and gene expression of multi-resistant

*Pseudomonas aeruginosa*wastewater isolates. Environ Sci Pollut Res 20:3539–3549Chambless JD, Hunt SM, Stewart PS (2006) A three-dimensional computer model of four hypothetical mechanisms protecting biofilms from antimicrobials. Appl Environ Microbiol 72(3):2005–2013

Chopp DL, Kirisits MJ, Parsek MR, Moran B (2002) A mathematical model of quorum sensing in a growing

*P. aeruginosa*biofilm. J Ind Microbiol Biotechnol 29(6):339–346Chopp DL, Kirisits MJ, Parsek MR, Moran B (2003) The dependence of quorum sensing on the depth of a growing biofilm. Bull Math Biol 65(6):1053–1079

Cogan NG, Cortez R, Fauci L (2005) Modeling physiological resistance in bacterial biofilms. Bull Math Biol 67(4):831–853

Cogan NG (2008) Two-fluid model of biofilm disinfection. Bull Math Biol 70:800–819

Cogan NG, Szomolay B, Dindos M (2013) Effect of periodic disinfection on persisters in a one-dimensional biofilm model. Bull Math Biol 75:94–123

Cogan NG, Rath H, Kommerein N, Stumpp SN, Stiesch M (2016) Theoretical and experimental evidence for eliminating persister bacteria by manipulating killing timing. FEMS Microbiol Lett 363(23):fnw264

Czaran T, Hoekstra RF (2009) Microbial communication, cooperation and cheating: Quorum sensing dries the evolution of cooperation in bacteria. PLoS One 4:e6655

Demaret L, Eberl HJ, Efendiev M, Lasser R (2008) Analysis and simulation of a meso-scale model of diffusive resistance of bacterial biofilms to penetration of antibiotics. Adv Math Sci Appls 18(1):269–304

Dillon R, Fauci L, Fogelson A, Gaver D (1996) Modelling biofilm processes using the immersed boundary method. J Comput Phys 129(1):57–73

Eberl HJ, Parker DF, Van Loosdrecht CM (2001) A new deterministic spatio-temporal continuum model for biofilm development. J Theor Med 3:161–175

Eberl HJ, Demaret L (2007) A finite difference scheme for a degenerated diffusion equation arising in microbial ecology. El J Differ Equ CS 15:77–95

Eberl HJ, Sudarsan R (2008) Exposure of biofilms to slow flow fields: the convective contribution to growth and disinfection. J Theor Biol 253:788–807

Eberl HJ, Collinson S (2009) A modelling and simulation study of siderophore mediated antagonsim in dual-species biofilms. Theor Biol Med Mod 6:30

Efendiev MA, Zelik SV, Eberl HJ (2009) Existence and longtime behavior of a biofilm model. Commun Pure Appl Anal 8(2):509–531

Emerenini B, Hense BA, Kuttler C, Eberl HJ (2015) A mathematical model of quorum sensing induced biofilm detachment. Plos ONE 10(7):e0132385

Englmann M, Fekete A, Kuttler C, Frommberger M, Li X, Gebefügi I, Fekete J, Schmitt-Kopplin P (2007) The hydrolysis of unsubstituted N-acylhomoserine lactones to their homoserine metabolites. Analytical approaches using ultra performance liquid chromatography. J Chromatogr A 1160(1–2):184–93

Frederick M, Kuttler C, Hense BA, Müller J, Eberl HJ (2010) A mathematical model of quorum sensing in patchy biofilm communities with slow background flow. Can Appl Math Q 18(3):267–298

Frederick MR, Kuttler C, Hense BA, Eberl HJ (2011) A mathematical model of quorum sensing regulated EPS production in biofilms. Theor Biol Med Mod 8:8

Fuqua W, Winans S, Greenberg E (1994) Quorum sensing in bacteria: the LuxR-LuxI family of cell density-responsive transcriptional regulators. J Bacteriol 176:269–275

Ghasemi M, Eberl HJ (2017) Extension of a regularization based time-adaptive numerical method for a degenerate diffusion–reaction biofilm growth model to systems involving quorum sensing. Proc Comput Sci 108:1893–1902

Ghasemi M, Eberl HJ (2018) Time adaptive numerical solution of a highly degenerate diffusion–reaction biofilm model based on regularisation. J Sci Comput 74:1060–1090

Hall-Stoodley L, Costerton JW, Stoodley P (2004) Bacterial biofilms: from the natural environment to infectious diseases. Nat Rev Microbiol 2(2):95–108

Hense BA, Kuttler C, Müller J, Rothballer M, Hartman A, Kreft JU (2007) Does efficiency sensing unify diffusion and quorum sensing? Nat Rev Microbiol 5:230–239

Hense BA, Schuster M (2015) Core principles of bacterial autoinducer systems. Microbiol Mol Biol Rev 79:153–169

Hentzer M, Givskov M (2003) Pharmacological inhibition of quorum sensing for the treatment of chronic bacterial infections. J Clin Invest 112(9):1300–1307

Imran M, Smith H (2014) A model of optimal dosing of antibiotic treatment in biofilm. Math Biosci Eng 11(3):547–571

Janakiraman V, Englert D, Jayaraman A, Baskaran H (2009) Modeling growth and quorum sensing in biofilms grown in microfluidic chamber. Ann Biomed Eng 37(6):1206–1216

Kalia VC (2013) Quorum sensing inhibitors: an overview. Biotechnol Adv 31(2):224–245

Khassehkhan H, Hillen T, Eberl HJ (2009a) A non-linear master equation for a degenerate diffusion model of biofilm growth. LNCS 5544:735–744

Khassehkhan H, Efendiev MA, Eberl HJ (2009b) Existence and simulations of solutions to a degenerate diffusion–reaction model of an amensalistic biofilm control system. Discrete Contin Dyn Syst B 12(2):371–388

Ladyženskaja OA, Solonnikov A, Ural’ceva NN (1968) Linear and quasi-linear equations of parabolic type. AMS, Providence

Lear G, Lewis GD (2012) Microbial biofilms: current research and applications. Caister Academic Press, Poole. ISBN 978-1-904455-96-7

Lin LH, Wang JH, Yo JL, Li YY, Liu GX (2013) Effects of Allicin on the formation of

*Pseudomonas aeruginosa*biofilm and the production of Quorum-sensing controlled virulence factors. Pol J Microbiol 62:243–251Machineni L, Rajapantul A, Nandamuri V, Pawar PD (2017) Influence of nutrient availability and Quorum sensing on the formation of metabolically inactive microcolonies within structurally heterogeneous bacterial biofilms: an individual-based 3D cellular automata model. Bull Math Biol 79(3):594–618

Mah TC, O’Toole GA (2001) Mechanisms of biofilm resistance to antimicrobial agents. Trends Microbiol 9(1):34–39

Martins dos Santos VA, Yakimov MM, Timmis KN, Golyshin PN (2008) Genomic insights into oil biodegradation in marine systems. In Diaz E (ed) Microbial biodegradation: genomics and molecular biology. Horizon Scientific Press, Poole, p 1971. ISBN: 978-1-904455-17-2

Matur MG, Müller J, Kuttler C, Hense BA (2015) An approximative approach for single cell spatial modeling of quorum sensing. J Comput Biol 22:227–235

Muhammad N, Eberl HJ (2011) Model parameter uncertainties in a dual-species biofilm competition model affect ecological output parameters much stronger than morphological ones. Math Biosci 233(1):1–18

Müller J, Hense BA, Fuchs TM, Utz M, Pötzsche C (2013) Bet-hedging in stochastically switching environments. J Theor Biol 336:144–157

Mund A, Kuttler C, Perez-Velazquez J, Hense BA (2016) An age-dependent model to analyse the evolutionary stability of bacterial quorum sensing. J Theor Biol 405:104–115

Nealson K, Platt T, Hastings J (1970) Cellular control of the synthesis and activity of the bacterial luminescence system. J Bacteriol 104:313–322

Ngamsaad W, Sunatai S (2016) Mechanically-driven spreading of bacterial populations. Commun Nonlinear Sci Num Simul 35:88–96

Nielsen EI, Cars O, Friberg LE (2011) Pharmacokinetic/pharmacodynamic (PK/PD) indices of antibiotics predicted by a semimechanistic PKPD model: a step toward model-based dose optimization. Antimicrob Agents Chemother 55(10):4619–4630

Poulsen LV (1999) Microbial biofilm in food processing. LWT Food Sci Technol 32(6):321–326

Prakash B, Veeregowda B, Krishnappa G (2003) Biofilms: a survival strategy of bacteria. Curr Sci India 85:1299–1307

Rang J (2013) Improved traditional Rosenbrock–Wanner methods for stiff ODEs and DAEs. ‘. Institute of Scientific Computing, Germany

Redfield RJ (2002) Is quorum sensing a side effect of diffusion sensing? Trends Microbial 10:365–370

Ross-Gillespie A, Kümmerli R (2014) Collective decision-making in microbes. Front Microbiol 5:54

Saad Y (1994) SPARSKIT: a basic tool for sparse matrix computations. http://www.users.cs.umn.edu/saad/software/SPARSKIT/sparskit.html

Schwermer CU, Lavik G, Abed RM et al (2008) Impact of nitrate on the structure and function of bacterial biofilm communities in pipelines used for injection of seawater into oil fields. Appl Environ Microbiol 74(9):2841–51

Sengupta S, Chattopadhyay MK, Grossart HP (2013) The multifaceted roles of antibiotics and antibiotic resistance in nature. Front Micriobiol 4:47

Sonner S, Efendiev MA, Eberl HJ (2011) On the well-posedness of a mathematical model of Quorum-sensing in patchy biofilm communities. Math Methods Appl Sci 34(13):1667–1684

Stewart PS, Raquepas JB (1995) Implications of reaction–diffusion theory for disinfection of microbial biofilms by reactive antimicrobial agents. Chem Eng Sci 50(19):3099–3104

Stewart PS, Hamilton MA, Goldstein BR, Schneider BT (1996) Modelling biocide action against biofilms. Biotech Bioeng 49:445–455

Stewart PS (1996) Theoretical aspects of antimicrobial diffusion into microbial biofilms. Antimicrob Agents Chemother 40(11):2517–2522

Stewart PS, Costerton JE (2001) Antibiotic resistance of bacteria in biofilms. Lancet 358:135–38

Stewart PS (2003) Diffusion in biofilms. J Bacteriol 185:1485–1491

Stewart PS, Davison WM, Steenbergen JN (2009) Daptomycin rapidly penetrates a staphylococcus epidermidis biofilm. Antimicrob Agents Chemother 53(8):3505–3507

Szomolay B, Klapper I, Dindos M (2010) Analysis of adaptive response to dosing protocols for biofilm control. SIAM J Appl Math 70:3175–3202

Ulitzur S (1989) The regulatory control of the bacterial luminescence system—a new view. J Biolumin Chemilum 4:317–325

van Loosdrecht MCM, Heijnen JJ, Eberl H, Kreft J, Picioreanu C (2002) Mathematical modelling of biofilm structures. Antonie van Leeuwenhoek 81(1):245–256

Vaughan BL, Smith BG, Chopp DL (2010) The Influence of fluid flow on modeling Quorum sensing in bacterial biofilms. Bull Math Biol 72:1143–1165

Van der Vorst HA (1992) Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for the solution of non-symmetric linear systems. SIAM J Sci Stat Comput 13(2):631–644

Wang LL, Zhang CL, Gong FY, Li HT, Xie XH, Xia C, Chen J, Song Y, Shen AX, Song JX (2013) Influence of

*Pseudomonas aeruginosa*pvdQ gene on altering antibiotic susceptibility under swarming conditions. Curr Microbiol 66:152–161Wanner O, Eberl HJ, Van Loosdrecht MCM, Morgenroth E, Noguera DR, Picioreanu C, Rittmann BE (2006) Mathematical modelling of biofilms. IWA Publishing, London

Ward JP, King JR, Koerber AJ, Croft JM, Sockett RE, Williams P (2003) Early development and quorum sensing in bacterial biofilms. J Math Biol 47:23–55

Watnick P, Kolter R (2000) Biofilm—city of microbes (minireview). J Bacteriol 182(10):2675–2679

Williams P, Winzer K, Chan W, Cámara M (2007) Look who’s talking: communication and quorum sensing in the bacterial world. Philos Trans R Soc B 362:1119–1134

Yates EA et al (2002) N-acylhomoserine lactones undergo lactonolysis in a pH-, temperature-, and acyl chain length-dependent manner during growth of

*Yersinia pseudotuberculosis*and*Pseudomonas aeruginosa*. Infect Immun 70:5635–5646Yurtsev EA, Chao HX, Datta MS, Artemova T, Gore J (2013) Bacterial cheating drives the population dynamics of cooperative antibiotic resistance plasmids. Mol Sys Biol 9:683

Zhao J, Wang Q (2017) Three-dimensional numerical simulations of biofilm dynamics with quorum sensing in a flow cell. Bull Math Biol 79(4):884–919

## Acknowledgements

This study was supported in parts by the Natural Science and Engineering Research Council of Canada (NSERC) with a Discovery Grant and a Research Tools and Infrastructure Grant awarded to HJE, and a postgraduate scholarship awarded to MG. A visit of MG to the TU Munich for collaborative work was supported in parts by the Technical University of Munich with an Entrepreneurial Award for Global Challenges in The Mathematical Sciences, awarded to CK.

## Author information

### Affiliations

### Corresponding author

## Appendices

### Appendix A Solution Theory

For the biological and physical interpretation of the model, it is important that biomass densities and substrate concentrations remain nonnegative throughout and that the volume fraction of particulate biomass is bounded from above by unity. That the solutions of (1) have these properties, in particular the latter, is not obvious. To elucidate this, solutions of (1) are understood as weak solutions in the sense of distributions. Existence of solutions for boundary data

and initial data

where \(0<\delta _0 <1\), can be proven using the same arguments as in Demaret et al. (2008) for a simpler model with antibiotics without QS up-regulation of resistance, and in Khassehkhan et al. (2009a, b) and Muhammad and Eberl (2011) for two multi-species biofilm models. Existence and uniqueness of solutions has been proven for a model of quorum sensing induction in biofilms in Sonner et al. (2011), where the existence part in Sonner et al. (2011) follows the same line of argumentation as Demaret et al. (2008), but the uniqueness result relies on special properties of the reaction terms of the QS model that was studied there. While the QS sub-model of system (1), i.e., the model without antibiotic exposure, has the same properties, the reaction terms of the antibiotic sub-model prevent the application of these ideas. We point out, however, that there is no proof of non-uniqueness either. We have

### Theorem A.1

The system (1) with boundary conditions (3) and initial conditions (4) possesses a solution in the sense of distributions in \(L^\infty (\mathbb {R}_+ \times \Omega ) \times L^\infty (\mathbb {R}_+ \times \Omega ) \times L^\infty (\mathbb {R}_+\times \Omega ) \times L^\infty (\mathbb {R}_+ \times \Omega ) \times L^\infty (\mathbb {R}_+ \times \Omega ) \times L^\infty (\mathbb {R}_+ \times \Omega )\). This solution satisfies almost everywhere \(A \ge 0\), \(B\ge 0\), \(I\ge 0\) and \(A+B+I<1\), as well as \(0\le C\le C_\infty \), \(0\le N\le N_\infty \), \(0\le S\).

We refer for the technical details of the procedure to Demaret et al. (2008), Khassehkhan et al. (2009a, b), Muhammad and Eberl (2011), Sonner et al. (2011). The key arguments in the formal proof are a regularization of the biomass diffusion coefficient, i.e., *D*(*M*) in (1) is first replaced by

for some \(\epsilon >0\). The resulting regularized system is a non-degenerate quasilinear parabolic system, the existence of solutions of which follows with standard arguments, e.g., those in Ladyženskaja et al. (1968). It can then be shown that the solutions converge as \(\epsilon \rightarrow 0\) and that this limit is indeed a weak solution in the distributional sense of (1). Boundedness of biomass, \(M=A+B+I<1\), is hereby obtained by constructing an upper bound that is a solution to the underlying single-species growth model, for which boundedness by unity has been shown by Efendiev et al. (2009). Positivity of *I*, *A*, *B* is shown with standard comparison theorems. Hence, each individual biomass component is bounded by unity as well. The upper and lower bounds on the dissolved substrates follow from standard comparison theorems.

The boundary conditions (3) differ from (2) for the components *I*, *A*, *B* at the top boundary \(y=H\). Note, however, that where the biofilm/water interface does not reach the boundary, solutions that satisfy the homogeneous Neumann condition also satisfy the homogeneous Dirichlet condition, i.e., the solution of (1) with (2) and the solutions of (1) with (3) are identical until the biofilm reaches the top boundary, where nutrients and antibiotics are added and autoinducers removed. The biofilm reaching the top boundary corresponds to the scenario where the domain is completely clogged up with biomass, at which point additional physical processes come into play that are not accounted for in the model, i.e., the model breaks down. Indeed, for the underlying biofilm growth model, it was shown by Efendiev et al. (2009) that under homogeneous Neumann conditions on the biomass everywhere, as in (2), solutions might only exist for \(0<t<T_\mathrm{end}\) for some finite \(T_\mathrm{end}\), after the interface has reached the boundary everywhere, and \(M(T_\mathrm{end},\cdot )=1\) almost everywhere. On the other hand, if somewhere along the boundary Dirichlet or Robin conditions are posed that keep the biomass density there below unity, then the solutions exist for all \(t>0\). Further to this, the simulation results will be overshadowed by non-physical boundary condition effects before this happens, as discussed in Eberl and Collinson 2009. Therefore, our numerical solutions will always stop long before this occurs, i.e., the existence proof for (1) with (3) covers our simulation periods, as our solutions satisfy both (3) and (2) simultaneously.

### Appendix B Numerical Treatment

The two nonlinear diffusion effects (i) and (ii) described in Sect. 2.2 lead to two challenges for the numerical simulation: The biomass gradient blow-up at the biofilm-water interface \(\Gamma (t)\) is prone to introduce interface smearing effects or spurious oscillations, whereas the super-diffusion singularity of the biomass diffusion coefficient forces time-steps to be small in order to not overshoot the singularity, even for implicit methods. The method that we apply in our simulations is an extension to the model at hand of the method presented by Ghasemi and Eberl (2017a) for a similar model of quorum sensing in a biofilm, but without antibiotics, which in turn is an extension of a numerical method for the underlying prototype biofilm growth model that was introduced and analyzed previously (Ghasemi and Eberl 2017b). Spatial discretization is performed by a Finite Volume method on a uniform grid, whereas an error-controlled, time-adaptive Rosenbrock–Wanner method is used for time integration. It has been shown in Ghasemi and Eberl (2017a, b) that this method is able to address the numerical challenges stemming from (i) and (ii) described above and that it overcomes the shortcomings of fixed time-step methods that have been used previously for similar models (Khassehkhan et al. 2009; Muhammad and Eberl 2011).

For the numerical treatment, we non-dimensionalize the model (1) with choices \(\tilde{x}=x/L,~T=t\mu _A\) for the independent variables, here *L* is a characteristic length scale of the computational domain. Under our boundary conditions, the height of the domain is a measure for diffusion length. \(\frac{1}{\mu _A}\) is the characteristic timescale for growth of biomass species *A*. The concentration variables *N*, *C* and *S* are non-dimensionalized as \(\tilde{N}=\frac{N}{N_{\infty }},~\tilde{C}=\frac{C}{C_0}\) and \(\tilde{S}=\frac{S}{\tau }\) where \(N_{\infty }\) is the bulk substrate concentration and \(C_0=\frac{\delta _A}{\mu _A}\). Note that the volume fractions *I*, *A* and *B* are already defined as dimensionless variables. With these choices, the new model parameters are defined as:

The non-dimensionalized equations read:

Although we use the non-dimensional formulation for the numerical treatment, we will make our choices of parameters based on the dimensional values, and describe our subsequent simulation experiments in terms of those, for easier biological and physical interpretation. Going forward, we will drop the tilde from our notation in (6).

For spatial discretization, we introduce a uniform grid of size \(N \times M\) for the rectangular domain \([0,1]\times [0,H/L]\) and integrate each equation of system (6) over each grid cell. For instance, integrating the equation for *A* over grid cell with index (*i*, *j*) and using the Divergence Theorem gives

for \(i=1,\ldots ,N\) and \(j=1,\ldots ,M\), where \(v_{i,j}\) denotes the domain of the grid cell, \(J_n = D(M)\partial _nA\) represents the outward normal flux across the grid cell boundary, and \(R_{1_{A}}(N,C,S)=\frac{N}{K_{N}+N}-\beta _A\frac{C^{n_1}}{K_{{C}}^{n_1}+C^{n_1}}-\omega \frac{S^{n_2} }{1+S^{n_2}}-k_A\) and \(R_{2_{A}}(S)=\psi \frac{1}{1+S^{n_2}}\) stand for the reaction terms.

To evaluate the area integrals in (7), we evaluate the dependent variables at the center of the grid cells,

for \(i=1,\ldots ,N\) and \(j=1,\ldots ,M\) with \(\Delta x=1/N=H/LM\); the value of all other dependent variables at the center of the grid cells can be evaluated in a similar way. We approximate the integrals in (7) by the midpoint rule, and similarly the line integral is evaluated by considering every edge of the grid cell separately. To this end, the diffusion coefficient *D*(*M*) in the midpoint of the cell edge is approximated by arithmetic averaging from the neighboring grid cell center points, and the derivative of *A* across the cell edge by a central finite difference. We get then for the biomass fraction *A* in cell (*i*, *j*) the ordinary differential equation

where \(R_{1_{A_{i,j}}}=\frac{N_{i,j}}{K_{N}+N_{i,j}}-\beta _A\frac{C_{i,j}^{n_1}}{K_{{C}}^{n_1}+C_{i,j}^{n_1}}-\omega \frac{S_{i,j}^{n_2} }{1+S_{i,j}^{n_2}}-k_A\) and \(R_{2_{A_{i,j}}}=\psi \frac{1}{1+S_{i,j}^{n_2}}\) and for the fluxes we have, accounting for the boundary conditions,

The spatial discretization of the equations for the other dependent variables *B*, *I*, *N*, *C*, *S* follows the same principle. The major difference is for the substrates for which we have a Robin boundary condition at the top of the domain instead of homogeneous Neumann condition. We refer the reader to Ghasemi and Eberl (2017b) for a detailed description of the changes to the discretization that this introduces.

Introducing the lexicographical grid ordering

and the vector notation \(\mathbf I \!=\!(I_1,\ldots ,I_{NM})\), \(\mathbf A =(A_1,\ldots ,A_{NM})\), \(\mathbf B =(B_1,\ldots ,B_{NM})\), \(\mathbf N =(N_1,\ldots ,N_{NM})\), \(\mathbf C =(C_1,\ldots ,C_{NM})\) and \(\mathbf S =(S_1,\ldots ,S_{NM})\) with \((I,A,B,N,C, S)_{p}:=(I,A,B,N,C,S)_{\pi (i,j)}=(I,A,B,N,C,S)_{i,j}\) for \(i=1,\ldots ,N\), \(j=1,\ldots ,M\), yield the coupled system of \(6 \cdot N \cdot M\) ordinary differential equations

### Remark B.1

In ODE system (15), the matrices \( D _{I,A,B,N,C,S}\) are block matrices of size \(NM \times NM\) that depend on the dependent variables \(\mathbf I,A,B \). They are symmetric, and weakly diagonally dominant with non-positive main diagonals and nonnegative off-diagonals and contain the spatial derivative terms of each equation. They can efficiently be stored in sparse diagonal format with 4 off-diagonals, with offsets \(\pm 1, \pm N\). The matrices \( R _{{1,2}_{(I,A,B,N,C,S)}}\) are diagonal matrices of size \(NM \times NM\) which contain the reaction terms of each equation. The vectors \(\mathbf b _{I,A,B,N,C,S}\) are of size *NM* and contain contributions from the imposed boundary conditions for each biomass species and dissolved substrates. By the posed boundary conditions for biomass species, the entries of \(\mathbf b _{I,A,B}\) are zero [note that this would be different if Robin conditions were posed instead, as in (3)]. For substrates due to the Robin boundary conditions at the top, the entries \(\mathbf b _{N,C,S}\) are zero for all grid cells (*i*, *j*) with \(j<M\), and for \(j=M\) are obtained from the boundary conditions. With these properties, it is easy to show with standard comparison theorems that the solutions of (15) are nonnegative if the initial data are, cf. Ghasemi and Eberl (2017a, b). Boundedness of the biomass densities, \(\mathbf M =\mathbf A +\mathbf B +\mathbf I <1\), can be shown for sufficiently small spatial step-sizes by constructing an upper bound using a discretized version of the underlying growth model, cf. Ghasemi and Eberl (2017a).

For the time integration of the semi-discrete system (15), we use the time-adaptive, error-controlled embedded Rosenbrock–Wanner method ROS3PRL, which was introduced by Rang (2013). It has been shown in Ghasemi and Eberl (2017a) that this method, together with the above Finite Volume discretization is able to proceed with efficient time-steps for models of this type. While intermittently very small time-steps might be required to meet accuracy requirements of the method, previous experience has shown that it will recover quickly, to values much larger than would be required by fixed time-step methods. It also has been demonstrated in Ghasemi and Eberl (2017a) that the spatial resolution that is required for a sufficiently accurate resolution of the interface, without numerical smearing, is acceptable. The spatial resolution chosen in our simulations is based on a grid refinement study (data not shown), and has been determined as an acceptable trade-off between accuracy and computational cost.

The embedded Rosenbrock–Wanner method requires in each time-step the solution of several sparse, large, non-symmetirc linear systems, for which we use the stabilized bi-conjugate gradient method (Saad 1994; Van der Vorst 1992). The numerical method has been implemented in Fortran 95. OpenMP was used for the parallelization of selected computationally expensive tasks, such as the linear solver and the formation of the Jacobian matrix within the Rosenbrock–Wanner method, and for the evaluation of nonlinear reaction terms. All simulations reported here have been carried out on custom built multi-core and multi-processor Linux desktop workstations.

## Rights and permissions

## About this article

### Cite this article

Ghasemi, M., Hense, B.A., Eberl, H.J. *et al.* Simulation-Based Exploration of Quorum Sensing Triggered Resistance of Biofilms to Antibiotics.
*Bull Math Biol* **80, **1736–1775 (2018). https://doi.org/10.1007/s11538-018-0433-3

Received:

Accepted:

Published:

Issue Date:

### Keywords

- Antibiotics
- Biofilm
- Mathematical model
- Nonlinear diffusion
- Quorum sensing
- Simulation

### Mathematics Subject Classification

- 92D25
- 35K65
- 65N08
- 34C60