In Silico Modeling of Biofilm Formation by Nontypeable Haemophilus influenzae In Vivo

Multiple respiratory illnesses are associated with formation of biofilms within the human airway by NTHI. However, a substantial amount of our understanding of the mechanisms that underlie NTHI biofilm formation is obtained from in vitro studies. Our in silico model that describes biofilm formation by NTHI within the middle ears of Chinchilla lanigera will help isolate processes potentially responsible for the differences between the morphologies of biofilms formed in vivo versus those formed in vitro. Thus, the in silico model can be used to glean mechanisms that underlie biofilm formation in vivo and connect those mechanisms to those obtained from in vitro experiments. The in silico model developed here can be extended to investigate potential roles of specific host responses (e.g., mucociliary clearance) on NTHI biofilm formation in vivo. The developed computational tools can also be used to analyze and describe biofilm formation by other bacterial species in vivo.

However, spatial morphologies of the biofilms formed in vivo and in vitro by the same bacterial pathogen can show distinct differences. Bjarnsholt et al. (14) analyzed the sizes of bacterial aggregates in biofilms formed in parts of the human body in several diseases such as cystic fibrosis, chronic wounds, otitis media, and rhinosinusitis and concluded that the biofilm sizes in chronic infections span a much smaller range (ϳ50 to 200 m) than the biofilms formed by these bacterial pathogens in vitro where the biofilms can extend over several centimeters. Furthermore, specific morphological features such as mushroom-like structures in biofilms formed by Pseudomonas aeruginosa in vitro were absent in biofilms formed by the bacterial pathogen within the host (14,15). It has been speculated that these differences arise due to the differences between the local environment in the in vitro culture and the host, e.g., nutrient-rich culture medium versus nutrient-restricted local host environment (14). Furthermore, the host responses such as mucociliary clearance (16), production of antimicrobial proteins as effectors of innate immunity (17), and nutritional immunity (18), among others partially eliminate bacterial cells in the biofilm. These host responses are absent in vitro and can contribute to specific structural differences in biofilms formed in vivo.
In order to characterize the spatial structures of biofilms formed by NTHI in vivo and determine potential mechanisms that underlie biofilm formation that occurs within a mammalian host, we first quantitatively characterized the spatial patterns of NTHI in confocal laser scanning microscopy (CLSM) images of biofilms formed by NTHI in the middle ears of chinchillas during experimental NTHI-induced OM and then developed a three-dimensional in silico agent-based model to describe formation of NTHI biofilms in vivo. The spatial patterns were analyzed by two-point pair correlation function, which has been widely used in statistical physics and material science to characterize spatial structures. The correlation length calculated from the correlation function provides an estimation of the characteristic size of the bacterial clusters in the biofilm. The agentbased model is a modification of our in silico model developed previously to describe formation of NTHI biofilms in vitro (7). The parameters in the agent-based model were estimated by comparing architectural features (as calculated from the pair correlation function) of the in silico biofilms with those of CLSM images obtained from NTHI biofilms formed in the chinchilla middle ears (4). A variety of in silico models based on agent-based models (19), partial differential equations (20)(21)(22), and particle-based simulations (23)(24)(25) have been developed previously to describe biofilm formation of different bacterial species such as P. aeruginosa and Vibrio cholerae. However, the majority of these models dealt with biofilm formation in vitro, and thus cannot be applied directly to describe biofilm formation in vivo.
Our analysis of the CLSM images for the NTHI biofilms formed in vivo found that NTHI bacterial cells are organized in fractal structures within the biofilms. This pattern is similar to that found in NTHI biofilms formed in vitro (7); however, the size of the NTHI clusters, as defined by the correlation length, in vivo are much smaller than their in vitro counterparts. The agent-based model statistically reproduced the spatial organization of NTHI bacterial cells in the biofilm formed in vivo. Our modeling results suggest that the dispersion of biofilm-associated NTHI bacterial cells induced by quorum sensing and the elimination of the planktonic or nonadherent NTHI bacterial cells by mechanisms not present in vitro such as the host response play key roles in giving rise to the smaller NTHI spatial clusters in vivo. The model parameter values estimated for NTHI biofilms 4 days and 11 days after challenge indicate changes in the rate at which planktonic bacteria are eliminated and NTHI quorum sensing during the course of the infection.

