Dynamic energy budget approach to evaluate antibiotic effects on biofilms

TL;DR

This paper introduces a dynamic energy budget model to simulate and analyze the effects of antibiotics on biofilm survival, capturing spatial distributions and resistance mechanisms to aid therapy development.

Contribution

It presents a novel computational framework combining cell metabolism and concentration fields to predict biofilm responses to antibiotics.

Findings

Higher survival rates near substratum with active cell-targeting antibiotics

Increased polymer production with antibiotic treatment

Model aligns qualitatively with experimental observations

Abstract

Quantifying the action of antibiotics on biofilms is essential to devise therapies against chronic infections. Biofilms are bacterial communities attached to moist surfaces, sheltered from external aggressions by a polymeric matrix. Coupling a dynamic energy budget based description of cell metabolism to surrounding concentration fields, we are able to approximate survival curves measured for different antibiotics. We reproduce numerically stratified distributions of cell types within the biofilm and introduce ways to incorporate different resistance mechanisms. Qualitative predictions follow that are in agreement with experimental observations, such as higher survival rates of cells close to the substratum when employing antibiotics targeting active cells or enhanced polymer production when antibiotics are administered. The current computational model enables validation and hypothesis testing when developing therapies.

Tables

Table 1. Table 1: Variables of the model. The state variables define the model. The auxiliary variables are introduced for simplicity. Non-dimensional variables are labeled ’n.d.’.

Symbol	Units	Variable
t	s	time
x	m	space
Symbol	Units	State variable	Equation
$e$	n.d.	Scaled energy density	Eq. (1)
$v$	n.d.	Dimensionless bacterial volume	Eq. (2)
$l$	n.d.	Dimensionless bacterial length	Eq. (2)
$v_e$	n.d.	Dimensionless EPS volume	Eq. (3)
$\epsilon$	n.d.	Environmental degradation	Eq. (8)
$a$	n.d.	Acclimation energy density	Eq. (9)
$h$	${\mathrm{hr}}^{-1}$	Hazard rate	Eq. (10)
$q$	${\mathrm{hr}}^{-2}$	Aging acceleration	Eq. (11)
$C_o$	$\frac{\mathrm{mg}}{\ell }$	Oxygen concentration	Eq. (4)
$C_a$	$\frac{\mathrm{mg}}{\ell }$	Antibiotic concentration	Eq. (5)
$\left[C_{IN}\right]$	$\frac{\mathrm{mg}}{\mathrm{Cmol}}$	Antibiotic cellular density	Eq. (6)
Symbol	Units	Auxiliary variable	Equation
$f$	n.d.	Scaled functional response	Eq. (1)
$r$	${\mathrm{hr}}^{-1}$	Bacterial production rate	Eq. (2)
${\nu }^{\prime }$	${\mathrm{hr}}^{-1}$	Conductance modified by exposure	Eq. (7)

Table 2. Table 2: Parameter values used in the simulations. Recall that 1 1 1 Cmol = 12 12 12 g.

Symbol	Units	Name	Value	Source
${\rho }_x$	$\frac{\mathrm{mg}}{\ell }$	cell density	$47000$	[35]
	$\mu$m	biofilm size	$1000\times 200$	chosen
$L$	$\mu$m	cell length	$1.67-2.75$	[34]
$2\rho$	$\mu$m	cell diameter	$0.5-0.73$	[34]
$C_{o,out}$	$\frac{\mathrm{mg}}{\ell }$	Oxygen concentration	$0.035$	[35]
$K_o$	$\frac{\mathrm{mg}}{\ell }$	Oxygen half saturation	$0.1$	[35]
$d_o$	$\frac{\mu m^2}{s}$	Oxygen diffusion	$2.2\times {10}^4$	[42]
${\mu }_{max}$	${\mathrm{hr}}^{-1}$	Growth rate with oxygen	$0.3$	[35]
$Y_{x/o}$	$\frac{\mathrm{mg}\mathrm{cell}}{\mathrm{mg}\mathrm{oxygen}}$	cell yield	$0.24$	[35]
$k$	$\frac{\mathrm{mg}\mathrm{polymer}}{\mathrm{mg}\mathrm{cell}}$	growth associated yield	$2.2371$	estimated
$k^{\prime }$	$\frac{\mathrm{mg}\mathrm{polymer}}{\mathrm{mg}\mathrm{cell}\mathrm{hour}}$	non growth associated yield	$0.29$	estimated
$Y_{x/o}$	$\frac{\mathrm{mg}\mathrm{cell}}{\mathrm{mg}\mathrm{oxygen}}$	cell yield	$0.34$	estimated
$Y_{p/o}$	$\frac{\mathrm{mg}\mathrm{polymer}}{\mathrm{mg}\mathrm{oxygen}}$	polymer yield	$0.56$	estimated
$C_{a,out}$	$\frac{\mathrm{mg}}{\ell }$	Antibiotic concentration	$0.78,0.20,3.13,1.56$	[17]
$d_a$	$\frac{\mu m^2}{s}$	Antibiotic diffusion	$0.5\times {10}^4$	[35]
$d_e$	$\frac{\mu m^2}{s}$	EPS diffusion	${10}^4$	estimated
$d_{\epsilon }$	$\frac{\mu m^2}{s}$	Degradation diffusion	$0.5\times {10}^4$	estimated
$\nu$	${\mathrm{hr}}^{-1}$	Energy conductance	0.84768	[21]
$m_{\kappa }$	${\mathrm{hr}}^{-1}$	Maintenance rate	0.1266	estimated
$g_{\kappa }$	n.d.	Investment ratio	0.9766	estimated
$s_q$	$\frac{\ell }{\mathrm{mg}}$	Multiplicative stress coeff.	$0.8921\times {10}^{-5}$	estimated
$h_a$	${\mathrm{hr}}^{-2}$	Weibull aging acceleration	$1.4192\times {10}^{-4}$	estimated
$h_0$	${\mathrm{hr}}^{-1}$	Initial hazard rate	0.4	estimated
$k_{qA}^I$	$\frac{\mathrm{Cmol}}{\mathrm{mg}{\mathrm{hour}}^3}$	Dissolved antibiotic toxicity	$\frac{s}{C_{a,out}}\frac{k_A^O}{k_A^I}$	estimated
$k_A^I$	$\frac{\ell }{\mathrm{hr}\mathrm{Cmol}}$	Antibiotic influx coeff.	$8.6\times {10}^{-6}$	[21]
$k_A^O$	${\mathrm{hr}}^{-1}$	Antibiotic efflux coeff.	0.17251	[21]
$K_V$	$\frac{\mathrm{mg}}{\ell }$	Noncompetitive inhibition coeff.	154.82	[21]
$a_M$	n.d.	Target acclimation energy	1.6703	[21]
${\nu }_{\epsilon }$	$\frac{\ell }{\mathrm{mg}}$	Environmental degradation coeff.	$0.23566/12000$	[21]
${\nu }_m$	n.d.	Maintenance respiratory coefficient	$0.054703$	[21]
${\gamma }_{\epsilon }$	n.d.	Environmental degradation effect coeff.	$1$	[21]

