Extracellular DNA and Type IV Pilus Expression Regulate the Structure and Kinetics of Biofilm Formation by Nontypeable Haemophilus influenzae

ABSTRACT Biofilms formed in the middle ear by nontypeable Haemophilus influenzae (NTHI) are central to the chronicity, recurrence, and refractive nature of otitis media (OM). However, mechanisms that underlie the emergence of specific NTHI biofilm structures are unclear. We combined computational analysis tools and in silico modeling rooted in statistical physics with confocal imaging of NTHI biofilms formed in vitro during static culture in order to identify mechanisms that give rise to distinguishing morphological features. Our analysis of confocal images of biofilms formed by NTHI strain 86-028NP using pair correlations of local bacterial densities within sequential planes parallel to the substrate showed the presence of fractal structures of short length scales (≤10 μm). The in silico modeling revealed that extracellular DNA (eDNA) and type IV pilus (Tfp) expression played important roles in giving rise to the fractal structures and allowed us to predict a substantial reduction of these structures for an isogenic mutant (ΔcomE) that was significantly compromised in its ability to release eDNA into the biofilm matrix and had impaired Tfp function. This prediction was confirmed by analysis of confocal images of in vitro ΔcomE strain biofilms. The fractal structures potentially generate niches for NTHI survival in the hostile middle ear microenvironment by dramatically increasing the contact area of the biofilm with the surrounding environment, facilitating nutrient exchange, and by generating spatial positive feedback to quorum signaling.

N ontypeable Haemophilus influenzae (NTHI) is a predominant bacterial pathogen of otitis media (OM), a common middle ear infection for children worldwide under the age of 15 years (1). Chronic and recurrent OM are both characterized by the formation of biofilms in the middle ear, often by NTHI (1,2). Biofilms are threedimensional complex structures formed by bacteria and are usually associated with either abiotic or biotic substrates, including mucosal surfaces (3). Biofilm formation contributes to bacterial resistance to host immune responses or antibiotics. Distinct spatial patterns, such as tower or mushroom structures, arise in biofilms formed by different bacterial species under specific experimental conditions. Such spatial patterns can be associated with survival of the bacteria within the biofilm. For example, mushroom-like structures separated by water-filled regions occurred in biofilms formed by Pseudomonas aeruginosa in flow chambers when glucose was used as a carbon source in the medium (4). The mushroom structures were proposed to create a circulatory system for efficient nutrient supply and waste removal (5,6). Similar morphological features may promote NTHI survival in the hostile host environment. Therefore, characterization of specific spatial patterns in NTHI biofilms and uncovering the mechanisms that give rise to those morphological features are important for understanding the establishment and pathogenesis of OM and for developing therapeutic strategies for treating OM.
Recent studies have identified several factors that play important roles in regulating the stability and morphology of NTHI biofilms (7)(8)(9)(10)(11)(12)(13)(14)(15)(16)(17). One such key factor is extracellular DNA (eDNA) (8,11), a major structural component of the extracellular polymeric substance (EPS) that comprises the biofilm matrix. The eDNA within an NTHI biofilm forms a tight meshwork of criss-crossing eDNA strands (9,11), which is stabilized by DNABII proteins bound to the vertices of crossed eDNA strands. The origin of eDNA in bacterial biofilms remains unclear, but it is often attributed to bacterial autolysis (18) or release of outer membrane vesicles. However, a recent study by Jurcisek et al. (19) described a novel mechanism: DNA and DNABII protein release occurs via a newly identified inner membrane pore complex (TraC and TraG) with homology to type IV secretion-like systems and through the ComE pore, which is part of the machinery for expression of type IV twitching pili and is located in the bacterial outer membrane. The role of the eDNA network in bacterial cell organization during the process of biofilm formation is not well understood. High-resolution time-lapse microscopy studies of Pseudomonas aeruginosa biofilms showed that the eDNA network facilitates mass movement of bacterial cells along "furrows" within the biofilm structure (20). This movement is mediated by type IV pili (Tfp), which execute a repeating pattern of extension, attachment, and retraction that results in the ratchet-like motion termed twitching motility. Mokrzan et al. (13) demonstrated that NTHI Tfp are required for the formation of "tower" structures characteristic of NTHI biofilms. Taken together, these data suggest that both the eDNA network and Tfp expression play important roles in NTHI biofilm structural development.
However, a quantitative characterization of spatial patterns found within biofilms formed by NTHI and its mutant strains, and how these patterns are regulated by the eDNA network and Tfp expression, are not well understood. We addressed the above issues by combining data analysis approaches rooted in statistical physics (21)(22)(23), confocal microscopy imaging of NTHI biofilms formed in vitro, and in silico mathematical modeling (24,25). Our analysis of the confocal images of biofilms formed in vitro by wild-type (WT) NTHI revealed formation of fractal structures in the organization of wild-type NTHI cells within the biofilms at short distances (Յ10 m). Emergence of fractal structures could potentially help bacteria survive in a hostile environment (e.g., the middle ear) by facilitating nutrient absorption or by locally trapping molecules that mediate quorum sensing. The agent-based in silico model showed that the emergence of such structures in wild-type NTHI biofilms is facilitated by the eDNA network and Tfp function. This was further confirmed by the agreement of model predictions regarding the large reduction in these fractal structures in biofilms formed in vitro by a ΔcomE mutant strain of NTHI. The ΔcomE mutant has impaired Tfp function, and the biofilms formed by the ΔcomE strain in vitro contain substantially less eDNA than those formed by the wild type. The developed in silico model can be used to test various hypotheses regarding biofilm formation by mutant NTHI strains and perturbation of biofilm structures by different drugs.