RESULTS
Quantitative characterization of NTHI biofilms formed in vivo. NTHI formed biofilms in the middle ears of Chinchilla lanigera. Sections of the middle ear epithelium and associated biofilms 4 and 11 days after challenge were cut, immunolabeled for NTHI and extracellular DNA (eDNA), and then imaged using confocal laser scanning microscopy (CLSM) (Fig. 1A and B). Further details regarding the experiments can be found in reference 4. The green and blue colors show the NTHI bacterial cells and the DNA strands in the images, respectively. On visual inspection, the NTHI aggregates appeared to be 1 to 3 m in diameter. The size of the NTHI aggregates is characterized quantitatively later in this section. The spatial organization of the NTHI bacterial cells within the biofilm was characterized by analyzing the pair correlation function C z (r). C z (r) was calculated using the densities (or {(x,y,z)}) of the NTHI bacterial cells at locations {(x,y,z)} within the biofilm. We assumed (x,y,z) to be proportional to the image intensity at the location (x,y, z) in the CLSM image. C z (r) is defined as where L is the length of the sample along the x or y direction. The above definition of C z (r) assumes translation symmetry and rotational invariance of the bacterial density, and 〈. . .〉 indicates an average over an ensemble of density configurations. The behavior of C z (r) at length scales larger than the scale of image resolution (a 0 Ϸ 0.29 m) and smaller than an intermediate scale (w ϳ several microns) can indicate the presence of specific spatial patterns within the range as we briefly mention below. Further details regarding the analysis can be found in reference 7. When the length scale r is within the range, a 0 Յ r Յ w, C z (r) can be approximated as C z (r)/C z (r ϭ 0) Ϸ 1 Ϫ ar (26,27). The less than unity values of exponent indicate the presence of clusters with fractal surfaces (27)(28)(29). The variation of C z (r)/C z (0) at distances r Յ 1 m for the NTHI biofilms at day 4 and day 11 postchallenge shows a behavior consistent with 1 Ϫ ar , where Ͻ 1 indicates the presence of fractal interfaces in the in vivo biofilm ( Fig. 1C and D). The fit to the pair correlation data [C z (r)/C z (r ϭ 0) versus r] at small values of r (r ϭ 0 to r ϭ 1.5 m) with the function 1 Ϫ ar ( Ͻ 1) was compared with a fit with a function 1 Ϫ ar using corrected Akaike information criterion (AIC c ) (see Text S1, Table S1A, and Fig. S1 in the supplemental material). In addition, the ability of an exponential decay function [exp(Ϫar)] to describe the data was compared with that of a nonexponential decay function exp(Ϫar ) using AIC c (Text S1, Table S1A, and Fig. S2). In both the above comparisons, the data convincingly supported the nonlinear and the nonexponential fits. A piece-wise fitting function with two different exponents ( ϭ Ͻ for r Ͻ r c and ϭ Ͼ for r Ͼ r c ) fitted the pair correlation over a longer range of r from r ϭ 0 to r ϭ 5 m (Fig. S2). This demonstrates that the behavior of C z (r) indicates the presence of fractal interfaces in the biofilm. C z (r) also provides estimates of characteristic length scales (e.g., size of a bacterial cluster or a void region) present in the spatial organization. One such characteristic length scale z (Fig. 1E) can be defined using C z (r ϭ z )/C z (0) ϭ 1 ⁄2. The z values for the images in Fig. 1A estimate the sizes of the NTHI clusters. The z values for the 4-day-old and 11-day-old biofilms ( Fig. 1C and D, insets) are less than 3 m, which are substantially smaller than the large z values (ϳ20 m) obtained in NTHI biofilms formed in vitro (7).
Development of an agent-based model. We developed a spatially resolved agentbased model to describe the formation of NTHI biofilms in the chinchilla middle ear. The model was based on the agent-based model that described the formation of NTHI biofilms in two dimensions in in vitro culture wells (7). The key elements of the model are described below. The processes involved in the model are shown schematically in Fig. 2. The parameter values are shown in Table 1, and further details regarding the model simulations are detailed in Materials and Methods and in the supplemental or is occupied with a planktonic (maroon) or a biofilm-resident (green) NTHI or extracellular DNA (blue). The planktonic (purple) or biofilm-resident bacteria (cyan) can coexist in a voxel with extracellular DNA. Each voxel can be occupied by at most one bacterial cell. The physical movements of the NTHI due to diffusion or Tfp-induced motility are shown with black arrows. The prohibited movements (e.g., diffusion of planktonic NTHI to a neighboring voxel occupied by another NTHI) are marked with null symbols. The processes that arise due to NTHI replication, NTHI dispersion, attachment (or detachment) of the NTHI to (or from) the ECM, and NTHI death are depicted in the rectangular box below. material. In the model, the space is discretized into cubic "voxels" of size l 0 (ϭ0.5 m) such that each voxel can possess at most one NTHI bacterium. The NTHI bacterial cells can be either in the biofilm-resident state or in the planktonic state. The biofilmresident (or planktonic) bacterial cells are assumed to be attached (or unattached) to extracellular DNA and other (not explicitly modeled) molecules that are part of the extracellular matrix (ECM) that comprises the biofilm. Epithelial cells are also not explicitly included in the model. Bacterial growth/division occurs at a rate proportional to the density of nutrient in the simulation box; each division event decreases the amount of nutrient. Nutrition is modeled at a global level where there is a source of new nutrient into the simulation box with a specified carrying capacity. In the simulation, the amount of nutrient in each voxel increases at every time step until it reaches a maximum value n max , for each bacterial division nutrient amount in each of the voxels is decreased by 1/V, where V is the volume of the simulation box. The extracellular DNA strands are assumed to generate a mesh-like structure that was fixed at the beginning of the simulation. Based on experimental images, a rough approximation for the eDNA was used: parallel DNA strands (going along the x direction) that are 1 m thick and spaced by 5 m, there are some random breaks in the strands (5% probability every 5 m) and there are random bridges between strands in the plus or minus y and z directions (5% probability for each bridge every 5 m). This assumption was implemented to investigate the role of the eDNA mesh structure in affecting the biofilm morphology. The assumption ignores the details regarding the formation of extracellular DNA mesh, as the biofilm forms within the host and the extracellular DNA strands are produced by NTHI. Since precise molecular details regarding the contribution of the extracellular DNA generated by the NTHI and the dead cells are not well-known, we avoided modeling these processes explicitly to protect the model from overparameterization. We modeled the following processes in our model. (i) NTHI movements. Three different processes produce movement of NTHI bacterial cells in our model. (i) When the number of NTHI in a voxel exceeds a maximum packing number (n thres ϭ 1 bacterial cell), the excess bacteria are moved to an adjacent compartment that can accommodate the excess bacteria without going over the threshold n thres . This rule represents passive movement of NTHI due to the mechanical forces exerted by the neighboring bacteria in a tightly packed region. (ii) In our simulation, the biofilm-resident bacteria on extracellular DNA can "twitch" to move along the extracellular DNA network. This rule was constructed motivated by observation of twitching motility of P. aeruginosa on extracellular DNA tracks in phase-contrast imaging experiments (30). However, such movements have not yet been directly observed for NTHI. Previous experiments showed that NTHI exhibits twitching motility on solid or biotic surfaces which is not included in the current model. Therefore, the movement of NTHI modeled here (twitching along eDNA strands) is a hypothesis used in the in silico model. This movement is velocity driven and nondiffusive (ballistic) on short time scales and always moves in a direction with extracellular DNA strands (30). (iii) Planktonic NTHI bacterial cells can move diffusively between the voxels.
(ii) NTHI dispersion. Quorum signaling brings the bacteria together and induces dispersion of the bacterial cells. NTHI dispersion changes biofilm-resident NTHI bacterial cells to planktonic bacterial cells. The dispersion is activated in a voxel when the number of NTHI bacterial cells in neighboring voxels (at a range of d quorum ) reaches a threshold (n quorum ). This step approximates and effectively models dispersion (at a rate k quorum ) arising due to attaining a threshold number of quorum-sensing-mediating molecules (modeled implicitly) that are secreted by NTHI in the local neighborhood.
(iii) Attachment/detachment of NTHI bacterial cells. Biofilm-resident (or planktonic) NTHI bacterial cells can detach (or attach) from (or to) the ECM to become planktonic (or biofilm-resident) bacterial cells.
(iv) Elimination of NTHI cells by the host response. The multifactored host response, including mucociliary clearance, antimicrobial proteins, neutrophil extracellular trap (NET) formation, nutritional immunity, and the adaptive immunity that eliminates the NTHI is not explicitly modeled. We represented the effect of the host response by removal of the planktonic NTHI bacterial cells at a higher rate compared to their replication rate. The elimination of the biofilm-resident NTHI cells by the host response was assumed to occur at a much lower rate and was ignored in the model.
(v) Simulation of the model with time. The model is progressed in time by executing the above processes with specific probabilities following a standard rejection kinetic Monte Carlo (KMC) algorithm (7). The parameter values used in the model are shown in Table 1. Further details are provided in Materials and Methods.
(vi) Estimation of model parameters. Parameters marked with asterisks in Table 1 are estimated using the C z (r) calculated from the CLSM images. We use a particle swarm optimization to find the best-fit parameter values. Details are provided in Materials and Methods.