Equations

$$\frac{de}{dt}={\nu }^{\prime }(f-e),f=\frac{C_o}{C_o+K_o},{\nu }^{\prime }={\nu }_A\nu ,$$

$$\frac{dv}{dt}=(r\frac{a}{a_M}-h)v,r={\left(\frac{{\nu }^{\prime }e-m_{\kappa }{g}_{\kappa }}{e+g_{\kappa }}\right)}^{+},$$

$$\frac{d{v}_e}{dt}=\frac{g_{\kappa }}{g_{e,\kappa }}(r\frac{a}{a_M}-h)v+\frac{m_{\kappa }}{g_{e,\kappa }}v=r_{e}v.$$

$$\frac{\partial C_o}{\partial t}=d_o\Delta C_o-\frac{{\mu }_{max}f}{Y_{x/o}}X-\frac{(k{\mu }_{max}f+k^{\prime })}{Y_{p/o}}X,$$

$$\frac{\partial C_a}{\partial t}=d_a\Delta C_a-(R_e+R)C_a,$$

$$\frac{d[C_{IN}]}{dt}=k_A^I{C}_a-k_A^O[C_{IN}],$$

$${\nu }^{\prime }={\nu }_A\nu ,{\nu }_A=e^{-{\gamma }_{\epsilon }\epsilon }{\left(1+\frac{C_a}{K_V}\right)}^{-1},$$

$$\frac{\partial \epsilon }{\partial t}=d_{\epsilon }\Delta \epsilon +{\nu }_{\epsilon }(R+{\nu }_m{m}_{\kappa })X.$$

$$\frac{da}{dt}=(r+r_e){\left(1-\frac{a}{a_M}\right)}^{+},$$

$$\frac{dh}{dt}=q-(r+r_e)h.$$

$$\frac{dq}{dt}=e(s_{q}Xq+h_a)({\nu }^{\prime }-r)+{\left(\frac{dq}{dt}\right)}_A-rq,$$

$${\left(\frac{dq}{dt}\right)}_A=k_{qA}^I[C_{IN}]-r_{e}q,$$

$$\frac{dp}{dt}=-hp,p(0)=1.$$

$$\frac{\partial C_e}{\partial t}=d_e\Delta C_e+(1-\alpha )R_{e}X,$$

$$\frac{d{X}_e}{dt}=[k{\mu }_{max}f+k^{\prime }]X,\frac{dX}{dt}={\mu }_{max}fX.$$

$$\frac{d{v}_e}{dt}=\left(\frac{g_{\kappa }}{g_{e,\kappa }}r+\frac{m_{\kappa }}{g_{e,\kappa }}\right)v,\frac{dv}{dt}=rv,r={\left(\frac{\nu e-m_{\kappa }{g}_{\kappa }}{e+g_{\kappa }}\right)}^{+}.$$

$$g_{\kappa}={1\over 2}\left(-{k\over k^{\prime}}\mu_{max}f_{max}+\sqrt{\big{(}{k\over k^{\prime}}\mu_{max}f_{max}\big{)}^{2}-4{k\over k^{\prime}}(\mu_{max}f_{max}-\nu)f_{max}}\right).$$

$$q(t)=\frac{c+s}{d-b}\left(e^{(d-b)t}-1\right),$$

$$h(t)=\left(a-\frac{c+s}{d-b}\left(\frac{1}{d}-\frac{1}{b}\right)\right)e^{-bt}+\frac{c+s}{d-b}\left(\frac{e^{(d-b)t}}{d}-\frac{1}{b}\right),$$

$$\displaystyle p(t)=exp\Bigg{(}{a\over b}\left(e^{-b\,t}-1\right)-{c+s\over d-b}\Bigg{(}\Big{(}{1\over b\,d}-{1\over b^{2}}\Big{)}e^{-b\,t}+$$

$$\displaystyle{e^{(d-b)t}\over d(d-b)}-{1\over b}t-{1\over b\,d}+{1\over b^{2}}-{1\over d(d-b)}\Bigg{)}\Bigg{)}.$$

$$C_{o;i,j}^{\ell +1}=C_{o;i,j}^{\ell }+\delta t{d}_o\frac{C_{o;i,j+1}^{\ell }+C_{o;i+1,j}^{\ell }-4C_{o;i,j}^{\ell }+C_{o;i,j-1}^{\ell }+C_{o;i-1,j}^{\ell }}{\delta x^2}+\delta tg(X_{i,j}),$$