Analysis of the pair correlation function in bacterial density reveals the presence of fractal structures of smaller length scales (<10 m) in older (>40 h) NTHI biofilms formed in vitro.
We analyzed spatial pair correlations between bacterial densities in confocal microscopy images for biofilms formed in static culture by NTHI strain 86-028NP in order to determine the presence of specific spatial patterns in the morphology. The local intensity in the images was assumed to be proportional to the local bacterial density (x,y,z,t), where {x,y,z} denotes the spatial location within the biofilm and t indicates the age of the biofilm (Fig. 1A). We calculated the pair correlation function C z (r,t) between the bacterial densities at two locations separated by distance r in a plane x-y parallel to and residing at a fixed distance (z) from the substrate (Fig. 1A) (see Materials and Methods for further details on our procedure). The behavior of C z (r,t) at short distances (r), characterized by exponent , was used to determine the nature of the spatial patterns for that length scale (see Materials and Methods). The value of exponent contains information regarding the nature of the spatial pattern for that length scale (21)(22)(23)26). For example, if the bacterial cluster in the two-dimensional x-y plane forms a scale-invariant structure, a mass fractal (Fig. 1B) or a surface fractal (a solid The three-dimensional density profile [(x,y,z,t*)] for NTHI from the confocal microscopy imaging data, which were obtained from a sample of an in vitro-grown NTHI biofilm at t* ϭ 88 h. The NTHI density at each z plane is shown in grayscale, with darker intensities indicating higher densities. The biofilm surface is shown in color. (Right) NTHI density (x,y,z*,t*) in a z plane (e.g., z* ϭ 14 m) was analyzed by computing the pair correlation between a pair of points ( r structure with a fractal interface) (Fig. 1C), or a solid structure with a sharp nonfractal interface (Fig. 1D), the values of will be Ͻ0, between 0 and 1, and equal to 1, respectively. Visual inspection of a two-dimensional density profile of the confocal images at a z plane for an NTHI biofilm grown for 88 h suggested the presence of fractal structures at the interfaces of the bacterial clusters with the surrounding liquid environment (Fig. 1E). To investigate the presence of such structures, we analyzed C z (r,t*) for the biofilms formed at t* ϭ 16 h, 40 h, 64 h, and 88 h after biofilm growth was initiated. C z (r,t*), calculated for the biofilm at many planes parallel to the substrate at t* ϭ 88 h, could be fitted with values smaller than unity for distances ( Fig. 2A and B; see also Fig. S1 in the supplemental material), with r between 0.27 and 10 m, indicating the presence of surface fractal structures at shorter distances (r Ϸ 10 m). Similar behavior in C z (r,t*) was present at other time points (e.g., t* ϭ 40 h, 64 h) as well. At 16 h, most of the bacterial clusters showed the presence of surface fractal structures at many z planes ( Fig. 2C and D) but also displayed solid structures with sharp but smooth interfaces ( Ϸ 1) ( Fig. 2C and D) at several z planes. Thus, the analysis showed the presence of fractal interfaces within short length scales (Յ10 m) in older (Ն40 h) biofilms.
Existence of multiple time and length scales in the biofilm structures. Next, we investigated whether there was a characteristic length and time scale associated with FIG 2 Surface fractal structures in NTHI biofilms. (A) Variation of C z (r,t*)/C z (r ϭ 0,t*), with r at a fixed t* (88 h) and at three different z planes within a replicate of the in vitro NTHI biofilm. The exponents obtained by fitting the function in equation 2 for the data at z ϭ 2, 6, and 12 m were 0.25, 0.28, and 0.33, respectively. The fits are shown with solid lines. The straight line shown in black gives an example for the case with ϭ 1. (B) Values of the estimated exponents at different z stacks, calculated by fitting equation 2 to C z (r,t*) for r Յ1.5 m, calculated for 8 different replicates for the 88-h NTHI biofilm. The different symbols indicate different replicates. All the estimated values for were smaller than 1, indicating the presence of scale-invariant or fractal structures at short length scales. The symbol colors (see the color bar on the right) display the 2 values, which indicate the quality of the fit (equation 3) used to estimate the values. A lower 2 value indicates a better fit of the data with equation 3. (C) Variation of C z (r,t*)/C z (r ϭ 0,t*) with r at a fixed t* (16 h) and at three different z planes within a replicate of the in vitro NTHI biofilm. The estimated exponent values at z ϭ 2, 6, and 8 m were 0.54, 0.71, and 0.72, respectively. The solid straight line in black shows the reference for the ϭ 1 case. (D) Estimated values for at different values of z for different replicates of the 16-h biofilm. Unlike the 88-h biofilm, the estimated value was closer to 1 for several z planes, suggesting the presence of sharp but smooth interfaces in those planes within the biofilm. The visualization scheme is the same as that for panel B.

Das et al.
® the organization of the bacterial cells within the biofilm. The presence of such a behavior points to the relevance of specific processes in determining the spatial structure in the system of interest. For example, theoretical models in statistical physics (27) describing growth of ordered magnetic domains display emergence of a characteristic length and time scale associated with the size of the ordered magnetic domains, where the domain growth is driven by the surface tension of the domains. In the presence of such a characteristic length scale z (t), C z (r,t) at a fixed z (ϭz*) but different times (e.g., t 1 and t 2 ) behaves as C z* (r,t 1 )/C z* (0,t 1 ) ϭ f z [r/ z* (t 1 )] and C z* (r,t 2 )/C z* (0,t 2 ) ϭ f z* [r/ z* (t 2 )]. Thus, C z* (r,t) calculated for density profiles collected at two different times overlaps on the same scaling function f z* (x), where r is scaled as r ¡ r/ z* (t). The correlation length z* (t) is calculated using the equation, C z* [r ϭ z* (t),t]/C z* (0,t) ϭ 1/2, i.e., z* (t) describes a length scale at which the normalized density correlation decreases to 1/2 and can be related to sizes of the NTHI-rich regions in the z ϭ z* plane. The analysis of the biofilm structures showed an absence of such a scaling form (Fig. 3A), revealing the existence of multiple length scales associated with different time scales. Next, we checked if the structures at different z planes at a fixed time (t ϭ t*) contained a characteristic length scale, which would be reflected in a scaling behavior as C z (r,t*)/ C z (0,t*) ϭ g t* [r/ z (t*)]. In the presence of such scaling, the graphs showing the variation of C z (r,t*) with r at multiple z planes at a fixed time t (ϭt*) will overlap on the same scaling function g t* (x) when r is scaled as r ¡ r/ z (t*). We found an absence of any scaling behavior (Fig. 3B), indicating the presence of multiple length scales at different z planes in the same biofilm. The variation of z (t*) with z showed a slow increase as the distance from the substrate increased and then a slight decrease near the top of the biofilms (Fig. 3B, inset). This behavior revealed the presence of a denser organization (bacterial cluster size of Յ1 m, void size of Յ1 m) of the bacteria near the substrate that segued into structures with larger clusters (Ն5 m) of bacteria separated by larger voids (Ն20 m) at intermediate distances (~30 m) from the substrate.
A predictive in silico model describes biofilm formation in vitro by wild-type NTHI. The above analysis of the in vitro NTHI biofilm morphologies showed the presence of specific spatial patterns (e.g., surface fractal structures) and the absence of any characteristic length scale, implying that the multiple processes spanning a range of time and length scales need to be considered when developing an in silico model for describing biofilm formation. We detail below a coarse-grained agent-based model that we were able to use to describe qualitatively the emergence of the above morphological patterns in the biofilms formed by wild-type NTHI and by a mutant strain of NTHI.
The agent-based model described the biofilm growth on a quasi-two-dimensional simulation box representing biofilm growth in the x-z plane on a substrate that lies in the z ϭ 0 line (Fig. 4A). The simulation box of size L x (ϭ128 m) by L z (ϭ40 m) by l 0 (ϭ1 m) is discretized on a regular lattice composed of small cubic compartments of size l 0 by l 0 by l 0 (l 0 ϭ 1 m) (Fig. 4A). The model describes the biofilm growth in terms of a set of agents, namely, NTHI bacterial cells associated with the biofilm, planktonic NTHI bacterial cells, eDNA strands, and nutrients. These agents change in time following a set of rules that represent key processes involved in NTHI biofilm growth in vitro. The time evolution in the model is implemented using a kinetic Monte Carlo scheme (25). The details of the simulation method are provided in Materials and Methods and Text S1, and details regarding the parameter values, model assumptions, and their connection with known experimental results are shown in Table 1 and also in the supplemental material. Briefly, processes that were modeled in the simulation were the following. (i) For modeling NTHI replication and movement in the biofilm, individual NTHI bacterial cells residing in a compartment consume nutrients within the compartment and replicate at a given rate. The replication of NTHI ceases if the nutrient concentration in a compartment falls below a threshold, c thres (ϭ1 molecule/m 3 ). When the number of biofilm-associated NTHI cells in a compartment exceeds a maximum packing number (n thres ), the excess bacteria are moved to an adjacent compartment that can accommodate the excess bacteria without exceeding the threshold n thres (Fig. 4A). This rule represents passive movement of NTHI due to the mechanical forces exerted by the neighboring bacteria in a tightly packed region. (ii) For modeling development of the eDNA network, eDNA strands can be present in the compartments. The eDNA strands are produced by a small fraction of NTHI cells at a fixed rate. The simulation is started with a small amount of eDNA present in the system. The eDNA strands that are in contact with the substrate become immediately connected to the substrate and initiate the formation of an eDNA network. The substratebound eDNA strands can irreversibly bind to eDNA strands in neighboring compartments and spread the network spatially. The eDNA strands in the eDNA network are considered immobile, and eDNA strands that are not connected to the eDNA network can diffuse into the neighboring compartments (Fig. 4A). Exact mechanisms behind the formation of the eDNA network are not known. DNABII proteins play a role in stabilizing the eDNA network (8); however, precise mechanisms for this stabilization are not clear. This rule is used to minimally develop an eDNA network in the model. (iii) To model NTHI dispersion, we consider that NTHI bacterial cells in the biofilm can disperse into the surrounding liquid medium. We hypothesized that the rate of dispersion is lower in a compartment that is a part of the eDNA network (Table 1). (iv) For modeling movement of NTHI within the eDNA network, we hypothesized that the NTHI bacterial cells move between adjacent compartments when the compartments belong to the eDNA network (Fig. 4A). This process represents a potential Tfp-driven movement of NTHI on eDNA strands. (v) To model nutrient injection and diffusion, a nutrient is introduced into the system twice in 24 h following the 16 h plus 8 h schedule used in the experiments. The nutrient diffusion is much faster than the other rates considered in the model. Therefore, we homogenized the nutrient concentration spatially at every time step in the simulation (see Text S1 for more details). (vi) In modeling diffusion of planktonic NTHI cells, the NTHI bacterial cells in the liquid medium are considered to be in the planktonic state and to move between the compartments via Brownian motion (Fig. 4A). (vii) For modeling removal of planktonic bacteria, we removed 80% of the planktonic bacteria each time the nutrient was added in the simulation. This was done to mimic the step in the experimental protocol where the older medium was removed by aspiration and replaced by fresh medium. Removal of the medium also removed the bacteria in the supernatant from the culture well. Values for the rates and other parameters used in the modeling, as well as the supporting references from the literature, are provided in Table 1.
The NTHI bacterial cells were inoculated at time t ϭ 0 into the compartments along the z ϭ 0 line (substrate) with a low concentration of bacterial cells per compartment (ϭ0.2 cells/m 3 ). We simulated the biofilm formation until 88 h. The NTHI cells  inoculated into the compartments replicated and moved to neighboring compartments once the number of NTHI cells in a compartment increased beyond the threshold number. The NTHI cells also spread due to their movements along the eDNA network. This resulted in spatial growth of NTHI cells on the substrate (Fig. 4B to D). eDNA strands were produced by a small subset of NTHI cells. These eDNA strands diffused into the neighboring compartments or became connected to the eDNA network (Fig. S2). The profiles of the biofilms generated in the in silico model were qualitatively similar to their in vitro counterparts (Fig. 4E to G) in the following ways.  (Fig. 4B to G) to generate smoother biofilm surfaces. In addition, at later times (Ն40 h), uneven structures at short length scales (Յ5 m) occurred in the biofilms formed in vitro and in silico, where regions with low and high NTHI bacterial cell numbers coexisted with each other (Fig. 4B to G). These similarities are further quantified in the section below. The dispersion of the NTHI cells to the planktonic state resulted in removal of NTHI from the system during the feeding times.
The slower dispersion of the bacteria from the compartments containing the eDNA network generated a spatially uneven organization of NTHI within the biofilm (Fig. 4B to D). Each model simulation produced a microscopically different biofilm configuration due to the statistical differences in the initial inoculation configurations and the stochastic nature of the dynamic processes (e.g., NTHI replication, eDNA production) involved in the biofilm growth. Next, we quantitatively compared the biofilm structures in the simulation with those formed in the in vitro culture. We validated the model simulations against the biofilm experiments by comparing, (i) the kinetics of the average biofilm height (h avg ), (ii) the kinetics of the coefficient of variation of the biofilm profile, (iii) the variation of C z* (r,t) with r at a fixed z plane (z ϭ z*) at different times (t) of biofilm growth, and (iv) variation of C z (r,t*) with r at multiple z planes for a fixed t (ϭt*). We focused on qualitatively comparing the morphological patterns between the biofilms formed in silico and in vitro, e.g., whether surface fractal structures were present or whether the values were less than unity in both these biofilms. This was more appropriate than a quantitative comparison in our case, because a quantitative comparison of the spatial patterns would compare both biofilms in terms of their associations with a specific type of a surface fractal or with a particular value of . Since the in silico model provides a reduced or coarse-grained description of the biofilms formed in vitro, many microscopic details responsible for generating the precise nature of the surface fractals are not explicitly present in the in silico model, and thus the in silico biofilm morphologies are unsuitable for quantitative comparisons.
We simulated NTHI growth for 88 h and calculated the above quantities based on the in silico morphology of the NTHI biofilm at regular time intervals. The above quantities (e.g., h avg ) were averaged over at least 1,000 in silico biofilm configurations obtained from independent simulation runs that were initiated with the same model parameters but with different random realizations of the initial NTHI inoculation. The biofilm height in the model was determined by finding a biofilm surface (the number of NTHI bacterial cells in a compartment that crossed a specific threshold value the first time) as the biofilm was approached from the top along the z axis (see Materials and Methods). The two-dimensional biofilm surface in confocal images of the biofilms was found similarly using a threshold in the image intensity. The coefficients of variation (CVs) for the biofilm surfaces were calculated using following definition, CV ϭ /h avg , where is the standard deviation for the fluctuations in the biofilm surface height across the average value h avg (see Materials and Methods and Text S1). The coefficient of variation is a dimensionless parameter that quantifies the magnitude of fluctuations in the biofilm surface height relative to the average height. A larger coefficient of variation is indicative of larger fluctuations in the biofilm height across the average value. The average biofilm height increased, but the rate of its increase slowed with time in the model (Fig. 5A). This behavior was in qualitative agreement with the results of the in vitro experiments (Fig. 5A). The coefficient of variation for the biofilm profile decreased with time, as the rate of decrease slowed down at later times (Fig. 5B). This behavior was also in qualitative agreement with results in the in vitro experiments (Fig. 5B).
Next we calculated C z (r,t) for the NTHI morphologies at different time points (16,40,64, and 88 h after the initial inoculation). Similar to the in vitro morphologies, C z (r,t) at low values (Յ5 m) of r, decreased as 1 Ϫ ar , where is Ͻ1 (Fig. 5C and D) for several z planes for biofilm growth. This demonstrated the presence of fractal structures in the interface between the NTHI-rich and -poor regions in the in silico model. We also (C) Variation of C z (r,t*)/C z (0,t*) with r at short distances (Յ5 m) for multiple z planes for the in silico biofilm at t* ϭ 88 h. The data were averaged over 30,000 initial configurations. The solid lines show the fits of equation 3 to the data. All the fits produces values of Ͻ1, in qualitative agreement with the in vitro data ( Fig. 2A and B). The dashed black line shows the ϭ 1 case as a reference. (D) Values of at different z planes obtained by fitting the C z (r,t*)/C z (0,t*) data for the in silico biofilm at t* ϭ 88 h. values are not shown where the C(r,t) decreased substantially in length scales of r Յ 1 m, comparable to the compartment size (l 0 ϭ 1 m). The values close to the biofilm surface also are not shown, due to large errors in the estimations due to the large errors in the data for C z (r,t*) near the biofilm surfaces. The symbol colors display the 2 values for the fits used to estimate the values. (E) Variation of C z* (r,t)/C z (0,t*) with r/ z* (t) at larger distances for multiple times (t* ϭ 40 h, 64 h, and 88 h) for the in silico biofilm at z* ϭ 12 m. The data show any absence of scaling, which is in agreement with the in vitro data (Fig. 3A). (F) Variation of the correlation length z (t*) (unit, m) with z at t* ϭ 88 h. The initial rise in z (t*) followed by a decrease as we approached the biofilm tip from the substrate was in qualitative agreement with the in vitro data (Fig. 3B, inset). calculated the correlation length z (t*) at a z plane by using C z (r,t*) for NTHI morphologies at a fixed time, t*. C z (r,t*) failed to scale with r/ z (t*) at different z planes at a fixed t (ϭt*) (Fig. S3), or C z* (r,t) failed to scale with r/ z* (t) at different values for different values of t at a fixed z plane (z ϭ z*) (Fig. 5E). This behavior was also in qualitative agreement with the NTHI profiles for biofilms formed in vitro. Similar to the morphologies formed in vitro (Fig. 3B, inset), z (t*) at a fixed t (ϭt*) first increased with z and then decreased as z approached the tip of the biofilm (Fig. 5F). There were few qualitative disagreements between the model and in vitro data; some of the disagreements arose due to the finite system size effects in the model. The early biofilm growth (16 h) in the model showed of Ͻ1 (Fig. S4), whereas the biofilm morphologies formed in vitro at a few z-planes generated values of Ͻ1 at t ϭ 16 h. This discrepancy could have arisen due to later (Ͼ16 h) induction of bacterial dispersion in the biofilms formed in vitro. Also, C z (r,t) becomes negative at larger values of r (Ͼ z (t)) in the model, whereas C z (r,t) when calculated from the confocal imaging data for the same values of r/ z (t) remains positive. This discrepancy between the model and experiments arises due to the influence of the finite size of the lateral (x-y) dimensions for the biofilms formed in vitro and in silico. The confocal images showed the presence of larger structures of sizes Ͼ40 m, comparable to the size (142 m by 142 m) of the biofilm area studied in the microscopy experiments. Therefore, the calculation of C z (r,t) at these larger length scales was influenced by effects due to finite sample sizes, e.g., incomplete statistical representation of spatial structures that were on the order of or larger than the biofilm sample. Resolution of this issue will require imaging of larger areas in the biofilm. The finite size of the system also affected the calculation of C z (r,t) in the in silico model at these length scales, which could be a possible reason behind the disagreement between the model and the data at larger length scales. Furthermore, many molecules, such as polysaccharides, that were not considered in the in silico model are likely to be present in the NTHI biofilms formed in vitro as a part of the EPS. The EPS can slow nutrient diffusion in in vitro culture (28) to produce nutrient gradients at larger length scales (Ն20 m). These effects could be responsible for the differences found between the biofilms formed in vitro and in silico at these length scales, e.g., with spikier patterns in biofilms formed in vitro compared to those formed in silico ( Fig. 4B and C).

Fractal structures in in vitro wild-type NTHI biofilms arise due to the presence of eDNA, as shown via model predictions with the ⌬comE mutant strain of NTHI.
We used the in silico model to determine mechanisms underlying the fractal structures present at short length scales (Յ5 m) in the biofilm morphologies. We found that in simulations where the eDNA network was limited and NTHI cells dispersed at a very low rate (~100 times lower than the WT) in all compartments, then the fractal nature of the interfaces in the biofilms was mostly absent. The absence of these fractal structures generated the short distance (Յ5 m) decay of C z (r,t) as 1 Ϫ ar , where the value is Ϸ1. In addition, the simulations under the above condition produced smoother biofilm surfaces (Fig. 6A) with smaller coefficients of variation than with the biofilms formed in silico by the wild-type NTHI. We compared the above predictions to experimentally generated values obtained from biofilms formed by NTHI strain ΔcomE, which potentially mimicked the above in silico conditions. Because the ComE pore is required for DNA release into the environment, biofilms formed by strain ΔcomE (Fig. 6B) contain very little eDNA compared to the parent strain (19). In addition, the ComE pore is also required for Tfp expression through the NTHI outer membrane. Thus, the ΔcomE mutant also lacks a fully functional Tfp, which could potentially give rise to a lower rate of dispersion than the WT strain. Biofilms formed in vitro were generated using the same protocol for biofilm seeding and maintenance as used for the wild-type NTHI strain. The predictions regarding the kinetics of the average height and the coefficients of variation for the biofilms generated using the in silico model for ΔcomE strain bacteria produced minimal amounts of eDNA, which was dispersed into the surrounding liquid medium at a lower rate and did not generate any twitching motility on eDNA strands in all the compartments, in qualitative agreement with the corresponding results for the in vitro ΔcomE mutant strain biofilms (Fig. 6C). In both the in silico model (Continued on next page) Structure and Kinetics of NTHI Biofilm Formation ® and the in vitro experiments, the older (e.g., 88 h) strain ΔcomE biofilms produced smoother biofilm surfaces (or lower coefficients of variation) (Fig. 6D). The calculation of C z (r,t*) for the NTHI densities for 88 h biofilms generated by the in silico simulations with the mutant NTHI strain generated bacterial clusters with smoother but sharper interfaces ( Fig. 6E and F). This was reflected in the values for the exponent, which for the majority of the z planes was close to 1 (Fig. 6F). This prediction agreed with the behavior of C z (r,t*) calculated for the strain ΔcomE biofilms formed in vitro at t* ϭ 88 h ( Fig. 6G and H). In the in silico simulations for the ΔcomE mutant strain, the lower rate of dispersion of bacteria as well as a substantial decrease of the motility of the bacterial cells along the eDNA network resulted in more homogeneous distributions of bacterial cells within the biofilms and smoother biofilm surfaces than with the wild-type NTHI biofilms formed in silico.

DISCUSSION
We quantitatively characterized the spatial organization of bacterial cells within NTHI biofilms formed in vitro by analyzing pair correlations between bacterial densities in confocal images of the biofilms. Our analysis revealed the presence of scale invariant or fractal structures at the interfaces separating the wild-type NTHI bacterial cells from the surrounding liquid. These interfaces occurred in shorter length scales (Յ10 m) ( Fig. 1E; Fig. S1). The presence of the fractal patterns can have several functional implications. The fractal interfaces increase the surface area of the NTHI biofilm, thereby exposing NTHI cells to the surrounding liquid, and this can facilitate greater absorption of nutrients from the environment. This could be particularly important for NTHI survival in the nutrient-poor environment of the middle ear. Furthermore, theoretical models investigating diffusion of particles on fractal landscapes (29,30) have shown the possibility of trapping the diffusing particles in a spatial location, because a fractal landscape can produce tortuous open spaces that can trap diffusing particles, like an "ant in the labyrinth" (31). Thus, the fractal interfaces in NTHI biofilms could trap molecules such as autoinducer-2 (AI-2), which is secreted by NTHI bacteria in a spatial location, and help mediate quorum signaling in that region. The localization of the AI-2 molecules could potentially generate spatial feedback for quorum sensing. However, the presence of the fractal surfaces also produces more regions with low bacterial cell densities, which could in turn lead to lower production of the quorum-sensing molecules. Thus, a trade-off between the above-described effects will determine how the fractal structures in a NTHI biofilm affect quorum sensing. While the effects described above can help in the survival of NTHI, the fractal structures can also increase the exposure of the biofilms to antimicrobial proteins, phagocytic cells, and other host immune effectors. Therefore, the beneficial or detrimental role of the presence of the fractal structures depends on how these opposing effects are balanced in a host environment. Our analysis also revealed the presence of multiple characteristic length and time scales in the biofilm morphologies (Fig. 3), implying the relevance of multiple dynamic processes (e.g., NTHI cell growth, NTHI cell mass movements and dispersion, formation of eDNA network) that span a range of length and time scales in regulating NTHI biofilm morphology. Therefore, altering specific features (e.g., changes in the fractal patterns) of NTHI biofilm architecture will likely require perturbing several processes simultaneously.  1 m). The values close to the biofilm surface also are not shown due to large errors in the estimation of C z (r,t*). (F) Values of at different z planes obtained by fitting the C z (r,t*)/C z (0,t*t) data for the in silico biofilm at t* ϭ 88 h. The 2 values, indicating the quality of the fits used to estimate the values, are shown by the color of the symbols. (G) Variation of C z (r,t*)/C z (0,t*) with r/ z (t*) at short distances (Յ5 m) for multiple z planes for the mutant ΔcomE biofilm at t* ϭ 88 h. Values of at different z planes were obtained by fitting the C z (r,t*)/C z (0,t*) data. The estimated values were in qualitative agreement with that for the model predictions in panels E and F. (H) Values of at different z planes obtained by fitting the C z (r,t*)/C z (0,t*) data for the mutant ΔcomE biofilm at t* ϭ 88 h. For the majority of the planes, was close to 1, as predicted by the in silico model (panel F). The 2 values indicating the quality of the fits used to estimate the values are shown by the colors of the symbols.

