Dynamic energy budget approach to evaluate antibiotic effects on biofilms

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.


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]. β-lactams (penicillins, cephalosporins, 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 β-lactamase, an enzyme that attacks β-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 β-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 β-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 β-lactam antibiotics into the bacteria provided the amount of chromosomal β-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.

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.

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 κ rule [23], see Figure 2. A κ 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 − κ fraction of the energy is used for EPS production.
Each alive bacterium in the biofilm evolves according to the following system of equations, setting κ = 1 for normal cells and 0 < κ < 1 for EPS producers.
1. Scaled energy density e(t): where C o is the oxygen limiting concentration, K o the limiting concentration half saturation value, f is the scaled functional response, ν is the energy conductance and ν takes into account toxic effects on conductance through ν A , defined in equation (7).
where r is the bacterial (cell biomass) production rate, m κ is the maintenance rate and g κ 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.

Dimensionless volume of EPS
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 , 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 ∼ 1 µmax or the conductance ν, 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.

Diffusion of oxygen and antibiotics
Let C o (x, t) denote the concentration of oxygen and C a (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.
1. Oxygen concentration diffusion inside the biofilm: where µ max is the maximum specific growth rate. k is the growth associated polymer formation rate coefficient and k is the non-growth associated polymer formation rate coefficient.   2. Antibiotic concentration diffusion inside the biofilm: 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 ∂Ca ∂n = 0 at the biofilm/substratum interface. 3. Antibiotic cellular density: where k I A is the antibiotic influx coefficient, and k O A 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.

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 where K V is the noncompetitive inhibition coefficient and γ ε is the environmental degradation effect coefficient. 2. Environmental degradation ε(x, t): ν ε is the environmental degradation coefficient and ν 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 ε feels the spatial variations. We impose no-flux boundary conditions. 3. Acclimation energy density a(t): where a M is the target acclimation energy and + stands for the positive part. 4. Hazard rate h(t): 5. Aging acceleration q(t): where h a is the Weibull aging acceleration and s q is a multiplicative stress coefficient. 6. Aging in acclimation due to dissolved antibiotic and EPS: where k I qA is the dissolved antibiotic toxicity.

Cell division, differentiation and death
Bacteria of the Pseudomonas Aeruginosa genus are rod-like. Their volume is defined by V = πρ 2 L, ρ being the radius and L the length. Increase in volume is equivalent to increase in length, setting ρ equal to their average radius. Reference [34] reports cell lengths L between 1.67 µm and 2.75 µm for P. Aeruginosa. Diameters 2ρ vary between 0.5-0.73 µm without clear relation to the length and with a preference for 2ρ = 0.6 µ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 = πρ 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 ≤ e ≤ 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 ≥ q eps for the cell to have the ability to become a EPS producer.
• For low energies e ≤ e eps1 , the cell just maintains itself.
• The hazard rate, given by (10), determines when cells die. Their survival probability p(t) is governed by the equation Now with the energies e, v, v e , oxygen and antibiotic concentrations C o ,C a , [C IN ] and the variables , 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.

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.
The computational strategy is the following: • For a fixed biofilm configuration, we let the diffusion problems for C o , C a and ε relax to stationary values, which occurs in a short time scale τ Each tile contains cells of certain volume and energy, with an associated EPS volume around it, whose dynamics is governed by the DEB equations. Here, we assume that each tile contains just one cell. In practice, all these items are stuck together in a more amorphous mass. This is a simple computational idealization allowing us to combine macroscopic and microscopic information. • 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 ∈ (0, 1) and kill the cell if p < n Na N , where N is the total number of cells and N a the current number of alive cells.
• 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 Each tile contains cells of certain volume and energy, with an associated EPS volume around it, whose dynamics is governed by the DEB equations. Here, we assume that each tile contains just one cell. In practice, all these items are stuck together in a more amorphous mass. This is a computational idealization allowing us to combine macroscopic and microscopic information.
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.
• Each tile contains one cell together with EPS and dissolved substances. The idea is schematically represented in Figure 5.
• 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).
• 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.
• A fraction α ∈ (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 , ε, we need to update the diffused concentration of EPS C e (x, t) too: with zero flux boundary conditions ∂Ce ∂n = 0 on the biofilm boundary. The equations in Sections 2.2-2.3 set α = 1. All the EPS produced by a cell remains attached to it. When α = 1, R e is replaced by αR e + R e in all these equations. α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 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 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 = 1 t 2 −t 1 ln( Ce(t 2 ) Ce(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.  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 , 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 based on an educated guess respecting the proportions observed for carbon [34].

Parameter calibration
Values for the maintenance rates m κ and investment ratios g κ are inferred as follows. The macroscopic equations for EPS and cell production are [34]: In absence of acclimation, toxicity and harzard corrections, that is a = a M , ν = ν, h = 0, equations (2) and (3) read Notice that the energy governed by (1) tends to an equilibrium e = f for each cell. We set f = f max = Co,out Co,out+Ko 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κ mκ by k k and νfmax−mκgκ fmax+gκ by µ max f max . Taking m κ = k k g κ , g κ is a positive solution of a second order equation, given by The values of the multiplicative stress coefficient s q , Weibull aging acceleration h a , dissolved antibiotic toxicity k I q,A , influx and efflux coefficients k I A , k O A , 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 ρ x and write down the equations for the hazard rate h, aging acceleration q and survival probability p set- .
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 )) ≈ 0.003, and set ν = ν. Then, c = e h a (ν − r) ≈ 3.61 × 10 −7 y d = e s q ρ x v(ν − r) ≈ 4.8 × 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 ∼ 0 for EPS cells, we set k = b and keep the ratio k k observed for carbon.

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 µm × 1 µm in Figure  5. Figure 8 shows the spatial structure of the oxygen distribution in a biofilm of maximum width 1000 µm and maximum height 200 µ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 dCo dn = 0 on the bottom interface. We take f (C o ) = Co Co+Ko and X = ρ 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 Figure 7: Comparison of theoretical predictions and experimental data. The predictions provided by formula (15) for p are compared to the data in Table 1 of Reference [17]. The parameter values for p(t) are a = 0.4, b = 0.29, c = 3.61 × 10 −7 and d = 4.8 × 10 −3 . The additional parameter s is adjusted for each antibiotic. Circles represent log 10 (p · N ) with N = 10 5.75 (magenta) and N = 10 5.75±0.10 (red). Asterisks represent the measurements in Table 1 of Reference [17] for the same antibiotic doses. As before, light (blue) colors correspond to the mean measured value, whereas dark (black) colors mark the error interval.  relaxes to a stationary configuration: where g(X) represents the source in equation (4). The spatial step δ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 δt, so that t +1 = t + δt = ( + 1)δt and C o;i,j ∼ C o;i,j (t ). For this scheme to be stable, we need δt do δx 2 ≤ 0.5. We advance in a time scale of seconds. Similar schemes are used to update the additional space dependent concentrations C a , C ε , C e when needed.
The cell energies are governed by DEB equations depending on f = Co Co+Ko , which varies through the biofilm. Figure 9 represents the equilibrium distribution e = f (C o ) at each tile. The maximum energy is e m = Co,out Co,out+Ko . In view of the energy distribution, we set threshold values e eps1 = 7 × 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 dCa 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 Figure 11: Simulated live cell counts compared to experimental measurements in Table 1 of Reference [17]: (a) in absence of antibiotics, (b) for ceftazidime. The continuous curve represents the simulated values. Asterisks depict the measurements in Table 1 of reference [17] for the same antibiotic doses, black lines represent the error bars.
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 0.1 median(e) k I qA . 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.

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 (a) (b) 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.