Full text

††journal: Commun. Nonlinear Sci. Numer. Simul.

1 Introduction

Biofilms are bacterial aggregates that grow on moist surfaces, encased in a self-produced polymeric matrix, see Figure 1. The matrix creates a favorable environment for their development, facilitating nutrient, oxygen and waste transport [30]. It also acts as a shield against external aggressions by flows, disinfectants and antibiotics. The minimal bactericidal concentration (MBC) and minimal inhibitory concentration (MIC) of antibiotics to bacteria in their biofilm habitat may be up to 100-1000 fold higher compared with planktonic bacteria [1, 16].

Implant associated infections typically involve biofilm growth on the surface of the implant [40]. They form on medical equipment and prostheses, such as pacemakers and endotracheal tubes, central lines, intravenous catheters, stents and artificial joints. Bloodstream infections, and many other hospital-acquired infections, may be caused by them. Biofilms may also spread on body surfaces such as heart valves (endocarditis), teeth, the lungs of cystic fibrosis patients (pneumonia), the middle ear and nose (otitis, rhinosinusitis), bones (osteomyelitis) or in chronic wounds [16]. The biofilm matrix hinders phagocytosis and other actions of the inmune system. Bacteria surviving standard antimicrobial therapies are able to reproduce, originating chronic infections [38]. To tackle this problem, we must be able to understand how antibiotics act and how resistance to antibiotics develops.

Antibiotics affect cells in diverse ways [19]. $\beta$-lactams (penicillins, cepha-losporins, carbapenems) and glycopeptides (vancomycin) inhibit cell-wall synthesis. Aminoglycosides (streptomycin, gentamycin, tobramycin) inhibit protein synthesis. Quinolones (ciprofloxacin, ofloxacin) inhibit DNA replication. Tetracyclines inhibit translation. Polymyxins such as colistin disrupt charge distributions in the outer cell membrane. Eventually, the damage caused to the cell produces its death. To exert their antibacterial action, antibiotics undergo a certain process. They need to penetrate the cells, remain stable and accumulate to reach inhibitory concentrations. Sometimes they have to take an active form. After detecting their target, they interact with it to exert their action. Interferences in either of these processes may result in cell resistance to the antibiotic. Resistance typically proceeds through efflux systems (the antibiotic is pumped out of the cell), chemical alterations of the antibiotic (cellular enzymes degrade it), mutations in antibiotic target molecules, and non-heritable resistance caused by environmental conditions [9, 16, 19, 36, 37]. The main resistance mechanisms for different types of antibiotics are summarized in [19].

The biofilm environment enhances bacterial resistance in a number of ways. Biofilm development is influenced by quorum sensing [8]. Through quorum sensing mechanisms, bacteria sense when a critical number of them are present in the environment. They respond by activating genes that produce exopolysaccharides [7]. The polysaccharide matrix surrounding the bacterial community delays diffusion of antibiotics inside the biofilm. Nevertheless, direct measurements suggest that some antibiotics equilibrate within the biofilm [6] after a waiting time. Pseudomonas aeruginosa tends to be the main source of gram-negative infections in intensive care units in developed countries [32]. P. Aeruginosa and other bacteria express $\beta$-lactamase, an enzyme that attacks $\beta$-lactams. An enzyme breaking the antibiotic at a rate at which it crosses the cell membrane combined with delayed diffusion might explain resistance to penicillins, but not to other antibiotics [6].

As mentioned above, as we penetrate from the outer biofilm surface towards the interface with the substratum, gradients of oxygen and nutrients develop. Oxygen depletes [5, 42]. These gradients result in increased doubling times for cell division and reduced bacterial metabolic activity. The intensity of metabolic processes is stratified: high activity in the outer layers and slow growth or no growth in the inner core. These dormant cells are partially responsible for tolerance to antibiotics. Popular monotherapies with $\beta$-lactams are only active against dividing cells [3], forcing combinations with antibiotics that are active against nondividing cells, such as colistin [16]. Oxygen limitation and metabolic rates are also important factors enhancing the tolerance of biofilms to ciprofloxacin and aminoglycosides [41].

In the biofilm, bacteria are exposed to oxidative stress, that causes hipermutability. Enhanced production of reactive oxygen species (ROS), either released in response to the infection or produced by the alterations in the DNA repair system of the bacteria, leads to an environment with low oxygen tension filled with oxygen radicals [16, 29]. Augmented $\beta$-lactamase synthesis, overexpression of efflux-pumps and increased EPS production follow [4, 16, 26]. Studies with toxicants have also shown increased EPS production and ability to adapt to the toxicant, repairing damage to the cell [13]. Quorum sensing inhibitors [15], efflux inhibitors [27, 28], antioxidants reducing the oxidative stress and mutations [29], together with enzymes able to dissolve the biofilm matrix [2], may provide strategies to overcome resistance mechanisms.

Mathematical modeling can assist in the design of therapies and the interpretation of experimental data. Early models were able to reproduce elementary qualitative behavior. Reference [35] uses coupled reaction-diffusion equations for the concentration of oxygen, antibiotics and the volume fractions of live and dead cells to predict survival profiles inside thick biofilms due to slow growth. Reference [33] predicts that the biofilm matrix can not prevent diffusion of $\beta$-lactam antibiotics into the bacteria provided the amount of chromosomal $\beta$-lactamase is low. The diversity of the mechanisms involved in biofilms resistance to antibiotics suggests the opportunity of adapting dynamic energy budget (DEB) frameworks to describe the effect of antibiotics on them. DEB models have already been exploited to describe the effects of toxicant exposure on populations of floating bacteria [21]. Like antibiotics, toxicants interfere with the metabolism of cells and increase the energy required for cell maintenance. The cell requires additional energy to expel the toxicant and repair damage caused by toxicant activity (DNA, RNA, protein repair).

