Computational and Mathematical Methods in Medicine

Volume 2017, Article ID 7246286, 11 pages

https://doi.org/10.1155/2017/7246286

## Mathematical Modeling of Biofilm Structures Using COMSTAT Data

^{1}Department of Clinical Pharmacy, School of Pharmacy, University of California San Francisco, San Francisco, CA, USA^{2}Novo Nordisk Foundation Center for Biosustainability, Technical University of Denmark, Kogle Alle 6, 2970 Hørsholm, Denmark^{3}Department of Civil and Environmental Engineering, James H. Clark Center, Stanford University, Rm E250, 318 Campus Drive, Stanford, CA 94305, USA

Correspondence should be addressed to Davide Verotta; ude.fscu@attorev.edivad

Received 10 August 2017; Revised 14 November 2017; Accepted 26 November 2017; Published 20 December 2017

Academic Editor: Michele Migliore

Copyright © 2017 Davide Verotta et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

Mathematical modeling holds great potential for quantitatively describing biofilm growth in presence or absence of chemical agents used to limit or promote biofilm growth. In this paper, we describe a general mathematical/statistical framework that allows for the characterization of complex data in terms of few parameters and the capability to (i) compare different experiments and exposures to different agents, (ii) test different hypotheses regarding biofilm growth and interaction with different agents, and (iii) simulate arbitrary administrations of agents. The mathematical framework is divided to submodels characterizing biofilm, including new models characterizing live biofilm growth and dead cell accumulation; the interaction with agents inhibiting or stimulating growth; the kinetics of the agents. The statistical framework can take into account measurement and interexperiment variation. We demonstrate the application of (some of) the models using confocal microscopy data obtained using the computer program COMSTAT.

#### 1. Introduction

Biofilms are structured communities of bacteria enclosed in an extracellular matrix composed of polysaccharides, proteins, and extracellular DNA adherent to a surface [1]. Unlike planktonic bacteria, biofilms exhibit differences in metabolism, antibiotic tolerance, and ability to evade the immune system, making infections due to biofilms difficult to treat [2]. Biofilms are a main cause of acute and chronic infections, including foreign-body infections, otitis media, and urinary tract infections.

When a population of microorganisms organized in a biofilm grows, it is likely to pass over different phases. The growth may follow a period of dormancy, if the environmental conditions before the beginning of the growth are not optimal. Eventually, cells start to divide and structure, and the biofilm grows into a period in which the overall rate of cell division prevails over that of their death. Under favorable conditions, the growth may be considered to be unlimited (hence exponential) for some time, but eventually [3] physiological and physical limits such as (i) exhaustion of available nutrients, (ii) accumulation of inhibitory metabolites or end products, and (iii) exhaustion of space intervene. As a result, the growth rate decreases and the colony reaches its maximum size. Repeated or continual exposure to these environmental or physiological stressors can then result in a decline in the biofilm size [4].

In vitro biofilms are often used in studies concerning therapy, dealing with the reactions of bacterial populations to various agents: drugs, for example, antibiotics, or mutagens [5–8]. Those are applied externally to the biofilm structures and change their environment or directly eliminate (kill) the bacteria or decrease their reproductive capacity.

A variety of mathematical models have been used to describe bacterial growth in investigation of dynamics in environment depending only on the activity of the bacteria [9, 10]. Similarly, a number of models have been proposed to describe the action of different agents, in particular drugs and their interaction [11–13]. The objective of this paper is to obtain an integrated framework that provides models describing the different stages of biofilm growth, the action of different agents, and the simultaneous modeling of resulting kinetics for live and dead biofilm. The statistical issues associated with the use of these data, in particular the treatment of different sources of variability, are also addressed.

We demonstrate the use of the resulting general models using data obtained by means of confocal microscopy and COMSTAT [14]. COMSTAT takes the image stacks created by the confocal microscope as source data and produces up to ten image analysis features for quantification of biofilm structures which are output as one or more text files. The models we describe in this paper apply to univariate measurements: total biomass, area in a specific layer, average thickness, and volumes of microcolonies identified at the substratum (COMSTAT also obtains multivariate data, such as thickness distribution, which can be used to quantify the three-dimensional structures in the biofilm. The modeling of such data is the subject of current modeling investigation and will be reported in future communications.).

#### 2. Methods

##### 2.1. Mathematical Modeling

Modeling of biofilm growth requires the specification of three components. The first, a function , describes the growth of the biofilm in absence of agents limiting or promoting growth. The second, a nonnegative function , describes the interaction with the agents. The third, , describes the temporal variation (kinetics) of the agents acting on the biofilm. The general model we consider expresses the rate of change of the biofilm as follows:where is time, is biofilm amount present in the system at time , is the initial conditions or biofilm amount at , and is the concentration of agent at time . The ± sign indicates that the model can describe either inhibition or stimulation.

###### 2.1.1. Biofilm Growth

