Biofilms are surface-attached, crowded communities of microbes that can grow in two directions: horizontally across the surface or vertically above the surface [1–3]. While horizontal growth and range expansion has been more extensively studied [4–7], a growing body of evidence suggests that vertical growth is a crucial aspect of bacterial behavior [8–12]. It has even been shown that vertical growth can impact horizontal growth, implying that colony expansion can only be fully understood by considering both phenomena [13]. Additionally, in some environments, confinement can prevent horizontal growth while allowing vertical growth to continue. Therefore, understanding biofilm growth requires a comprehensive examination of vertical growth dynamics.

Our understanding of vertical growth has been limited in part by its complexity. Nutrients diffuse into the colony through its interfaces, and then are uptaken by microbes. Thus, in a single vertical direction, bacteria experience many different microenvironments [14, 15]. There are many approaches to modeling such systems, cellular automata-based models [16–18], individual based models [4, 19–21], or to treat the microbes as an active fluid. Models can even mix aspects of these approaches [22, 23]. In active fluid models, microbes are treated as growing and energy-consuming particles characterized by viscous interactions, allowing for analytical modeling of biofilms based on fluid principles [24, 25]. Prior research has modeled the expansion of biofilms via nutrient diffusion, uptake, and growth as an active fluid [26–32].

A recent paper, Bravo et al [15], proposed and experimentally validated a simple heuristic model which we will refer to as the interface model (see figure 1):

Figure 1. (A) Experimental data from Bravo et al [15] depicting the growth of A. veronii colony height over time, measured using an interferometer. The interface model was applied using equation (1) and the active fluid model using equation (27). Fit parameters for the active fluid model are , , and ; for the interface model, , , and .(B) This panel shows the slope of the height over time curve as a function of height, with experimental data again sourced from Bravo et al [15]. The active fluid model is represented by the green line and the interface model by the orange line, as in panel A.

Download figure:

While this model captured the empirical behavior of nine different species, including prokaryotes and eukaryotes, gram negative and gram positive bacteria, and cells with a wide range of shapes and sizes, it was proposed based on empirical data, and justified term-by-term, rather than derived from fundamental dynamics. Thus, the underlying biophysics behind vertical growth remains unclear.

Further, the interface model was unable to explain multiple empirical observations. In particular, the vertical growth rate consistently reached a maximum at biofilm heights greater than the characteristic growing region of length L. This effect has also been observed at the edge of growing colonies [13]. However, without a pure theoretical foundation on which to develop hypotheses, it is challenging to investigate the deviations between the model and empirical observations.

Here, we show that the interface model from Bravo et al [15], arises directly from a commonly used active fluid model of biofilm growth. We derive the interface model from analytical expressions of the active fluid model, and then we fit the experimental data from Bravo, et al to both the active fluid model and the interface model, showing that the heuristic model is nearly as good (see figure 1). Moreover, we identify that the discrepancies between the interface model and the observed maximum growth rates are effectively captured by the active fluid model. Thus, we conclude that vertical growth can be accurately modeled using the simpler heuristic framework proposed by Bravo et al [15], with minimal deviation.

Here, we develop a one-dimensional model for vertical biofilm growth, with the goal of producing a model for the height over time in alignment with experimental observations. We initiate this process by deriving a model for a growing ‘bio-fluid’ based on the principles of active fluid dynamics. The derivation here closely follows the methodologies established in previous works [22, 27, 28, 30, 33]. We consider a biofilm, consisting of a single species of cells, growing vertically upward. We let be the cell mass density at a distance, z, above the base of the biofilm at time t. We also consider the mass concentration of a generic diffusible resource, , (e.g. a nutrient) which is consumed by the cells with rate . The mass continuity equation is then:

where u is the velocity field of the system and g is the growth function; this will be a function the resources R. Note, diffusion of the bio-fluids is assumed to be negligible [24]. The resources diffuse, move under the same velocity field as the bio-fluids, and are uptaken by each mass species in the environment:

where DR is the diffusion coefficient and there is some functional form of consumption for the system’s mass species. The resources themselves are small molecules and thus we do not model them as having volume in the same way as the bio-fluid.

As the bio-fluid can be treated as incompressible (they are mostly made of water [34]), we can re-frame the equations away from mass density, instead replacing them with volume fraction where ρ is the specific density of the cells themselves. This leaves N as the volume occupied by our species divided by the total volume at a given position. We also replace the growth function with volume growth function :

Note this one dimensional, incompressible, single species approach also allows for a ‘no voids’ constraint, meaning that at every position within the biofilm, the volume fraction of ‘bio-fluid’ must fill all available space:

We can now take the volume density continuity equation (4) and apply the constraint equation (5) to arrive at an expression relating the velocity field and growth:

It is next assumed that this bio-fluid can be treated as homogeneous, viscous, and incompressible with constant density flowing through a porous environment, satisfying Darcy’s law; this allows us to model the global velocity field with a scalar potential, i.e. the pressure field p:

Finally, the growth function must be defined. The terms in the growth function indicate the increase or decrease of the bio-fluid volume; hence, a positive term signifies volume creation, while a negative term indicates volume loss. Our species grows in proportion to both its concentration and the availability of resources at a specific location (R; it also experiences decay at rate β. Thus, our growth function is:

where equation (5) allows us to replace N with 1 (equation (8)). Note, the decay rate β signifies the rate at which the bio-fluid volume decreases, not the disappearance of the mass itself, but rather the reduction in volume occupied by the bio-mass. We assume that cells degrade on a timescale that is short compared to the timescale of growth; hence degradation is treated as instantaneous. Additionally, while alternative common growth functional forms related to resource consumption, such as the Monod equation, could certainly be utilized, the form we have adopted facilitates the solution of the system. The use of an unbounded maximum growth rate as R increases could potentially introduce complications; however, we will implement a Dirichlet boundary condition on the resources at on the bottom layer of the colony (see equation (11)). is then a prescribed concentration of resources at the bottom layer, stabilizing R to be at most (so long as the colony does not drop suddenly in height leading to a compressed R level). Thus, while growth is typically defined explicitly by a predetermined maximum growth rate α, we can still impose such a maximum by aligning it with realistic values through the use of constant c , such that . Thus, we replace c with α.

Now, we must define the resource consumption term which is now a function of N not m as we re-framed away from mass density. As seen above, the rate at which resources are consumed and made into the bio-fluid is controlled by the constant c, which has units of one over mass-density time. This leads to the following consumption term where ε is a yield constant, with units of volume fraction of bio-fluid produced per concentration of resources consumed [35]:

To continue towards our goal of generating an expression for the height of the colony over time, we next define the rate at which the height changes over time. As it is solely a function of the velocity u at the top of the colony, we can write it as such:

We now define our boundary conditions. The gradient of pressure must be zero at z = 0 (which imposes that the biomass velocity field must be zero at z = 0, thus the biomass cannot flow into the ground). We imposed a Dirichlet boundary (p = 0) at the top of the colony (z = h) in line with previous work [24, 33] as pressure is only unique up to an additive constant. For the resources, we assume that the agar underneath the colony acts as source (at z = 0) such that its concentration is constant (see equation (11)), and the gradient of R is zero at the top of the colony (see equation (12)). Thus, resources flow freely from the bottom of the colony, but they cannot flux out of the top of it:

At this point we note that we have not defined any initial conditions for our system. This is intentional as our goal is to find an equation for the height of the colony over time which we will then fit to the experimental data from [15] in the same manor as the interface model.

We now derive a solution for the height of the colony over time from the active fluid model from section 2. To do so, we must find the velocity field at the top of the colony (equation (10)), which means we must solve for the pressure gradient (equation (7)), which means we must know the solution of the growth functions (equation (6)), which means we must find the solution for R (equations (17) and (9)). This expression has a time dependency and the unknown velocity field u. To address this, we nondimensionalise the system with the following scalings:

Where the advection time is based on cell-movement and thus the growth rate (), and h0 is the initial height of the colony (≈ 1 cell length). This also allows the introduction of the following dimensionless parameters

where is the growth rate scaled by the diffusion time and Pe is the Péclet number (the ratio of diffusion time to advection time) This leaves us with the following expression for the time dependent resource equation:

Using values from [15] Pe is , implying that diffusion is significantly faster than advection and growth in our system. This ‘rapid’ diffusion of resources is an important characteristic of such systems, and it allows us to neglect the advection and time dependent terms from (17). This leaves us with the quasic static approximation in which we return to dimensional units:

And the general solution:

Using the boundary conditions defined in equations (11) and (12), we arrive at the following expressions for our unknown constants. We find a characteristic length scale L (equation (20)); it is proportional to the diffusion of resources, their concentration, and inversely to the growth rate. The functional form of this length scale is similar to characteristic lengths in horizontal growth models [36]. It will become clear later that this L is the same L as in the interface model:

Thus, as h grows, the contribution of the positive exponential term decays; however it is significant in early times when h is low, see figure 2(B). The combination of these terms results in higher concentrations of resources when the height is low; this is illustrated in figure 2(A) and show in figure 5(A). It is now more convenient to proceed with a hyperbolic solution rather than exponential; thus, our resource curve is the following:

Figure 2. (A) When the total height of the colony is small, the no-flux top surface of the biofilm allows for more growth (orange cells) than when the colony is large, and the cells beyond L can no longer grow (blue cells). (B) The magnitude of the coefficient C1 is plotted over height of the colony where is the coefficient of the exponentially growing term in equation (19) while C2 is the coefficient of the exponentially decaying term. Thus, we see that at small heights C1 has some impact, while C1 has little-to-no impact at larger heights. (C) A schematic of our system outlining the coordinate transformation to z′, growth rate α, death rate β, and diffusion of nutrients.

Download figure:

It is further convenient to undergo a coordinate transformation according to equation (24). This does not change any of the results we already have derived so far, but it will change our derivative functions going forward according to the equations below. This transformation allows us to rigorously define the domain . Similar coordinate transformations were used in [27, 28, 33]:

With this coordinate transformation we can express the growth of our bio-fluid on a defined domain, see equation (25):

Now we can combine our expression for the growth of the bio-fluid with our expressions for velocity and pressure (equations (6) and (7)) which yields the following expression for the second derivative of pressure (see equation (26)). The extra factor of h2 appears from the coordinate transformation:

As the pressure is in not a function of itself and only spatially defined by z′, we can integrate across the whole domain, and then use the boundary conditions for pressure to find the pressure gradient at the top of the colony. We then use this pressure gradient and equation (7) in the modified coordinates (which cancels out another factor of h) to find the velocity at the top of the colony. Finally, with equation (10), we have arrived at an expression for the rate that the height changes over time:

Thus, we recover the interface model described in Bravo et al. The active growth regime of length L is the characteristic length of diffusing resources as suggested in the study. At , we are left with:

and at :

Using the expression in (27) allows us to find the height at which the growth rate is maximized by taking the derivative with respect to h and setting the expression to zero:

We now compare the interface model (as outlined in Bravo et al [15]) and the full analytical expression (equation (27)), referred to as the active fluid model, to experimental data (originally published in [15]). Figure 1(A) illustrates the best fits of both the interface and active fluid models against the height-over-time curve of an A. veronii colony (for this and all fitting in this paper, we use the first height recorded as the initial condition) , revealing that the two models are virtually indistinguishable. Furthermore, figure 1(B) highlights the similarity in form between the fluid and interface models, though the rationale behind the active fluid model’s form is more transparent, and it improves upon the interface model’s insights. To comprehensively demonstrate that the active fluid model is as robust as the interface model, we fitted the height-over-time curves of all nine strains outlined in [15], as presented in figure 3. As we observed previously, the fits between the two models are indistinguishable. This similarity is expected, given that both models are grounded in the same fundamental forms and share the same number of parameters. We also find that the fit parameters are all very similar to each other as seen in figure 1(A). This makes sense for α and β, as they represent growth and death, and both models treat growth and death similarly; we discuss L further below.

Figure 3. Displayed is the active fluid model fit to the same experimental data from Bravo et al [15] using equation (27). The fit is as precise as the interface model detailed in (1) with similar differences in fit parameters as shown in figure 1.

Download figure:

Our objective is not to present a model that better fits the experimental data, but rather to validate that the broad functional form adopted is appropriate and captures a previously unaddressed observation. In particular, the interface model proposes a characteristic height, L, which marks the transition from the actively growing regime to the nutrient-depleted regime, a transition from mostly growth to mostly decay. This suggests that in every cellular layer above height L, the rate of decay exceeds the rate of growth. Thus, height L should represent the level at which biofilm growth reaches its maximum rate. However, experiments consistently show that the growth rate peaks (denoted ) at a biofilm height ; specifically, as illustrated in Figures 1(B) and 4(A), L does not align with .

Figure 4. (A) The height at which is maximized, as calculated from equation (30), is plotted in green against the experimental peak in versus h from all data shown in figure 3. In orange, the L fit from the interface model is plotted against the same experimental peaks. (B) The residuals, representing the distance from the experimental peak in A and each model’s , are displayed for each strain to illustrate how the active fluid model more accurately captures this phenomena in all but one strain.

Download figure:

The active fluid model provides a framework to resolve this discrepancy. The best fit L values in both models are notably similar, which can be understood through the work presented in the previous section which shows that they measure the same phenomena. While L represents the characteristic length at which resources are depleted in both models, equation (28) shows that L only represents the size of the growing region when the colony achieves sufficient height (); if , then the thickness of the active growing region is greater than L. This phenomenon arises from the boundary condition on the upper surface of the colony where nutrients cannot flux outside of the biofilm; this prevents a simple characteristic decay of nutrients, as demonstrated by equations (21) and (22). The nutrients are ‘reflected’ back, so at low h the cells have an abundance of nutrients (see figures 5(A) and 2(B)). Thus, as the biofilm grows, the top surface becomes farther from the nutrient source (with increasing h), the taller biofilm develops a distinct characteristic depth L of metabolically active cells. This phenomenon is qualitatively illustrated in figure 2(A). Figure 5(A) demonstrates that when the colony initially attains height L, the nutrient concentration at z = L can be approximately 70% higher due to the boundary effect compared to later where h exceeds L, and the two converge. This interaction results in a more gradual transition in our final expression, and it facilitates the theoretical determination of through equation (30). Experimentally measured values are much more aligned with the active fluid model predictions than with the interface model L (figures 4(A) and (B)).

Figure 5. (A) The available nutrient concentration is shown at three different z distances as the height of the colony grows; we see that as the colony grows the nutrient concentration decreases. The black dashed line represents the steady state concentration at z = L, so we see that at low h the growing region is greater than L. (B) The ratio , where is the experimental height at which is maximized and L is the fit parameter from equation (27), is plotted against the fit parameters as green circles. The dashed line represents equation (30). We observe that the strains from Bravo et al [15] occupy a closely aligned regime within the parameter space, and the active fluid model closely approximates the predicted .

Download figure:

An alternative explanation could be reached from figure 4(A) which demonstrates a nearly linear positive correlation between L and experimental , this relationship might suggest there is residual growth beyond the characteristic nutrient decay length L, consistent with the persistent nature of exponential decay curves. In essence, L and may only be a factor apart, a