Dynamic energy budget models relate biomolecular processes to individual physiology and population dynamics [23]. They have been successfully used to study scaling behaviors of all sorts of living beings, from plants and animals to cells [14]. Basic DEB models describe acquisition of biomass and energy, as well as energy allocation for cell maintenance, growth and division. A standard description of aging and death processes is also available [23]. Surplus reactive oxygen species (ROS) causes irreparable damage to the cell. Damage components become ‘damage inducing compounds’, resulting in cell malfunctioning and further damage, raising the death probability. Death and damage are represented through ‘hazard’ and ‘aging acceleration’ variables. Antibiotic impact, like toxicant consequences, must be described modifying standard fluxes and rates. To represent the effect of toxicants on floating bacteria, Reference [21] included acclimation to the toxicant effect, environmental degradation due to cell products, toxicant induced ROS production, and their influence on assimilation, growth, aging and harzard rates. The only variable states characterizing the environment are substrate and toxicant concentrations.

In this paper, we analyze antibiotic effects on Pseudomonas aureginosa biofilms by means of an adapted DEB model. Unlike typical DEB models for cells [21], we introduce here spatial variations in the concentration equations to study the local influence of the biofilm shape and composition. Spatial gradients in the concentrations of oxygen and antibiotics are relevant to the bacterial survival profile in these communities. We couple numerically the metabolism of individual cells to the diffusion equations describing concentrations and some additional field. Our model accounts for EPS production, identified as a cause of antibiotic resistance. The antibiotic acts on growth, maintenance, aging and hazard rate in a similar way to toxicants. EPS production increases acclimation, while reducing growth, damage and hazard rates, see Figure 2.

We have fitted the model for a quantitative agreement with the death rates reported for P. aureginosa under different antibiotics [17] and explored the spatial variations in the dead/alive cell distribution. Our simulations reproduce expected qualitative trends. Augmenting the dose, the number of dead cells increases. Increasing the EPS extent, the number of dead cells decreases. Antibiotic presence enhances EPS production. When we apply antibiotics targeting active cells a necrotic region progressing from the outer biofilm layers to the inner biofilm core appears. Instead, antibiotics targeting slowly growing cells would destroy first the inner dormant core.

The paper is organized as follows. Section 2 describes our spatially varying DEB model. Section 3 presents the computational setting. Parameters are fitted to measured death rates in Section 4. Section 5 discusses numerical simulations illustrating the above mentioned qualitative behaviors. Finally, Section 6 states our conclusions.

2 Basic DEB model for a biofilm including EPS production and antibiotics

In this section we propose the basic equations for a dynamic energy budget theory (DEB) of a biofilm including polymer (EPS) production, coupled to diffusion of oxygen and antibiotic concentrations. We also introduce equations describing the effect of the antibiotic on the bacteria. The framework is similar to that developed in [21] for the study of the effect of toxicants on a homogeneous population of floating bacteria. Differently, in a biofilm we must incorporate EPS production, as well as diffusion processes in a number of spatially varying magnitudes, taking into account the specific action mechanisms of antibiotics. Also, the population of bacteria differentiates in several types. Only some of them become EPS producers. The remaining bacteria grow and divide normally, unless resources are so scarce to trigger deactivation or damage is large enough to kill them.

2.1 The basic energies

DEB is essentially a scaling theory of different types of energies characterizing the biofilm. The basic differential equations describe the fluxes of these energies. We apply a variant of the $\kappa$ rule [23], see Figure 2. A $\kappa$ fraction of the energy is used for bacterial growth and division. Notice that bacteria do not mature and age, just grow and divide. A $1-\kappa$ fraction of the energy is used for EPS production.

Each alive bacterium in the biofilm evolves according to the following system of equations, setting $\kappa=1$ for normal cells and $0<\kappa<1$ for EPS producers.

Scaled energy density $e(t)$:

[TABLE]

where $C_{o}$ is the oxygen limiting concentration, $K_{o}$ the limiting concentration half saturation value, $f$ is the scaled functional response, $\nu$ is the energy conductance and $\nu^{\prime}$ takes into account toxic effects on conductance through $\nu_{A}$, defined in equation (7). 2. 2.

Dimensionless cell volume $v(t)$:

[TABLE]

where $r$ is the bacterial (cell biomass) production rate, $m_{\kappa}$ is the maintenance rate and $g_{\kappa}$ is the investment ratio. The remaining parameters and magnitudes are linked to environmental toxicity, in our case antibiotic concentration, see below. $h$ is the hazard rate, $a$ is the acclimation energy density and $a_{M}$ is the target acclimation energy. The symbol + stands for positive part, which becomes zero for negative expressions.

Bacteria are rod-shaped, and grow in length. If we fix an average radius, equation (2) provides an equation for the time evolution of the cell length $l$. 3. 3.

Dimensionless volume of EPS $v_{e}(t)$:

[TABLE]

This equation defines a rate of EPS production $r_{e}$. In absence of acclimation and hazard effects ($a=a_{M}$, $h=0$), we set the reference rate $r_{e}=kr+k^{\prime},$ that can be fitted to experiments, as we explain in Section 5.