A large number of models have been proposed to describe bacterial growth, for example, [15–17]. For the purpose of this paper, we only describe the simplest model for unlimited growth (exponential) and three semiempirical models for limited growth. Exponential biofilm growth assumes that the rate of growth is proportional to the amount of cells present in the system, and that there is no limitation to growth; thuswhere is the growth rate of the biofilm. The analytic solution of (2) is an exponential growth:with doubling time, , being equal to .

In general, exponential growth can only describe the early stages of biofilm growth. A number of semiempirical models can be used to take into account the decrease in proliferation that results from limitations of nutrient supply or mechanical constraints or metabolites accumulation. We consider the Logistic, Gompertz, and Bertalanffy models. The Logistic [18] takes the following form:where is the maximum biofilm level that can be reached*. *The corresponding analytic solution is as follows:The Gompertz model [19–21] takes the following form:with analytic solution:Finally the Bertalanffy model [22] assumes that growth occurs proportionally to biofilm surface area, while biofilm loss is proportional to biofilm:where is the death rate for biofilm. The solution of (8) is also sigmoidal in shape and tends to an asymptote as time increases, where the birth and death term balance each other out. A general version of the Bertalanffy model takes the following form [22]:where . The Logistic equation is a special case of this model with .

###### 2.1.2. Agents Interaction with Biofilm

The simplest model describing the interaction of bacterial biofilm with an agent assumes that the action of the agent is proportional to the product agent and biofilm:where now quantifies the interaction. For an inhibitory agent combining (10) with exponential growth (2) obtainsAccording to this model if the agent concentration is kept at constant level , the growth rate is zero. That concentration is the minimal concentration of agent killing bacteria or BIC. The BIC is frequently used to characterize bacterial growth in planktonic, in vitro experiments, and clinical settings (e.g., [23]). Note that if the growth rate is not exponential, the BIC is not a constant but depends on the current amount of biofilm. For example, for the Logistic model, incorporating an inhibitory agent action takes the following form:and the BIC is given by the following equation:that shows how the BIC decreases as a function of . Similarly, for the Gompertz model,More complex models for the interaction agent/biofilm can be used, for example, threshold models, according to which the agent is effective if its concentration reaches a threshold,where the parameter is the threshold; or models following saturation kinetics similar to what is used in, for example, pharmacodynamics [13],where now, for large concentrations of agent, the killing rate asymptotes to . In presence of two agents, and , models (10) and (16) increase in complexity, as the agents can affect the killing rate in different ways. For example, the linear model (10) can be generalized as follows:orwhere is a killing rate for the agents’ interaction. When saturation kinetics are present, as is the case for the single agent model (16), a number of possibilities arise, generating additive, synergistic, and antagonistic models often applied in pharmacodynamic, enzymology, or binding experiments [12]. For example, incorporating a model for competitive interaction between two agents yieldswhere if the model shows additivity, and when it reduces to the competitive antagonism model and yields

###### 2.1.3. Agent Kinetics

A mathematical representation of the agents’ kinetics is needed to incorporate their action on the biofilm. This in general does not present particular difficulties. Analytical solutions or sets of differential equations can be used to do so [24, 25]. As an alternative one can use “model independent” representations of the data, such as smoothing splines [26], or, when the measurement error in the kinetics data is low an even simpler linear interpolant as we do in the examples reported below.

###### 2.1.4. Modeling Dead Biofilm

Confocal microscopy data allow the measurement of both live and dead cells (see Data section). Indicating by the dead biofilm present in the system at time , its rate of generation can be directly obtained from the growth equations reported above by identifying the loss or gain in live cells due to the presence of the agent. For the Logistic growth rate equation (4), this yieldsand, for the generalized Bertalanffy model equation (8), an inhibitory agent obtainsFor the Gompertz model, the growth equation does not lead directly to express a rate of cell loss since it only expresses a decrease of the growth rate inversely proportional to biofilm growth. A possibility is to express dead biofilm as follows:the difference between what would result from exponential growth and the actual biofilm level.

###### 2.1.5. Modeling Post-Plateau Biofilm Decrease