Agent-based model describes the organization of NTHI bacterial cells in vivo statistically.
The agent-based model was used to generate spatial organization of NTHI bacteria in three dimensions at different days postchallenge. The simulation was initiated by the introduction of a few (ϳ2) NTHI in two voxels within the simulation box of size 50 ϫ 50 ϫ 50 m 3 . These bacteria replicate and migrate following the processes described above and generate clusters of biofilm-resident NTHI bacterial cells. The size of the cluster increases initially. However, as the size of the cluster increases, the effect of the quorum-sensing mechanism is activated in the model, increasing the dispersion of biofilm-associated bacteria to planktonic bacteria. Most of the planktonic bacteria are removed by the simulated effective host response. However, a small fraction of these bacteria become reattached to the ECM and generate new clusters of biofilm-resident NTHI. The hypothesized twitching movements of the NTHI on the strands of extracellular DNA mesh also spread out the biofilm-resident bacterial cells. Therefore, the size of the NTHI clusters start decreasing after initial growth, as many smaller clusters form and the system is no longer dominated by a single large cluster ( Fig. 3A to C). The biofilm-resident NTHI bacterial cells settle into such configurations for a very long time, and the system appears to have reached a more spatially inhomogeneous steady state ( Fig. S4 and Movie S1). The above kinetics can be characterized by the initial increase followed by a decrease and then saturation to a fixed value in the correlation length, z (t) (Fig. 3D). The values of the cluster sizes [or z (t)], the timing of its peak value, and the cluster size at long times depend on the values of the parameters in Table 1, though the qualitative nature of above-described changes remain the same for wide ranges of these parameter values.
We conducted a simulation of the model without any DNA strands. Since in the model we stipulated that NTHI twitch only along the DNA strands, this twitching motility will be absent in the model when DNA strands are not present. In the simulation that did not contain any extracellular DNA, the bacterial cells were less spread out compared to the simulation that contained extracellular DNA strands (Fig. S3). There are two important points to note here. (i) We reiterate that previous work has shown that NTHI exhibits twitching motility on both biotic and abiotic  Table 1, and Table 2 shows the best-fit parameters to the CLSM image obtained at 11 days postchallenge. (D) Change in the correlation length z calculated from a slice though the center of the box Ϯ1 m (where the bacteria were initialized), with time for the model parameter values associated with the configurations in panels A to C. surfaces, and there is no biological evidence that NTHI twitching motility requires the presence of extracellular DNA. The twitching motility along the extracellular DNA in the model is based on data from studies with P. aeruginosa, but such movements have not yet been directly observed for NTHI. (ii) In the simulations without any DNA strands, we found that the bacterial cells did not spread out and form smaller clusters in large numbers as can be seen in the simulation with DNA strands. Note that both the DNA strands are exported and that the Tfp are expressed through the ComE-comprised pore in the outer membrane of NTHI, and thereby, a comE mutant of NTHI will lack both the Tfp-induced twitching activity and the ability to export DNA strands (8). Thus, the modification of the model where DNA is absent from the biofilm and the bacteria are unable to execute any twitching movements (which occurs as a consequence of the model assumption stated above) might appear to describe biofilm formation by the comE mutant as a consequence of these linked biological activities. However, we have not yet tested the model predictions against any in vivo experiments.
Estimation of parameters not available from experiment. We estimated parameter values for eight different parameters (indicated with an asterisk in Table 1) whose values are not known from laboratory measurements by fitting the pair correlation function C z (r) for a specific z value at a particular day postchallenge (e.g., 11 days) that is calculated from the NTHI bacterial cell configurations in the agent-based model with the C z (r) calculated from the corresponding CLSM image. Since there was no reference to match the z value between the in silico and in vivo configurations, we arbitrarily chose to fit the NTHI bacterial cell configurations in the center z-stack (e.g., z ϭ L/2) of the CLSM image. The arbitrariness in choosing the z value for comparison does not affect the general conclusions drawn from the parameter estimations, as z (t) values were quite similar between different z-stacks in the C z (r) calculated from the CLSM images (Fig. 1C, inset). The particle swarm optimization method obtained a low value for the cost function (E) in equation 2, indicating an excellent fit for the NTHI bacterial cell organization at 11 days postchallenge. Since we fit the pair correlation function with the in vivo images, the agreement refers to statistically similar NTHI bacterial cell configurations between the model and the NTHI bacterial cell configuration in the chinchilla middle ears. The representative in silico and in vivo configurations are shown in Fig. 4A and B. The fits to the C z (r) are shown in Fig. 4C.
The estimated parameter values and the standard deviations of those values are shown in Table 2.
The rates in Table 2 show changes in the simulation parameters as the infection progresses from day 4 to day 11. The rate of elimination of the planktonic NTHI (or k death ) by the host response as well as the rate of dispersion (or k quorum ) induced by quorum sensing increase from day 4 to day 11, increasing the resistance of the host response on NTHI growth on the later days. This mechanism is consistent with the slight decrease (from ϳ1 m to ϳ0.5 m) of the NTHI spatial cluster sizes as the NTHI biofilm age increases from day 4 to day 11 (Fig. 1). We note the following caveat in the above parameter estimation. The precise values of the parameters appear to depend on the choice of the initial bacterial organization. For example, when a biofilm configuration obtained at day 4 for a model that fit the 4-day-old biofilm data was further simulated using the parameters obtained for the fit to the 11-day-old biofilm (Table 2), the biofilm configurations in the simulation at day 11 (Fig. 3B) agreed qualitatively with the data at day 11 (Fig. S5) in that the clusters of bacterial clusters become more diffuse, and the difference in the amount of bacteria between Fig. 3B and Fig. S5 is due to the different initial conditions. However, the behavior of the pair correlation function did not match with that obtained from the data well (Fig. S6A). The calculation of two-time autocorrelation function (31), A(t 1 ,t 2 ) ϭ〈(r,t 1 )(r,t 2 )〉 , where t 2 Ͼ t 1 , for these simulations seemed to indicate a power law decay (A(t 1 ,t 2 )ϰ 1/(t 2 ) f(t 1 ,t 2 )) with time (Fig. S6B). The power law decay of the autocorrelation function points to a strong dependence of the system kinetics on the initial condition which has been observed in many models of ordering kinetics in statistical mechanics (31). However, the mechanisms that underlie the initial condition dependence of the spatial kinetics in our model are currently unclear and need further investigation.