The dimensionless magnitudes $e,v$ and $v_{e}$ represent the fundamental energies characterizing the cell (bacterium) and its EPS production. Their fluxes are written here for each individual cell, $e$ being the total energy available for both growth and reproduction and for EPS production. $v$ keeps track of the volume of the cell as it grows, this is the energy available for reproduction by division. $v_{e}$ keeps track of the energy available to form the matrix that constitutes the biofilm and protects the bacterium. They may be related to dimensional volumes $V$ and $V_{e}$ multiplying by a characteristic cell volume $V_{m}$. The dynamics of these variables is coupled to the evolution of other magnitudes, governed by equations we introduce next. Table 1 collects all the variables appearing in the model. Table 2 gathers the parameters, specifying units and selected values. The discussion in Section 4 suggests that a characteristic time scale for equations (1)-(3) is set by the bacterial doubling time $\sim{1\over\mu_{max}}$ or the conductance $\nu$, in a time scale of hours. Diffusion processes for oxygen and antibiotic concentrations, however, evolve in a time scale of seconds, as concluded from the equations presented next and the parameter values in Table 2.

2.2 Diffusion of oxygen and antibiotics

Let $C_{o}(\mathbf{x},t)$ denote the concentration of oxygen and $C_{a}(\mathbf{x},t)$ the concentration of antibiotic. These quantities are described by two diffusion equations. The antibiotic accumulated inside the cell is described using the antibiotic cellular density $[C_{IN}](t)$, as discussed next.

Oxygen concentration diffusion inside the biofilm:

[TABLE]

where $\mu_{max}$ is the maximum specific growth rate. $k$ is the growth associated polymer formation rate coefficient and $k^{\prime}$ is the non-growth associated polymer formation rate coefficient. $Y_{x/o}$ is the cellular yield coefficient of $C_{o}$ and $Y_{p/o}$ is the polymer yield coefficient of $C_{o}$. $d_{o}$ is the oxygen diffusion coefficient inside the biofilm. $X$ is the cellular structure concentration computed in a control volume $V_{T}$ containing $N$ cells: $X=\rho_{x}\sum_{i=1}^{N}V_{i}/V_{T}$, $\rho_{x}$ being the density of a cell. When the volume control contains just one cell, we set $X=\rho_{x}v$, where $v$ is the dimensionless cell volume. Equation (4) is supplemented with the boundary condition $C_{o}=C_{o,out}$ at the interface with the oxygen providing fluid and the no-flux condition ${\partial C_{o}\over\partial\mathbf{n}}=0$ at the biofilm/substratum interface. 2. 2.

Antibiotic concentration diffusion inside the biofilm:

[TABLE]

where $R_{e}$ and $R$ are evaluated averaging $r$ and $r_{e}$ in the control volumen $V_{T}$. $d_{a}$ is the antibiotic diffusion coefficient. Equation (5) is supplemented with boundary condition $C_{a}=C_{a,out}$ at the interface with the antibiotic providing fluid and the no-flux condition ${\partial C_{a}\over\partial\mathbf{n}}=0$ at the biofilm/substratum interface. 3. 3.

Antibiotic cellular density:

[TABLE]

where $k_{A}^{I}$ is the antibiotic influx coefficient, and $k_{A}^{O}$ is the antibiotic efflux coefficient. $C_{a}$ should be the antibiotic concentration outside the cell, that we take equal to $C_{a}$ given by (5) in the control volume that contains the cell. We may set $[C_{IN}(0)]=0$ initially.

In view of the parameter values in Table 2, the chemical concentrations $C_{o}$ and $C_{a}$ will evolve much faster than the other magnitudes related to cell behavior, that is, in a time scale of seconds, not hours.

We have to complete the above equations by describing how the antibiotic and oxygen concentration act on the cell and influence its energy fluxes. We mostly follow [22, 21] considering the antibiotic to be a toxicant. Besides, we introduce several new DEB variables and make small modifications to account for spatial variations.

2.3 The effect on cells of oxygen and antibiotics

The previous equations take into account the decline in the organism’s capacity to acquire and use energy due to respiration and to the presence of antibiotics through a diminished conductance, reduced growth due to acclimation to the antibiotic and the hazard rate, see Figure 3.

Conductance modified by exposure $\nu^{\prime}$:

[TABLE]

where $K_{V}$ is the noncompetitive inhibition coefficient and $\gamma_{\varepsilon}$ is the environmental degradation effect coefficient. 2. 2.

Environmental degradation $\varepsilon(\mathbf{x},t)$:

[TABLE]

$\nu_{\varepsilon}$ is the environmental degradation coefficient and $\nu_{m}$ is the maintenance respiratory coefficient. $X$ is the cellular structure concentration computed in the control volume $V_{T}$ containing $N$ cells. $R$ is similarly computed averaging $r$ over this control volume. The diffusive term accounts for the fact that $\varepsilon$ feels the spatial variations. We impose no-flux boundary conditions. 3. 3.

Acclimation energy density $a(t)$:

[TABLE]

where $a_{M}$ is the target acclimation energy and + stands for the positive part. 4. 4.

Hazard rate $h(t)$:

[TABLE] 5. 5.

Aging acceleration $q(t)$:

[TABLE]

where $h_{a}$ is the Weibull aging acceleration and $s_{q}$ is a multiplicative stress coefficient. 6. 6.

Aging in acclimation due to dissolved antibiotic and EPS:

[TABLE]

where $k_{qA}^{I}$ is the dissolved antibiotic toxicity.

2.4 Cell division, differentiation and death

Bacteria of the Pseudomonas Aeruginosa genus are rod-like. Their volume is defined by $V=\pi\rho^{2}L$, $\rho$ being the radius and $L$ the length. Increase in volume is equivalent to increase in length, setting $\rho$ equal to their average radius. Reference [34] reports cell lengths $L$ between $1.67$ $\mu$m and $2.75$ $\mu$m for P. Aeruginosa. Diameters $2\rho$ vary between $0.5$-$0.73$ $\mu$m without clear relation to the length and with a preference for $2\rho=0.6$ $\mu$m. Notice that the equation for length is formally the same as the equation for volume (2), up to a factor. We may use $V_{m}=\pi\rho^{2}L_{m}$, $L_{m}$ being the maximum length, to nondimensionalize.