To account for effects leading to a post-plateau decline in the biofilm size [4], one can introduce a hypothetical endogenous variable, , that is reduced proportionally to the amount of biofilm present [27, 28].where is the rate of reduction. is arbitrarily set to 1 at time zero; it is positive as increases, but eventually it becomes negative (as it can be seen by the relationship ). The growth rate of biofilm takes the following form:where ( is given by (4) or (6). For example, for Logistic growth,The right side of equations (25) and (26) becomes negative as time progresses and leads to asymptote to zero ( asymptotes to a negative value). The main limitation of the model is that for any value of the biofilm always asymptotes to 0. To avoid this the model can be modified as follows:where and are production and elimination rates of the endogenous substance. Equations for dead cells are obtained as before.

We remark that especially in in vivo systems one can observe pre- or post-plateau appearance of cyclical growth, associated with seeding and dispersal [29, 30]. To model these situations one could use a time-varying that changes value after a pure time-delay or depending on a threshold biofilm value.

###### 2.1.6. Modeling Dormancy

Biofilm growth can follow an initial period of dormancy that can be expressed easily using a (pure) lag time. The growth rate becomeswhere is the lag time.

##### 2.2. Statistical Modeling and Model Selection

Biofilm data usually consist of several measurements made on a number of growing “cells” on different occasions. The measurements present different levels of random variation: among measurements within a given cell (intracell variation), among cells (intercell variation), and interoccasion. A hierarchical nonlinear mixed effect model [31] can be used to represent such data. According to this model the th observation from the th cell and th experiment, , takes the following form:where is the corresponding sampling time, is a nonlinear function describing the relationship between time, parameters, and predicted response, and describes the parameters values for each cell, which depend on , a fixed effect mean parameter vector, a vector of intercell random effect parameters, , and a vector of intercell random effect parameters, , with mean zero and a certain, typically multivariate normal, distribution with variance-covariance matrix . The random effect indicates intracell variability (e.g., measurement error) typically assumed (multivariate) normally distributed with mean zero and covariance matrix , which could depend on (we omit details for simplicity) and an unknown parameter vector . Popular methods used to fit mixed effect models to data are implemented in the computer programs such as NONMEM [32] and NLME [33], which allow the computation of corresponding second-order statistics for the estimates, in particular the large sample variance-covariance matrix of the parameters , individual parameters as empirical-Bayes estimates, and corresponding predictions, conditional on the estimates for .

To select between different models one can use statistical model selection criteria together with the usual graphical displays based on predictions, observations, and residuals. Some statistical model selection criteria are the Hannan-Quinn (HQ) [34] and Akaike (AKA) [35]; these criteria penalize the likelihood of the model proportionally to the number of parameters in the model, thus penalizing the selection of models with ‘too many’ parameters. The HQ is the most conservative; it uses twice the log of the number of observations times the number of parameters as penalty, while the AKA uses twice the number of parameters.

##### 2.3. Data

*Pseudomonas aeruginosa* PAO1 tagged with green fluorescent protein (GFP) were studied in flow-chamber experiments described in [36]. In brief, biofilms were grown at 30°C in flow chambers. Each flow chamber was inoculated with 250 *μ*L of an overnight culture of PAO1 diluted to an OD_{600} of 0.05 and left without flow for one hour. After one hour, flow was started with minimal media at a flow rate of 20 ml/h using a peristaltic pump (Watson Marlow 205 S). After cultivation for 24 or 72 hours, flow was stopped and minimal media were replaced with an antibiotic flask containing the desired concentration of either Meropenem (MEM) or Tobramycin (TOB) where administered as an intermittent bolus. This antibiotic flask is connected to bubble traps and the flow chambers containing the cultivated PAO1 biofilms. Flow was restarted and minimal media were pumped from the dilution flask to the antibiotic flask to the flow chambers at a constant rate calculated to mimic the elimination rate constant of the antibiotic. MEM and TOB were obtained from the pharmacy of the University of California San Francisco Medical Center [37]. Concentration-time profiles were based on previously described PK parameters of MEM and TOB from healthy volunteers and patients with cystic fibrosis [38, 39]. The target MEM peak concentration based on human population values was computed to be 107.53 mg/L with a h. The target TOB peak concentration, based on a dose of 10 mg/kg in a 70 kg adult, was 32.79 mg/L with an associated h. Samples were analyzed by liquid chromatography-tandem mass spectrometry (LC-MS/MS) [40].

Microscopic observations of the flow cells were completed using a Leica TCS SP2 confocal laser scanning microscopy (CLSM) with an argon/krypton laser and detectors and filter sets for simultaneous monitoring of GFP (excitation, 488 nm; emission, 517 nm) for live cell staining and propidium iodide (excitation, 543 nm; emission, 565 nm) for dead cell staining. Images were acquired at approximately 1 *μ*m intervals in the direction down through the biofilm. Each channel of the flow chamber was randomly imaged at two separate locations per time point. CLSM images were analyzed using COMSTAT [14].

#### 3. Results

*Biofilm Growth Models: Simulation*. Figure 1(a) shows a plot of the Gompertz and Logistic models (see (7) and (5)) with parameters values , , , and time in arbitrary units. Figure 1(b) shows their corresponding time derivative (i.e., the biofilm overall rate of change), and Figure 1(c) the biofilm doubling time. Doubling time is not a constant, as for the exponential model equation (3), but is instead expressed by the following relationships:for the Gompertz, andfor the Logistic model. Note how the doubling time is only defined up to the time for which , where it tends to infinity. Figure 1(d) shows the BIC for the two models: note how the BIC is a decreasing function of , because intrinsic growth rate decreases as reaches the asymptote .