DISCUSSION
The analysis of spatial organization of NTHI in biofilms formed in chinchilla middle ears using pair correlation functions revealed the presence of surface fractal patterns at short length scales (ϳ1 m). Similar fractal patterns were also found in the biofilms formed by NTHI in in vitro cultures (7). The presence of such surface fractal patterns indicates scale invariant organization of NTHI within the biofilm, which can potentially facilitate absorption of nutrients by increasing the surface area of the interfaces between the bacterial cells and the local environment in the nutrient-restricted middle ear. The increased surface area could also further expose the bacterial cells to antimicrobial proteins or other immune effectors. However, the ECM in the biofilm has been found to provide protection to biofilm-resident NTHI against the host immune response (9), and therefore, the survival of the NTHI in vivo could depend on finding an optimal pattern of spatial organization of NTHI within the biofilm that balances the above trade-offs. Our analysis also showed that the sizes of the NTHI spatial clusters are much smaller (ϳ1 m) compared to their in vitro counterparts which can be more than 20 m in diameter (7). As we discuss below, our modeling showed that these smaller size clusters could arise due to the elimination of planktonic NTHI bacterial cells due to the host immune response (or other means that eliminate planktonic bacteria). A similar decrease in the size of aggregates of bacterial cells in biofilm patches formed within the host by bacterial pathogens in chronic infections ranging from cystic fibrosis to rhinosinusitis compared to biofilms formed in vitro has also been observed (14).
We developed an agent-based model including basic processes such as NTHI replication, quorum sensing, dispersion of biofilm-resident NTHI to planktonic NTHI, and elimination of planktonic NTHI by the host immune response to describe formation of NTHI biofilms in three dimensions in the chinchilla middle ears. The model was a modification of our previously developed model in two dimensions (7) that described formation of NTHI in culture. The agent-based model incorporated most of the processes that were present in our previous agent-based model except the ones that were associated with the feeding protocol. The key modifications made in the agent-based model to describe biofilm formation in vivo were inclusion of the elimination of planktonic NTHI by the host immune response, no dependence of the nutrient production on physical orientations of biofilm, and creation of a nutrient-restricted environment. The model was able to reproduce the surface fractal patterns and the small sizes of NTHI clusters in the biofilms formed in vivo (Fig. 4). The perturbations of the model parameters showed that in the absence of the elimination of the planktonic NTHI, the NTHI bacterial cells formed large clusters and filled up the simulation box (see Fig. S7 in the supplemental material), demonstrating the key role elimination of planktonic cells (such as by the host response) plays in limiting the NTHI cluster sizes in the model. This behavior suggests the presence of the host response, which can be comprised of mucociliary clearance, antimicrobial proteins, neutrophil extracellular traps (NETs), and nutritional immunity, among others, each of which is likely important in limiting the size of the bacterial clusters in biofilms formed within the host. The agent-based model does not explicitly include the layer of epithelial cells in the middle ear, which produce antimicrobial proteins or induce additional protective host responses (e.g., mucociliary movements) that impact relative survivability of NTHI. Analyzing the roles of these processes in shaping the organization of NTHI within the biofilm can provide key mechanistic insights regarding the transport of antimicrobial proteins and antibody molecules induced by vaccine candidates. These processes will be included in future agent-based models.
We estimated several parameters in the agent-based model using the pair correlation functions calculated from the CLSM images of the NTHI biofilms formed in the chinchilla middle ears. The experimental measurements for some of these parameters were unavailable, and the rest of these parameters were associated with processes that approximated detailed molecular interactions and thus are inaccessible to direct measurements. Similar parameter estimations are common in well-mixed or spatially homogeneous models of biological processes (32,33); however, such estimations are rare for spatial models of biofilm formation because these calculations can be computationally intensive. A recent work estimated parameters describing pairwise interactions between bacterial cells residing in biofilms formed by V. cholerae in vitro (25). We took advantage of parallel computation and used an optimization method (particle swarm) that can be easily parallelized to deal with this problem. Our estimations showed most of these parameters can be well estimated from the spatial information contained in C(r). The comparison of the parameters for the day 4 and day 11 CLSM images for the NTHI biofilms showed that the rate of elimination increases from day 4 to day 11; thus, the smaller size NTHI clusters at day 11 compared to those at day 4 could arise due to this effect.
In conclusion, we developed an in silico model for describing formation of NTHI biofilms in the chinchilla middle ears. The in silico model can be used to probe changes in the biofilm morphologies formed by NTHI over time in vivo and relate those changes to biologically relevant parameters (e.g., immune response, quorum sensing). In addition, the model can be used to test therapeutic strategies for disrupting biofilms and/or preventing their formation, e.g., increase the rate of NTHI elimination by increasing the NTHI dispersion rate from the biofilm or preventing them from entering into biofilm residence.