In dimensionless variables, cell division happens provided a threshold maintenance energy $e_{d}$ is surpassed when a threshold length $l_{d}$ is reached, see Figure 4. Then the cell divides and creates daughter cells.

Nondividing cells behave in different ways according to their lengths, energies, aging and hazard rates, see Figure 4:

For intermediate energies $e_{eps1}\leq e\leq e_{eps2}$ the cells produce EPS according to equation (3), provided they are not newly born. A way to keep track of their age is the aging rate $q$, governed by (11). It is set initially equal to zero for newborn cells. We require $q\geq q_{eps}$ for the cell to have the ability to become a EPS producer.

• 2.

For low energies $e\leq e_{eps1}$, the cell just maintains itself.

• 3.

The hazard rate, given by (10), determines when cells die. Their survival probability $p(t)$ is governed by the equation

[TABLE]

Now with the energies $e,v,v_{e}$, oxygen and antibiotic concentrations $C_{o}$,$C_{a}$, $[C_{IN}]$ and the variables $\epsilon,a,h,q$, which measure how such concentrations influence the energies, we have a complete set of variables and of equations. These additional rules define the behavior of each cell.

3 Computational framework

For computational purposes, the biofilm is identified with a slab containing cells. We divide the slab in a grid of cubic control volumes, that we use to discretize the concentration fields evaluating their spatial changes. A control volume $V_{T}$ may contain several cells, or just one, see Figure 5.

Variables $e$, $v$, $v_{e}$, $[C_{IN}]$, $a$, $h$, $q$, $p$ are assigned to each cell, governed by equations (1), (2), (3), (6), (9), (10), (11) and (13). These equations use background values of concentrations and environmental degradation, that are continuous fields defined on the biofilm, governed by the diffusion problems (4), (5) and (8).

The computational strategy is the following:

For a fixed biofilm configuration, we let the diffusion problems for $C_{o}$, $C_{a}$ and $\varepsilon$ relax to stationary values, which occurs in a short time scale $\tau$ (seconds). These variables are then defined and fixed in all the control volumes.

• 2.

For each cell in the biofilm, we update the values of the variables $e$, $v$, $v_{e}$, $[C_{IN}]$, $a$, $h$, $q$, $p$ solving equations (1), (2), (3), (6), (9), (10), (11) and (13) in a slower time scale $T$ (minutes) using the background values of $C_{o}$, $C_{a}$ and $\varepsilon$.

• 3.

We revise the status of all cells in the biofilm, in a random ordering. Cells differentiate and divide according to the rules stated in Section 2.4. To decide whether a cell dies, we generate a random number $n\in(0,1)$ and kill the cell if $p<n{N_{a}\over N},$ where $N$ is the total number of cells and $N_{a}$ the current number of alive cells.

• 4.

We update the concentration fields and repeat the procedure.

Allocating newborn cells, reabsorbing dead cells and allocating produced EPS in a 2D or 3D biofilm geometry are challenging issues for which many ideas have been proposed depending on the biofilm types [10, 11, 12, 24, 25, 39]. We do not intend to address these issues here. Instead, we work in a simplified geometry, implementing simple rules that allow us to check the performance of the DEB dynamics in a spatially varying film. In the numerical simulations presented here we make the following assumptions:

The biofilm is a 2D region, divided in cubic tiles for computational purposes.

• 2.

Each tile contains one cell together with EPS and dissolved substances. The idea is schematically represented in Figure 5.

• 3.

When a cell divides, the length and energy split in two. The daughter cell shifts other cells in the direction of minimal mechanical resistance (closest to the surface or to a dead cell).

• 4.

When a cell dies, it is reabsorbed by its neighbors, provided there are enough neighbors alive. A neighbor or newborn cell occupies its place in that case. Otherwise, its stays there and we get a necrotic region.

• 5.

A fraction $\alpha\in(0,1)$ of the EPS produced by equation (3) remains attached to the cell. The rest diffuses through the whole biofilm. Therefore, when updating the concentrations $C_{a}$, $C_{o}$, $\varepsilon$, we need to update the diffused concentration of EPS $C_{e}(\mathbf{x},t)$ too:

[TABLE]

with zero flux boundary conditions ${\partial C_{e}\over\partial{\bf n}}=0$ on the biofilm boundary. The equations in Sections 2.2-2.3 set $\alpha=1$. All the EPS produced by a cell remains attached to it. When $\alpha\neq 1$, $R_{e}$ is replaced by $\alpha R_{e}+R_{e}^{\prime}$ in all these equations. $\alpha R_{e}$ represents the fraction of EPS produced at each site that remains in it. It may be zero when there are no EPS producers in that tile. $R_{e}^{\prime}$ represents the fraction of EPS accumulated due to global diffusion processes governed by (14). It will always be non zero by diffusion. The way we evaluate $R_{e}^{\prime}$ is the following. If at a certain control volume we start with an initial concentration $C_{e}(t_{1})$ to reach a final concentration $C_{e}(t_{2})$ this defines a rate $R_{e}^{\prime}={1\over t_{2}-t_{1}}\ln({C_{e}(t_{2})\over C_{e}(t_{1})}).$ The overall computational chart is represented in Figure 6.

The physical justification for this treatment of EPS is that the EPS matrix is formed by polymers. Monomers will diffuse easily, but a fraction of them will form polymeric chains of increasing size attached to the cells. An empirical justification is that if the produced EPS remains attached to producers, the active cells not producing EPS in outer layers die too fast and form unphysical necrotic layers in absence of any toxicants. In practice, the EPS matrix envelops and shelters all the biofilm cells.