Das et al.
How unique are the fractal structures formed by NTHI biofilms? The existence of fractal structures within bacterial biofilms has been reported for a few biofilm morphologies, e.g., biofilms formed in vitro by activated sludge (32). The majority of the analyses of confocal images of bacterial biofilms in the recent literature were performed with the software package COMSTAT2 (33) and characterized the biofilms in terms of biomass, average height, and surface roughness. Thus, it is unclear whether fractal structures at short length scales are present in these studies. To address this question, we analyzed biofilms formed by another major OM pathogen, Moraxella catarrhalis, using C z (r,t) values calculated at multiple z planes. We found that for smaller length scales (Ͻ1.5 m), the exponent is close to unity for most of the z planes, which indicated a sharp nonfractal interface (Fig. 1D) between the M. catarrhalis cells and the surrounding medium (Fig. S5). Thus, fractal patterns at short length scales along the biofilm-liquid interface are not universal, and understanding how these structures arise in biofilms formed by different bacterial species under different experimental conditions is an interesting future direction of research.
We developed an in silico agent-based model to investigate the mechanistic roles of NTHI replication, eDNA, and Tfp expression in giving rise to fractal structures in NTHI biofilms. Recent studies have identified several candidates, including eDNA (8,9,11), quorum signals (7), and Tfp (13,14,34), as important regulators of the morphology of NTHI biofilms formed in vitro. The processes considered in the agent-based model were based on the known information about NTHI biofilm biology. We studied the roles of eDNA and Tfp in influencing the growth and structure of NTHI biofilms by explicitly considering potential eDNA-NTHI interactions in the in silico model. We hypothesized two different types of interactions between NTHI cells and the eDNA. (i) NTHI adheres to and ratchets along eDNA strands via twitching motility generated by Tfp. (ii) NTHI dispersion from the biofilm to surrounding liquid, induced by Tfp and quorum signals, occurs at a lower rate in the presence of the eDNA network. The growth and structure of the in silico biofilms agreed qualitatively with that in the NTHI biofilms formed in vitro. The role of the above interactions between the NTHI and eDNA in giving rise to the fractal patterns in wild-type NTHI biofilms was investigated further by analyzing spatial patterns in in silico biofilms for an NTHI mutant that produced substantially less eDNA within the biofilm matrix and showed a low dispersion rate and impaired motility along the eDNA network. The biofilm morphologies for the in silico mutant were compared with that of biofilms formed in vitro by a ΔcomE mutant of NTHI. The biofilms formed by the ΔcomE mutant contained substantially less eDNA than the wild type, and the Tfp function was impaired in the ΔcomE strain. As predicted by the in silico simulation, the ⌬comE strain biofilms displayed a considerable loss of the fractal patterns found in the wild-type biofilms and smoother biofilm surfaces (Fig. 6). The biofilms formed by the pilA mutant of NTHI, in which NTHI cells are devoid of Tfp function and thus are unlikely to possess the NTHI-eDNA interactions hypothesized above, also showed similar loss of the fractal patterns in the biofilm morphologies (Fig. S6). These findings support the hypothesized interactions between NTHI and eDNA in biofilms.
We used an agent-based model to describe the kinetics of NTHI biofilm formation. To the best of our knowledge, this is the first in silico model that describes the formation of NTHI biofilms. The agent-based model described the biofilm formation in terms of agents, such as NTHI bacterial cells, eDNA strands, and nutrients, that interact with each other with specific rules representing relevant biological processes, e.g., NTHI replication. An advantage of agent-based modeling is the flexibility to introduce rules that can be used as hypotheses to probe specific interactions between the components (e.g., NTHI bacteria and eDNA) of interest in a biological system. Agent-based models have been widely used in modeling generic biofilm growth (24,(35)(36)(37) as well as in a wide range of host-pathogen systems (38) where spatial heterogeneity plays an important role. However, these biofilm models did not explicitly investigate the role of eDNA in regulating biofilm morphologies. Biofilm models based on other methods, such as partial differential equations (39)(40)(41) and Brownian dynamics (42), have also been used previously to describe generic biofilm formation or biofilm formation in specific bacterial species, such as P. aeruginosa (41) and Bacillus subtilis (40,43). These models usually entail microscopic details pertaining to specific bacterial species or contain simplifying assumptions; therefore, they cannot be readily used for describing NTHI biofilm formation.
How accurate is the in silico model in describing biofilm formation by other bacterial species? The agent-based model includes distinct hypotheses pertaining to NTHI biofilm formation, such as Tfp-induced movements on the eDNA network and the influence of Tfp and eDNA on NTHI dispersion. Though eDNA is found in the EPS of biofilms formed by many other bacterial species, such as M. catarrhalis, and many other bacterial species (e.g., P. aeruginosa) express Tfp, it is unclear if the specific hypotheses regarding the interactions between eDNA and Tfp in NTHI biofilms can be generalized to other bacterial species. Therefore, based on this study's findings, we anticipate that the in silico model can be used to analyze biofilms formed by the wild-type or mutant NTHI strains but the model as it stands cannot be readily applied to describe biofilm formation in another bacterial species, such as M. catarrhalis, S. aureus, or P. aeruginosa. However, the in silico model still contains the flexibility of an agent-based modeling framework (38), where the existing rules in a model can be modified and new rules can be added easily; therefore, the current model can be extended to describe biofilm formation by other species.
What do the results obtained here mean for biofilms formed in vivo? Since the middle ear presents NTHI with a nutrient-poor and hostile environment, we expect that the survival of NTHI would require formation of these fractal structures. The in silico analysis and the agent-based model can be easily extended to analyze NTHI biofilm formation in vivo. The irregular three-dimensional shape and larger size (~100 mm) (44) of the middle ear cavity in Chinchilla lanigera, an animal model for experimental OM, will likely play an important role in determining formation and kinetics of NTHI biofilms in vivo. Therefore, an extension of the current in silico model will require simulations of biofilm growth in three dimensions in larger (Ͼ10 times) system sizes and should also incorporate explicit nutrient diffusion, as the current model assumption of instantaneous equilibration of the nutrient density will not be appropriate at larger distances. Another useful improvement for the in silico model will be to include a more detailed model for eDNA network formation based on microscopic details of the eDNA network informed by in vivo imaging experiments (11). Such an effort is likely to offer exciting mechanistic insights regarding the role of the eDNA network and Tfp expression in regulating NTHI biofilm structures in the middle ear.

Analysis of in vitro NTHI biofilm morphology based on the pair correlation function.
We calculated the pair correlation function between the bacterial densities at two locations, r → 1 ϵ͑x, y, z͒ and r → 2 ϵ͑x ϩ r x , y ϩ r y , z͒ in a plane (x-y) parallel to and residing at a fixed distance (z) from the substrate (Fig. 1A). The pair correlation function C z (r,t) is widely used in statistical physics (22,23,26) and materials science (21) for characterizing spatial structures. C z (r,t) is defined as shown in equation 1: where r ϭ | r → 1 Ϫ r → 2 | ϭ ͱ ͑r x 2 ϩr y 2 ͒ denotes the distance between the two locations in the x-y plane and L is the length of the sample along the x or y direction. In the above expression, it has been assumed that the bacterial density in the x-y plane depends on the magnitude of the relative distance between the vectors r → 1 and r → 2 (translational symmetry) and does not depend on the direction of the vector ͑ r → 1 Ϫ r → 2 ͒(rotational invariance). The angular brackets indicate an average over an ensemble of density configurations. The behavior of C z (r,t) at short distances or small values of r can be used to determine the nature of the spatial patterns within that length scale. When r is larger than a microscopic length scale a 0 (e.g., spatial resolution of the confocal images or length scales Ϸ0.27 m) but smaller than an intermediate scale w, i.e., a 0 Յ r Յ w, C z (r,t) can be approximated by using equation 2 (21,22): The value of the exponent contains information regarding the nature of the spatial pattern in that length scale (21)(22)(23)26). For example, if the bacterial cluster, embedded in a space of dimension d, forms a scale-invariant structure or a mass fractal (Fig. 1B) When the bacterial cluster forms a regular solid structure associated with a fractal (Fig. 1C) or a sharp (Fig. 1D) interface separating regions rich and poor in bacterial density, then the equation ϭ d Ϫ d s , and ϭ d Ϫ (d Ϫ 1) (Porod's law [22,27]), respectively. d s (d Ͼ d s Ͼ d Ϫ 1) denotes the dimension of the fractal surface. In our analysis, d ϭ 2, therefore, of Ͻ1 will indicate the presence of bacterial density in a surface fractal pattern. If ϭ 1, it will indicate the presence of a sharp but nonfractal interface, and a value of Ͻ0 will indicate the presence of bacterial clusters with mass fractal structures. We found that the in vitro NTHI biofilm morphologies produced different types of surface fractals at different length scales that were characterized by different values of the exponent, i.e., the surface fractal crosses over from one type to another type as the distance r increases beyond the length scale r c (see Fig. S1) (22). This behavior can be described in terms of C z (r,t) as shown in equation 3: when r Ͻ r c , and when r Ͼ r c , equation 3 is as follows: The changes in equation 3 are well described by a crossover function (22) (Fig.   S1), as shown in equation 4: where g(r,r c ) ϭ r m /(r m ϩ r c m ) and m ϭ 5. We also attempted to fit C z (r,t*) for t* ϭ 88 h for a larger range of r values (0 to 69.6 m) via a function described by a combination of different power law functions in three segments: 0 Յ r Յ r 1 , r 1 Ͻ r Յ r 2 , and r 2 Ͻ r Յ 69.6 m. The fitting function was constructed by considering crossover behavior similar to that in equation 4, resulting in equation 5: Ϫ g(r, r 1 )͔a 1 r 1 Ϫ ͓g(r, r 1 ) Ϫ g(r, r 2 )͔a 2 r 2 ϭ g(r, r 2 )a 3 r 3 (5) where g(r,r 1 ) ϭ r m /(r m ϩ r i m ), i ϭ {1,2}, r 1 Ͼ r 2 , and m ϭ 5. Equation 5 fit reasonably well to C z (r,t* ϭ 88 h) ( Fig. 3B; Fig. S1B). This suggests that the spatial patterns in the biofilms formed in vitro by wild-type NTHI possess different fractal structures extending to large (r Ͼ 10 m) length scales.

Agent-based model simulations. (i) Kinetic Monte Carlo simulation.
We implemented a kinetic Monte Carlo (kMC) scheme to simulate NTHI biofilm growth, using the rules described in Table 1. kMC methods have been used in simulating stochastic kinetics in diverse spatially resolved systems. We followed the algorithm described in reference 25. Each MC step (MCS) involved carrying out N X ϫ N Z MC trials, where N X ϫ N Z is the total number of compartments considered in the model. Execution of an MCS evolved the system by a time step of Δt/n, where n (ϭ 9) is the total number of rules. At each MC trial, the following steps were implemented. All the random numbers used in the following steps were drawn from a uniform distribution with the specified range.
(1) Draw two random integers, 1 Յ x Յ N X , and, 1 Յ z Յ N Z , to choose the compartment at the location ( x , z ) for execution of a rule to be decided in the next step.
(2) Call a random number 1 between 0 and n. When i Ϫ 1 Ͻ 1 Յ i (i Յ n), the process given by rule no. i is implemented with rate k i by following step 3.
(3) Draw a random number 2 between 0 and 1, and if 2 Յ p i , then process or rule 1 is implemented. p i is the probability for the ith process to occur and is given by the equation k i ϫ Δt, where k i is the rate for process i.
(ii) Boundary conditions. We imposed periodic boundary conditions at x ϭ 0 and x ϭ N x and no-flux fixed boundary conditions at z ϭ 0 and z ϭ N z .
Determination of the biofilm surface and calculation of h avg and the coefficient of variation. (i) Confocal image data. The two-dimensional biofilm surface z(x,y) was determined by approaching the biofilm from the top along the z axis and identifying the voxel located at {x, y, z(x,y)} where the intensity crossed a threshold for the first time. The threshold intensity was taken to be 1/10 the highest voxel intensity in an image. The average height h avg in the confocal images was calculated using the formula h avg ϭ 1/N voxel Α i z i (x i ,y i ), where z i (x i ,y i ) denotes the biofilm height at voxel i on the biofilm surface and N voxel is the total number of voxels in the computed biofilm surface. The standard deviation is calculated using the equation 2 ϭ 1/N voxel Α i [z i (x i ,y i )] 2 Ϫ (h avg ) 2 and the coefficient of variation is calculated as /h avg .
(ii) In silico model. The same procedure as that used for the confocal imaging data was followed for the in silico model. In this case, the biofilm surface, z(x), resides in one dimension. The threshold is chosen to be 1 NTHI bacterial cell per in silico compartment.
Evaluation of the exponent. C(r,t) versus r data for r values of Յ1.5 m or r/(t) values of Ͻ1 were fitted with a function, 1 Ϫ ar , where the fitting parameters a and were determined using the Levenberg-Marquardt algorithm (45). The calculations were performed using MatLab.
Formation of NTHI biofilms in vitro. Several microbial in vitro biofilm models have been developed in order to study biofilm formation, structure, and physiology under specific laboratory conditions (1,2). We generated biofilms formed in vitro by nontypeable Haemophilus influenzae strain 86-028NP or its isogenic ΔpilA or ΔcomE mutants on eight-well chambered coverglasses under static conditions according to a method described previously by Jurcisek et al. (51). Briefly, bacteria from an overnight chocolate agar culture were suspended in brain heart infusion broth that had been supplemented with 2 g/ml each of heme (Sigma-Aldrich, St. Louis, MO) and ␤-nicotinamide adenine dinucleotide (Fisher Scientific, Pittsburgh, PA) (sBHI) to an optical density at 490 nm of 0.65. The bacterial suspension was diluted 1:6 in medium and incubated for 3 h at 37°C, 5% CO 2 , under static conditions. The culture was further diluted 1:2,500 with sBHI, and 200 l of the diluted bacterial suspension, equivalent to 4e4 CFU NTHI, were inoculated into each well of an eight-well chambered coverglass (Fisher Scientific). Biofilms were grown under static conditions at 37°C, 5% CO 2 in a humidified environment. The medium was replaced with fresh sBHI at intervals of 16 h (night) and 8 h (day). Biofilms were processed for imaging at 16, 40, 64, and 88 h postseeding.
Biofilm staining and imaging. At the final time points, biofilms were washed twice with sterile saline solution and stained using the Live/Dead BacLight bacterial viability kit (Molecular Probes, Eugene, OR) according to the manufacturer's instructions. Biofilms were then fixed as described previously (8). z-stack images were collected with a Zeiss 510 confocal laser scanning microscope (CLSM; Carl Zeiss, Thornwood, NY) with a 63ϫ, 1.2 numerical aperture objective. All biofilm assays were performed in duplicate for a minimum of three times.