MATERIALS AND METHODS
Simulation of the agent-based model. The model was progressed in time using a threedimensional generalization of the standard rejection kinetic Monte Carlo (KMC) scheme described in reference 7. The time constant for the rejection KMC algorithm is based on the maximum rate in the system, k max , which in this case is k plank . Each KMC trial time is advanced by ⌬t ϭ Ϫ log͑u͒⁄4N bact k max where u is a uniform random number in the interval (0,1] that is generated each MC trial. The factor 4 in the denominator of the above expression results from the number of possible processes (movement, growth, death/turning, attachment). For the chosen parameters above, on average, 0.0521 s passes every time all bacteria in the system attempt one (randomly chosen) process one time. At each KMC trial, the nutrient density in a voxel goes up by Δn, given by ⌬n ϭ k nutrientͩ 1 Ϫ n n max ͪ ⌬t and each time a bacterium divides (at a rate of nk div ), the nutrient in all the voxels decrease by 1/V, where V is the volume of the box.
(i) Initial condition. All the simulations are started by having two biofilm-attached NTHI bacterial cells in the center of the domain. The extracellular DNA network is created as described above and is held fixed throughout the simulation.
(ii) Boundary conditions. Periodic boundary conditions along all the three spatial directions are used.
(iii) Particle swarm optimization. Parameters listed in Table 1 that are marked with an asterisk were optimized (on a logarithmic scale when range given as powers of 10) using that standard particle swarm optimization algorithm (34,35) with constrains as implemented in pyswarm which is available at the link https://github.com/tisimst/pyswarm. The parameters used were: particle inertia coefficient of ϭ 0.8, coefficient weighting a particle's best-known position of p ϭ 2.5, and coefficient weighting the swarm's best-known position of p ϭ 1.5. The swarm size was 200 particles for the 11-day fit and 400 particles for the 4-day fit (due to the shorter simulation time), and the algorithm was run for 40 iterations.
The optimization was done to match the parameters of a fit of the pair correlation function to the following equation: C z͑ r ͒ ⁄ C z͑ r ϭ 0 ͒ ϭ 1 Ϫ ar (1) and z is the correlation length defined as C z ͑r ϭ z ͒ ϭ 1⁄2.
The cost function used was E ϭ 1 2 iϭx,y,z where the "exp" subscript indicates an experimental value fit from a single z-stack of an 11-day-old biofilm, and the in silico fits were from a 5-voxel-wide slice down the middle of the simulation box in the x, y, or z direction, as indicated. All three directions in the model were used for the comparison, as the orientation of the CLSM image was unknown. The sizes of the bacterial clusters in our simulations are smaller or of the same order of the smallest spacing in the eDNA mesh. The pair correlation function decays substantially across that length scale and thus does not depend on the direction. The fit parameters are shown in Table 2. The errors in the parameter values were calculated using standard deviations of points in solution space sampled with E Ͻ 0.01. Using the clustering algorithm in reference 36 (normalizing all parameters to a unit hypercube, and using a cutoff of 0.4), only a single cluster was found, therefore, we concluded that all points in that cluster were in the same basin. The ranges in Table 2 are the standard deviations of the locations of these points in solution space since these points give a rough approximation of the dimensions of this basin. This allowed us to distinguish between parameters that are similar and dissimilar. Experiments. Chinchilla lanigera were challenged transbullarly with 2,500 CFU NTHI strain 86-028NP as described previously (4). Four or 11 days later, the bullae, which contained biofilms formed by NTHI, were removed, embedded in OCT, and snap-frozen over liquid nitrogen. After removal of the exterior bone, we cut serial 10-m sections of the middle ear epithelium and associated biofilms, immunolabeled them for NTHI with anti-NTHI outer membrane protein (OMP), and counterstained with 4=,6=-diamidino-2-phenylindole (DAPI) to label eDNA. Samples were imaged by CLSM. Further details are provided in reference 4.