4 Parameter calibration

Table 2 lists the parameters used in the simulations. The diffusion coefficient for oxygen $d_{o}$ is selected so that oxygen penetration inside the biofilm ranges in the values reported (about $50$ $\mu m$) [42]. Oxygen penetration depth is defined as the distance into the biofilm at which the first derivative of the oxygen concentration reaches $5\%$ of its maximum value. As said before, oxygen is considered to be the limiting factor [16] for biofilm development in medical environments. We lack measurements of $k$, $k^{\prime}$, $Y_{x/o}$ and $Y_{p/o}$ in this case, but we may propose values for $Y_{x/o}$ and $Y_{p/o}$, $k$ and $k^{\prime}$ based on an educated guess respecting the proportions observed for carbon [34].

Values for the maintenance rates $m_{\kappa}$ and investment ratios $g_{\kappa}$ are inferred as follows. The macroscopic equations for EPS and cell production are [34]:

[TABLE]

In absence of acclimation, toxicity and harzard corrections, that is $a=a_{M}$, $\nu^{\prime}=\nu$, $h=0$, equations (2) and (3) read

[TABLE]

Notice that the energy governed by (1) tends to an equilibrium $e=f$ for each cell. We set $f=f_{max}={C_{o,out}\over C_{o,out}+K_{o}}$ and $e=f_{max}$. Dimensions can be restored in these equations by scaling $v$ and $v_{e}$. Comparing with the equations for $X_{e}$ and $X$, we approximate ${g_{\kappa}\over m_{\kappa}}$ by ${k\over k^{\prime}}$ and ${\nu f_{max}-m_{\kappa}g_{\kappa}\over f_{max}+g_{\kappa}}$ by $\mu_{\max}f_{max}$. Taking $m_{\kappa}={k^{\prime}\over k}g_{\kappa}$, $g_{\kappa}$ is a positive solution of a second order equation, given by

[TABLE]

The values of the multiplicative stress coefficient $s_{q}$, Weibull aging acceleration $h_{a}$, dissolved antibiotic toxicity $k_{q,A}^{I}$, influx and efflux coefficients $k_{A}^{I}$, $k_{A}^{O}$, and the initial hazard rate $h_{0}$ in equations (10), (11), (12) and (13) are fitted to data in reference [17]. To do so, we consider a single cell of volume $v$, energy $e$ with density $\rho_{x}$ and write down the equations for the hazard rate $h$, aging acceleration $q$ and survival probability $p$ setting $[C_{IN}]=C_{a,out}{k_{A}^{I}\over k_{A}^{O}},$ the equilibrium value. Taking $a=h_{0}$, $b=r+r_{e}$, $c=e\,h_{a}(\nu^{\prime}-r)$, $d=e\,s_{q}\,\rho_{x}v(\nu^{\prime}-r)$ and $s=k_{qA}^{I}[C_{IN}]$, the solutions are:

[TABLE]

[TABLE]

[TABLE]

The way we find values for $a$, $b$, $c$, $d$, $s$ is the following. We first remove the antibiotic and seek to adjust the viable cell counts in Table 1 of reference [17] for the control population. Figure 7 compares the bacterial counts at times $3$, $6$, $24$ hours with our fittings of $log_{10}(pN)$ for a population of the initial size $N=10^{5.75}$ CFU/m. Once these parameters are fitted, we successfully fit the remaining parameters to the approximate number of viable cells for four different antibiotics: levofloxacin, ciprofloxacin, gentamicin, ceftazidime.

We approximate the data in Table 1 of [17] in absence of antibiotics dividing by $10$ the values $s_{q}$ and $h_{a}$ from reference [21]. We choose average values for $e$ and $v$, $v=0.5$, $e=median(f(C_{o,out}))\approx 0.003$, and set $\nu^{\prime}=\nu$. Then, $c=e\,h_{a}(\nu^{\prime}-r)\approx 3.61\times 10^{-7}$ y $d=e\,s_{q}\,\rho_{x}v(\nu^{\prime}-r)\approx 4.8\times 10^{-3}$. We vary $a$ and $b$ so as to approximate the data in Table 1 of [17]: $a=h_{0}=0.4$ and $b=0.29$. These values are chosen because the fit still holds when we replace single cells by a biofilm formed by a few hundred layers. Since $r\sim 0$ for EPS cells, we set $k^{\prime}=b$ and keep the ratio ${k\over k^{\prime}}$ observed for carbon.

5 Numerical results

The purpose of this section is to investigate the influence of the spatial variations of the concentrations in cell death. To simplify, in the small scale 2D simulations presented here we consider tiles of size $1\,\mu m\times 1\,\mu m$ in Figure 5. Figure 8 shows the spatial structure of the oxygen distribution in a biofilm of maximum width $1000\,\mu m$ and maximum height $200\,\mu m.$ We solve equation (4) discretized in the grid, with Dirichlet boundary condition $C_{o}=C_{o,out}$ on the interface with air and Neumann boundary condition ${dC_{o}\over dn}=0$ on the bottom interface. We take $f(C_{o})={C_{o}\over C_{o}+K_{o}}$ and $X=\rho_{x}v$, where $v$ is the dimensionless volume of the cell contained in the tile. Figure 8 sets $v=0.5$. Initially, $C_{o}=0$ inside the biofilm. We solve using a explicit scheme until $C_{o}$ relaxes to a stationary configuration:

[TABLE]

where $g(X)$ represents the source in equation (4). The spatial step $\delta x$ is the side of a tile. Tiles in the spatial grid are labelled using a couple of indices $(i,j)$ to indicate their position. Time is discretized with time step $\delta t$, so that $t_{\ell+1}=t_{\ell}+\delta t=(\ell+1)\delta t$ and $C_{o;i,j}^{\ell}\sim C_{o;i,j}(t_{\ell})$. For this scheme to be stable, we need ${\delta t\,d_{o}\over\delta x^{2}}\leq 0.5$. We advance in a time scale of seconds. Similar schemes are used to update the additional space dependent concentrations $C_{a}$, $C_{\varepsilon}$, $C_{e}$ when needed.

The cell energies are governed by DEB equations depending on $f={C_{o}\over C_{o}+K_{o}}$, which varies through the biofilm. Figure 9 represents the equilibrium distribution $e=f(C_{o})$ at each tile. The maximum energy is $e_{m}={C_{o,out}\over C_{o,out}+K_{o}}.$ In view of the energy distribution, we set threshold values $e_{eps1}=7\times 10^{-6}$, $e_{eps2}=0.02$ and $e_{d}=0.01$ to produce EPS or divide.

Starting from a uniform distribution of cell types (non EPS producers, $v_{e}=0$) with volumes $v(0)=0.5$ we compute the evolution of their energy and size to check the evolution of the cell type distribution with time. We choose as initial energy a perturbation of the equilibrium energy distribution at each location, and fix $q>q_{eps}=10^{-8}.$ By solving the equations for $e$, $v$, $v_{e}$ for each cell tile ignoring aging and hazard, we see that cells evolve to be EPS producers in an intermediate region. The outer layers are normal active cells, likely to divide. A thin inner core may be formed by cells with little energy than do not divide and do not produce EPS. The equilibrium distribution of types is depicted in Figure 10(a). If we increase the initial energy, we may find intermediate states such those in Figure 10(b), where that thin inner layer vanishes. Repeating the simulation with initial volumes $v$ randomly distributed between $(0.45,0.6)$ and random perturbations of the equilibrium energy we find similar results with a smoother transition between regions, reflecting the noise.

To study the distribution and number of dead cells, as well as the effect of antibiotics, we activate the equations for aging and hazard rates. We compute a reference concentration of antibiotic $C_{a}$ in the biofilm solving (5) with Dirichlet boundary condition $C_{a}=C_{a,out}$ on the interface biofilm/air and Neumann boundary condition ${dC_{a}\over dn}=0$ on the bottom interface. We assume that all cells are undifferentiated in this reference computation ($R_{e}=0$) and employ the equilibrium oxygen and energy distributions to evaluate $R$. The concentration of antibiotics inside the biofilm evolved towards a constant, and we will take it to be approximately constant and equal to $C_{a,out}$ in the biofilm. Figure 11 shows that the fittings of the time evolution of the survival rates performed in Section 3 persist when all the biofilm is taken into account. A low energy cell population is considered in that figure, forbidding cell division for simplification.

The model yields information on qualitative behavior that agrees with some experimentally observed tendencies and may be useful to understand the causes or to predict new behaviors. Panels (a) and (c) in Figure 12 reveal that for an antibiotic acting mostly on active cells like Ceftazidime (it inhibits cell wall synthesis), cells mostly die in the outer biofilm layer and the extent of the dead region increases with the dose of antibiotic. Panel (b) shows that it also increases with time. We are assuming that the antibiotic has diffused inside the biofilm so that the spatial distribution of dead regions depends on what type of cells are preferential targets of the antibiotic under study. For antibiotics like Ceftazidime, acting on active cells, we use toxicity coefficients depending on the cell activity, of the form ${e\over 0.1\,{\rm median}(e)}k_{qA}^{I}$. For other antibiotic like colistin, targeting the dormant core cells, we may use toxicity coefficients decreasing with the cell activity instead. Mutations increasing the efflux might be accounted for raising the value of that coefficient. Augmented resistance to deactivating enzymes may be incorporated modifying the toxicity.

In absence of antibiotics cells die in the inner regions in a scattered way. Active cells in the outer layers are still more likely to die, a standard fact in the DEB theory, that assumes organisms with more energy more likely to suffer damage. That effect would be diminished increasing the difussivity of produced EPS $d_{e}$, or allowing differentiation into EPS producers for larger energies. Cells that die in the outer layers in absence of antibiotics may also be eroded by external agents or absorbed by newly born cells. Unlike the previous case, we do not get an expanding outer necrotic region. Analyzing the amount of EPS produced with and without antibiotics, we find that the presence of antibiotics enhances EPS production, see Figure 13. This fact is also responsible for the small difference between panels (a) and (c) of Figure 12. The number of dead cells is larger for panel (c), but not as large as could be expected from the dose increase due to enhanced EPS production that neutralizes that fact.

6 Conclusions

Coupling dynamic energy budget dynamics for the metabolism of individual cells to diffusion equations for the evolution of chemical fields we have reproduced qualitative trends experimentally observed in bacterial biofilms under the action of antibiotics. We calibrate the parameters for a quantitative agreement with the death rates measured for P. aureginosa under a number of antibiotics and study spatial variations in the distribution of dead and alive cells. Our numerical simulations show that the presence of antibiotics enhances production of polymeric matrix, as expected from true measurements. For antibiotics targeting active cells, we observe the formation of a necrotic outer layer, that expands deeper in the biofilm as time goes on and the antibiotic dose increases. Our tests suggest ways to handle antibiotics targeting dormant cells and resistance mechanisms such as efflux pumps or enzymes. New experimental measurements [31] would be required to calibrate the mechanisms to practical cases. The current computational model enables validation and hypothesis testing when developing therapies to handle chronic infections caused by biofilms. Determining the true relevance of increasing the dose and whether periodic or continuous infusion of antibiotics is more effective [20] would be essential for practical applications.

Acknowledgements. Research partially supported by the NILS Mobility Program and MINECO grant No. MTM2014-56948-C2-1-P.