2.1 Classic Keller-Segel model, also incorporating growth

Two forms of the Keller-Segel model have been explored in the prior literature to describe two distinct biophysical processes: one describes cellular aggregation and pattern formation in response to chemoattractant produced by the cells themselves [44–46], while the other describes cellular spreading in response to an exogenous chemoattractant that is not produced, but just consumed, by the cells [26, 27, 29–33]. Here, we focus on the latter case. Before considering confinement, we first describe how chemotactic spreading is typically modeled using this form of the classic one-dimensional Keller-Segel model—which does not incorporate the influences of growth and confinement, but can successfully capture the key features of experiments on dilute populations of bacteria in bulk liquid [25–27, 29–34]. We also introduce growth into this model.

To directly connect the model to many experiments [24, 26, 28, 34], we consider a sole nutrient that also acts as the chemoattractant—as is conventionally done [24, 29, 32]—represented by the continuum variable c(x, t), where x is the position coordinate and t is time. The number density of bacteria, in turn, is given by the continuum variable b(x, t). Furthermore, given the experimental conditions, we assume that the cells do not excrete their own chemoattractant or other diffusible signals, as is sometimes the case in low-nutrient conditions and for specific strains. Recent extensions of this model have also considered the case in which nutrient and attractant are separate chemical species, which leads to fundamentally different behavior that would be interesting to explore using our framework in future work [27, 47].

As the nutrient diffuses through space with thermal diffusivity D_c, it is consumed by the cells at a rate bκg(c); here, κ is the maximum consumption rate per cell and the Monod function g(c) ≡ c/(c + c_char), with the characteristic concentration c_char, quantifies the reduction in consumption rate when nutrient is sparse [27, 34, 48–52]. Therefore, the nutrient dynamics are given by (1)

The bacterial dynamics have two contributions: a motility-driven flux and cellular proliferation. The flux arises from the combination of the undirected spreading of cells, a diffusive process with an active diffusivity D_b0 [53], and directed chemotaxis with a drift velocity that quantifies the ability of the bacteria to logarithmically respond to the local nutrient gradient [29–31]. The well-established function f(c) ≡ log [(1 + c/c₋)/(1 + c/c₊)] quantifies the ability of the cells to sense nutrient with characteristic bounds c₋ and c₊ [26, 27, 54–62], while the chemotactic coefficient χ₀ quantifies the ability of the cells to bias their motion in response to a sensed nutrient gradient. Therefore, the motility-driven flux . Proliferation, on the other hand, is given by bγg(c), where γ is the maximal growth rate per cell and g(c) reflects the reduction in growth rate when nutrient is sparse—circumventing the need to introduce an ad hoc “carrying capacity” of a logistically-growing population, as is sometimes done. However, as we detail further in §2.3, the added feature of confinement does introduce a maximum cell density at which purely growth-driven spreading dominates. Here, in the absence of confinement, the bacterial dynamics are therefore given by (2)

Together, Eqs 1 and 2 represent the classic Keller-Segel model that describes the coupled dynamics of nutrient and bacteria, also including the added influence of cellular growth. In particular, they successfully capture the key features of chemotactic spreading in unconfined liquid, in which cells collectively generate a local gradient of nutrient that they in turn bias their motion along—leading to the formation of a coherent pulse of bacteria that continually propagates, sustained by its continued consumption of the surrounding attractant [23–26].

2.2 Characteristic dimensionless parameters

Non-dimensionalizing Eqs 1 and 2 yields useful dimensionless parameters for characterizing population spreading. We rescale {c, b, t, x} by the characteristic quantities {c_∞, b₀, t_c,0, ζ}, where c_∞ is the initial nutrient concentration taken to be constant everywhere, b₀ is the maximal initial cell density, t_c,0 ≡ c_∞/(b₀κ) is a characteristic time scale of nutrient consumption, and is the characteristic extent of cellular diffusion over the duration t_c,0. This process yields the non-dimensional equations (3) (4) where the tildes indicate non-dimensionalized variables. Three dimensionless parameters emerge:

• The diffusion parameter δ₀ ≡ D_c/D_b0 compares the thermal diffusion of nutrient to the active diffusion of bacteria. When δ₀ ≪ 1, variations in nutrient are localized to the leading edge of the bacterial population, whereas when δ₀ ≫ 1, nutrient levels vary over large spatial extents.
• The directedness parameter α₀ ≡ χ₀/D_b0 compares the influence of chemotaxis to active diffusion in driving cellular spreading. When α₀ ≪ 1, diffusion dominates and cells do not appreciably direct their motion in response to a nutrient gradient, whereas when α₀ ≫ 1, motile cells strongly direct their motion in response to a gradient.
• The yield parameter β₀ ≡ γ/(b₀κ/c_∞) compares the rates of cell growth and nutrient consumption, γ and , respectively. It therefore quantifies the yield of new cells produced as a population consumes nutrient. When β₀ ≪ 1, nutrient consumption is much faster than proliferation, whereas when β₀ ≫ 1, many new cells are produced in the time required to consume the available nutrient.

The quantity Λ₀ ≡ α₀/(β₀δ₀) = γ⁻¹ ⋅ χ₀/(D_ct_c,0) therefore characterizes the interplay between chemotatic and growth-driven spreading of bacterial populations. In particular, [χ₀/(D_ct_c,0)]⁻¹ is a characteristic time scale needed to spread via chemotaxis over the nutrient diffusion length , while γ⁻¹ is the time scale over which cells grow. Previous studies in bulk liquid focused solely on chemotactic spreading, which is characterized by the limit Λ₀ ≫ 1 [25, 26, 29, 31, 33]. Other studies of non-chemotactic cells focused solely on growth-driven spreading, characterized by the opposite limit Λ₀ = 0 [41, 42, 63–66]. However, experiments performed in semi-solid agar [34] as well as in defined packings of particles [28] indicate that confinement in such crowded media introduces new cell-cell and cell-medium interactions that are not incorporated in the classic Keller-Segel model. Hence, in this paper, we describe a first step toward incorporating these complexities, which not only tune Λ₀ over a broad range, but also fundamentally alter spreading dynamics—as described hereafter.

2.3 Keller-Segel model incorporating confinement

As a model system, we consider bacterial populations confined in media with close-packed, rigid, and immovable obstacles surrounding a free space that is sufficiently large for cells to move through. This form of confinement alters bacterial spreading dynamics in three ways:

1. (i). Collisions with the surroundings impede cellular spreading [35–38], reducing the transport parameters D_b0 and χ₀, as quantified in recent experiments in 3D porous media [28, 38] as well as in semi-solid agar [34];
2. (ii). The presence of surrounding obstacles reduces the free space available to cells to move through, increasing cellular crowding and promoting cell-cell collisions that further truncate the motility parameters, observed experimentally using in situ microscopy [28];
3. (iii). When the number density of cells is sufficiently high, this reduction in free space causes the cells to be jammed; hence, they are able to spread only through proliferation, which pushes cells outward, as quantified in experiments using single cell visualization [39, 40].

Notably, (ii)-(iii) are absent from the classic Keller-Segel model, which treats cells as non-contacting, and require modifications beyond simply changing the transport parameters D_b0 and χ₀.

1. (i). Impeded spreading of isolated cells. Bacterial spreading is typically modeled as a random walk with directed steps of characteristic length ℓ and characteristic duration τ that are punctuated by reorientation events [53]. Consequently, both transport parameters D_b0, which describes the unbiased component of the random walk, and χ₀, which describes the biased component, are set by ∼ ℓ²/τ. In bulk liquid, the directed steps are known as runs, which extend along straight-line paths ℓ ∼ 40 μm long, punctuated by rapid tumbles. In tight confinement, however, a cell collides with an obstacle and becomes transiently trapped well before it completes such a run. Therefore, as established in recent experiments [37, 38], runs are truncated by collisions with surrounding obstacles, and the directed steps of the random walk are instead set by the geometry of the available free space; thus, for isolated cells, ℓ ∼ ℓ_c, the mean length of straight line chords [67] that fit in the free space [37]. Moreover, because the trapping process induced by collisions with obstacles occurs over a duration τ_t that is longer than that of the truncated runs, τ ≈ τ_t. As a result, for cells confined in tight media, both transport parameters D_b0 and χ₀ are instead —and because increasing confinement both decreases ℓ_c and increases τ_t [37], it concomitantly decreases both D_b0 and χ₀, as confirmed experimentally [28]. Within the context of prior work investigating diffusion in porous media [68], we note that while the volume fraction of free space (porosity) ϕ is known to influence diffusion, it alone does not determine the diffusion coefficient because the geometry of the free space plays a key role as well. Thus, in our model, ϕ influences cellular transport (active diffusion and chemotaxis) indirectly through its effect on the chord length ℓ_c, which characterizes the length scale associated with straight paths that fit within the free space—and therefore determines the length scale over which cells can move in a directed manner.
2. (ii). Crowding-induced collisions between cells. Confinement also reduces the free space available to cells. Our definition of the number density of bacteria b quantifies the number of cells per unit total volume of space, which includes the volume of surrounding obstacles; hence, the local density of cells is given by b/ϕ, where ϕ < 1 is the volume fraction of free space that is reduced by the presence of obstacles. This increase in the local density of cells increases the propensity of neighboring cells to collide as they move, further truncating ℓ. Single-cell imaging in a porous medium confirms this expectation [28]: when the available free space is so tight that multiple cells cannot fit side-by-side, cells are necessarily restricted to end-on collisions between each other as they move, also inducing reorientations akin to those induced by collisions with surrounding obstacles. Therefore, as a first step toward incorporating this behavior into the model described in §2.1, we adopt a mean-field treatment of cell-cell interactions in which cells truncate each other’s directed steps in a density dependent manner, inducing transient trapping events again of duration τ_t akin to collisions with obstacles.
   In particular, wherever the local density b/ϕ is larger than a threshold value b*/ϕ such that the mean separation between the surfaces of neighboring cells, ℓ_cell, decreases below the mean chord length ℓ_c, we expect that cell-cell collisions truncate ℓ from ℓ_c to ℓ_cell (schematized in the middle inset of Fig 1). Because the diffusion and chemotactic coefficients both vary as ∝ ℓ², we therefore multiply both density-independent parameters D_b0 and χ₀ that characterize isolated cells by the density-dependent correction factor μ_crowd(b) = (ℓ_cell/ℓ_c)², where the cell separation is approximated as the mean value ℓ_cell ≡ (3ϕ/4πb)^1/3 − d; here, d ≈ 1 μm is the characteristic size of a cell, and therefore b* ≡ 3ϕ/[4π(ℓ_c + d)³]. As b increases further, it eventually reaches the jamming density b_jammed ≡ 3ϕ/(4πd³) at which cells cannot move at all, and ℓ_cell = 0; in this case, both transport parameters are zero, and the bacterial population can only spread via growth. Therefore, in Eq 2, D_b0 and χ₀ are replaced by the corrected values (5) (6) where the crowding correction factor μ_crowd(b) is piecewise defined as (7) as shown in Fig 1; the limits b* and b_jammed are indicated by the left and right vertical dashed lines, respectively. This way of correcting the diffusive term in Eq 2 represents a simplifying approximation; strictly speaking, one cannot simply commute the divergence operator with the diffusion coefficient, given that D_b depends on b(x, t) via the crowding correction factor. However, as we describe further in S1 Text, this simplification does not appreciably influence our results and conclusions. We term cases with low cell density (b < b*) the obstacle collisions limited regime described in (i) above; cases with intermediate cell density (b* ≤ b < b_jammed) the cell collisions limited regime; and cases with the highest possible density of cells (b = b_jammed) the jammed growth spreading regime described in (iii) below. Because we take the cells and surrounding obstacles to be incompressible, b cannot exceed b_jammed.
3. (iii). Jammed growth spreading. When they are jammed, cells form a contact network that holds them in place and prevents motion by active propulsion. However, these cells can continue to proliferate if supplied with nutrient; based on the experiments in [28], we assume that the maximal growth rate γ is not affected by confinement. Thus, in this case, their high body stiffness enables growing cells to push outward on their neighbors; the bacterial population can then be treated as an incompressible “fluid” in which the added stress due to cellular growth relaxes rapidly via spreading, as is conventionally done in models of growing immotile populations [41, 69–71] and supported by experiments [39, 40]. Because we treat the obstacles comprising the medium as being rigid and immovable, and the interstitial free space large enough for cells to move through without being deformed, this process leads to jammed growth spreading. We incorporate this behavior into the Keller-Segel model following previous work modeling the growth of immotile biofilms [42]. In particular, at each time step δt, we first identify the smallest x_i at which b(x_i, t + δt) exceeds b_jammed; we then set b(x_i, t + δt) = b_jammed and instead relocate the newly-formed cells δb(x_i) ≡ b(x_i, t + δt) − b_jammed to the nearest location x_j > x_i at which b(x_j, t) < b_jammed. We then repeat this process for all successive positions x > x_i such that at time t + δt, the upper limit on cell density b_jammed is globally satisfied.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Summary of the cell density-dependent crowding correction μ_crowd, which we use in the model to incorporate the influence of confinement on cell-cell collisions.

In particular, the cellular transport parameters are multiplied by μ_crowd, this case shown for the prototypical case of intermediate confinement. (Left) At low densities, spreading of cells (green) is impeded only by collisions with surrounding solid obstacles (grey), not with neighboring cells, so μ_crowd = 1. This impeded spreading is quantified by the transport parameters D_b0 and χ₀, whose values are regulated by the characteristic chord length ℓ_c characterizing the amount of free space between obstacles. (Middle) When the local density of cells is so large that the characteristic separation between neighboring cells ℓ_cell is less than the characteristic chord length ℓ_c, cell-cell collisions further truncate the transport parameters. This effect is quantified by μ_crowd < 1. (Right) At the maximal density b = b_jammed, the cells are jammed and have no free space to move. Therefore, μ_crowd = 0, and the population spreads solely through growth and division of cells. Note that our definition of the number density of bacteria b is as the number of cells per unit total volume of space, which includes the volume of surrounding obstacles.

https://doi.org/10.1371/journal.pcbi.1010063.g001

2.4 Implementation of numerical simulations

To explore the influence of confinement, we perform numerical simulations of Eqs 1 and 2, modified as described in §2.3. Motivated by its simplicity and amenability to the addition of our discrete jammed expansion rule, we implement a forward Euler method to solve these equations; specifically, we discretize the spatial coordinate x using a forward difference form for first derivatives and a central difference form for second derivatives. The update equations for nutrient concentration and bacterial cell density, corresponding to Eqs 1 and 2 respectively, are then: where time points advance in discrete steps of δt and are indexed by n, and spatial positions are separated by discrete steps of δx and are indexed by i. The spatial resolution δx is 10 μm and the time step δt is 0.01 s; as shown in S3 Fig, these choices are sufficiently fine so that our results are not sensitive to the choice of resolution. We note that implicit methods (such as backwards Euler) or semi-implicit methods would likely improve the efficiency of the numerical simulations—an important consideration for those seeking to extend our work e.g., to higher spatial dimensions.

To connect our results to an experimental system, we use input parameters and initial conditions that mimic the experiments described in [28], which explored the chemotactic spreading of E. coli populations in 3D porous media composed of densely-packed hydrogel particles. We use a Cartesian rectilinear coordinate system extending to a maximum distance of 1.75 × 10⁴ μm, matching the length of the experimental system. Because our system is one-dimensional, vectors (e.g. fluxes) oriented in the + or −x directions are represented by positive or negative quantities, respectively, with the vector notation suppressed. Both boundaries have no flux conditions. In these experiments, L-serine was considered to act as the primary nutrient and chemoattractant for cells. Because the hydrogel particles are polymer networks swollen in liquid, they are permeable to the nutrient, similar to many other naturally-occurring media such as biological gels and microporous clays/soils. Therefore, we take the nutrient diffusivity D_c to be equal to its value in bulk liquid, 800 μm² s⁻¹ [72], and the nutrient is initially saturated at c_∞ = 10 mM throughout the simulation domain. For all the simulations, we use direct measurements of individual cells [27, 28, 34] to choose fixed values of the cellular parameters c₋, c₊, and γ given by 1 μM, 30 μM, and 0.69 h⁻¹, respectively; furthermore, as detailed in S2 Text, we use the data from experiments on spreading populations [28] to directly determine c_char and κ, given by 10 μM and 1.3 × 10⁻¹² mM (cells/mL)⁻¹ s⁻¹, respectively.

Each experiment used a long 3D-printed cylinder of close-packed cells not containing hydrogel particles (ϕ = 1) as the initial inoculum, embedded within and surrounded by the 3D porous medium. The cells then continued to spread radially outward through the pore space. Thus, as the initial condition in all the simulations, we consider a Gaussian profile of b(x, t = 0) centered at x = 0 with a full width at half maximum of 100 μm and a peak number density of b₀ = b_max ≡ 3/(4πd³) = 2.4 × 10¹¹ cells/mL, where b_max is defined as the number density of close-packed cells and is therefore the maximum possible value of b₀—with the exception of the lower-density simulations that employ a lower value of b₀, as detailed further in §3.1.3. Hence, for all simulations except those in §3.1.3, the initial inoculum has a maximal cell density b_max > b_jammed, where b_jammed instead corresponds to the maximal possible density of cells in confinement (ϕ < 1). For simplicity, wherever b(x, t = 0) > b_jammed, we still apply the jammed growth spreading rule described in §2.3(iii), but with b_jammed replaced by b(x, t = 0).

The experiments tuned cellular confinement by using porous media with varying porosities ϕ and mean chord lengths ℓ_c [37], resulting in varying values of the transport parameters D_b0 and χ₀ [28, 38]. In particular, as determined from the experiments, D_b0 and χ₀ both decrease with increasing confinement as cellular mobility is increasingly hindered. Hence, in our simulations, we tune confinement by varying these parameters, using the values of D_b0 obtained from single-cell imaging [28] and extracting χ₀ from experimental measurements of population spreading, as detailed in S2 Text. The confinement-dependent parameters are summarized in Table 1.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Parameters used to describe bacteria in weak, intermediate, and strong confinement, as defined in the text.

All parameters are defined in the text and their values are obtained from experiments as detailed in §2.4 and S2 Text, with the exception of , which is determined directly from the simulation. We note that the values of α₀ are taken directly from the experiments in [28], which indicate that this parameter surprisingly decreases with increasing confinement. Thus, while we expect that both transport parameters χ₀ and D_b0 are tuned by confinement in a similar way, with both proportional to , it appears from the experiments that the ratio of the proportionality constants for each is also confinement-dependent. That is, experiments suggest that confinement more strongly hinders directed (quantified by χ₀) versus undirected (quantified by D_b0) spreading. While more work needs to be done to fully unravel why this is the case, in absence of a theoretical model for the confinement-dependence of α₀, we directly use the experimental values in our work. Furthermore, in the absence of any experimental data assessing the influence of cell density on α, we make the simplest possible assumption that cell-cell collisions hinder both χ₀ and D_b0 through the same crowding correction factor μ_crowd(b), which quantifies the reduction in free space available to the cells. Thus, we take α = α₀. Future experiments could further probe this density dependence and motivate the introduction of additional extensions to our model.

https://doi.org/10.1371/journal.pcbi.1010063.t001

The corresponding dimensionless parameters characterizing the Keller-Segel model (§2.2) are also summarized in Table 1:

• The diffusion parameter δ₀ ≡ D_c/D_b0 increases with confinement as cellular mobility is increasingly hindered. For all conditions tested here, however, δ₀ is always much greater than one, reflecting fast diffusion of nutrient; thus, we expect that nutrient levels vary over large spatial extents, as confirmed in the simulations that follow.
• For all conditions tested here, the directedness parameter α₀ ≡ χ₀/D_b0 is always much greater than one, indicating that motile cells strongly direct their motion in response to the nutrient gradient established through consumption. Intriguingly, the α₀ determined from the experimental parameters decreases with increasing confinement, indicating that confinement more strongly hinders directed versus undirected motion—consistent with previous reports that confinement fundamentally alters the mechanism by which cells perform chemotaxis [28]. Further investigating the determinants of α₀ in confinement will be a useful direction for future experiments.
• Because the maximal growth rate is not affected by confinement [28], the yield parameter β₀ ≡ γ/(b₀κ/c_∞) is independent of confinement for all of our simulations. For all simulations employing b₀ = b_max, β₀ is much less than one, reflecting the fact that nutrient consumption by a maximally dense population is faster than cellular proliferation; conversely, for the lower-density simulations presented in Figs 4 and 5, β₀ is much greater than one, indicating the dominant role of proliferation in this case.

Therefore, for our simulations testing the influence of confinement on bacterial spreading, the parameter Λ₀ ≡ α₀/(β₀δ₀) varies over a broad range, decreasing over nearly three orders of magnitude as confinement increases. We note that because the different parameters δ₀, α₀, β₀ do not incorporate the influence of density-dependent cellular crowding, we do not expect this transition to occur precisely at Λ₀ ≈ 1. We therefore define a new version of this parameter, Λ ≡ α/(βδ), where now , , and (Table 1); is defined as the long-time mean cell density within each propagating pulse, and is directly calculated from each simulation as described further below. Thus, the newly-defined Λ explicitly incorporates density-dependent crowding. As summarized in Table 1, our simulations explore the transition from weak confinement (Λ = 0.95) to strong confinement (Λ = 3.6 × 10⁻⁴); consistent with our expectation, this range reflects a transition from chemotactic to growth-driven spreading, as demonstrated directly by the simulations presented below.

3.1 Intermediate confinement

As a prototypical starting case, we first examine bacterial spreading from a dense-packed Gaussian-shaped inoculum under intermediate confinement (Λ = 0.037), shown by the initial profile for t = 0 in Fig 2A. The simulation incorporates both motility and growth. The cells rapidly deplete nutrient locally via consumption (S1 Video) over a time scale ∼ c_∞/(κb_max) ≈ 30 s, establishing a steep nutrient gradient at the leading edge of the population. This gradient extends over a large distance ahead of the population (Fig 2B and inset)—as expected from our calculation of the diffusion parameter δ₀ ≫ 1. Cells at this leading edge then continue to grow outward as a jammed front with b = b_jammed, shown by the flat region at t = 1.8 h in Fig 2A. Eventually, a lower-density, coherent pulse of cells detaches from this jammed region (t = 3.7 h), continues to propagate the nutrient gradient along with it, and thus continues to spread outward (t > 3.7 h), as shown by the outward-moving peak in Fig 2A.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Results from a numerical simulation of population spreading in intermediate confinement.

The simulation incorporates both motility and growth. (A) shows the dynamics of the cells while (B) shows the corresponding dynamics of the nutrient, quantified by the normalized density b/b_max and concentration c/c_∞, respectively. As noted in §2.4, the initial inoculum is composed purely of dense-packed cells with liquid between them (ϕ = 1), with the entire inoculum surrounded by the obstacle-filled medium (ϕ = 0.17); hence, the initial inoculum has b = b_max, which is larger than b_jammed, the jamming density of cells in confinement. Different colors indicate different times as listed. The dense inoculum initially centered about the origin spreads outward, first as a jammed front (jamming density shown by the dashed grey line in A), then detaching as a coherent lower-density pulse that propagates continually via chemotaxis. At long times, this pulse appears to approach an unchanging shape and speed, as suggested by the collapse of the profiles in the upper inset (showing the same data, but shifted horizontally to center the peaks). The cellular dynamics arise in response to consumption of the nutrient, which is initially saturated everywhere, but is rapidly depleted and forms a gradient that is propagated with the pulse (inset shows the same data but with both axes zoomed out). In B, the three dashed grey lines show the characteristic concentrations of sensing c₊ and c₋ and the characteristic Monod concentration c_char; the corresponding positions are shown by the pluses, diamonds, and triangles, respectively, in A-B. An animated form of this Figure is shown in S1 Video, and a sample MATLAB code to implement this simulation is provided in S1 File.

https://doi.org/10.1371/journal.pcbi.1010063.g002

Indeed, this pulse spans the extent over which nutrient varies between the upper and lower bounds of sensing, c₊ and c₋ (pluses and diamonds shown for the t = 11 h profiles, respectively)—reflecting the central role of chemotaxis in driving its propagation. The forward face of the pulse is also exposed to sufficient nutrient for cells to proliferate (with c ≥ c_char, the characteristic Monod concentration, shown by the triangles on the t = 11 h profiles)—suggesting that cellular growth contributes to population spreading over long time scales, as well. The overall width of this pulse, W ≈ 200 μm, is set by the length scale over which nutrient is depleted by consumption; at its rear, the nutrient concentration and nutrient gradient are both low, causing both growth and chemotaxis to be hindered. As a result, cells are shed at a near-constant density b_trailing ≈ 0.02b/b_max (see 0.5 mm < x < 1.2 mm in Fig 2A). This coherent pulse of cells continues to move apparently without an appreciable change in shape, as suggested by the inset to Fig 2A, at a speed v_pulse ≈ 0.15 mm/h. However, given the limited duration of the simulation, our results do not enable us to definitively conclude that the simulated pulse develops into a traveling wave with an unchanging shape; building on our simulation to explore longer length and time scales to test this possibility will be a useful direction for future work. The nutrient profile concomitantly propagates with the pulse, as shown in Fig 2B. Notably, similar spreading behavior was observed in experiments [28]; as shown in S4 Fig, in both simulation and experiment, we observe similarly-shaped bacterial pulses with comparable widths, trailing densities (compared to the peak densities), and final positions at t ≈ 11 h.

3.1.1 Initial dynamics.

To further characterize these spreading dynamics, we track the position x_l of the leading edge of the population over time t, as shown in Fig 3A. Specifically, motivated by a similar definition used in prior experiments [28], we define x_l as the position at which b falls below the threshold value b/b_max = 10⁻⁴. Initially, population spreading is hindered (x_l ∼ t^ν with ν ≪ 1 e.g., red point), but as the coherent pulse forms and propagates, it eventually approaches constant speed spreading (ν ≈ 1 e.g., blue point). A similar transition from hindered to constant speed spreading was observed in experiments [28], although the underlying reason has thus far remained unclear. Here, we use our model to clarify the origin of this transition.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Population dynamics, morphology, and fluxes driving spreading for the simulation of bacteria in intermediate confinement (Fig 2).

The simulation incorporates both motility and growth. (A) Increase in the position of the leading edge of the population is initially hindered (red), but approaches constant-speed motion indicated by the triangle at long times (blue). The corresponding instantaneous speed v_pulse is shown in the inset. (B) At a short time corresponding to the red point in A, the population expands as a jammed front (top panel). Lower panel shows that cellular growth and diffusion are the primary contributors to the expansion of this front. (C) At a long time corresponding to the blue point in A, the population spreads as a coherent pulse (top panel). Lower panel shows that chemotaxis is the primary contributor to pulse propagation. Positions of the maximal diffusive and chemotactic fluxes are indicated by the circles and squares, respectively, in B-C; note the slight upward kink in the diffusive flux in C indicated by the circle. (D) Variation of the maximal diffusive and chemotactic fluxes, indicated by the circles and squares in B-C, over time. The initial population dynamics are dominated by cellular diffusion (circles), while at longer times chemotaxis dominates (squares). To illustrate the role played by cellular collisions, we show the same data with and without the crowding correction μ_crowd in the upper (grey) and lower (black) datasets; crowding hinders population spreading, as shown by the vertical offset in the curves, but plays a less appreciable role at long times, as shown by the curves approaching each other.

https://doi.org/10.1371/journal.pcbi.1010063.g003

In particular, we examine the two different contributions to the motility-driven flux of cells—active diffusion and chemotaxis—for the population at early and late times (Fig 3B and 3C, respectively); for simplicity, we do not consider the added influence of growth, which only plays an appreciable role for long times t ≫ γ⁻¹, until the next subsection. The magnitude of the active diffusive flux −D_b(b)∇b = −D_b,0μ_crowd(b)∇b as it varies across the population is shown by the dashed lines in the bottom panels of Fig 3B and 3C, while the magnitude of the chemotactic flux bv_c = bχ₀μ_crowd(b)∇f(c) is shown by the dash-dotted lines instead. At early times, the gradient in cell density is steep, as set by the sharp initial profile of cells and the limited extent of subsequent population spreading (Fig 3B, top). As a result, spreading is primarily due to active diffusion, which dominates over chemotaxis, as shown in the lower panel of Fig 3B. By contrast, as cells spread outward, the gradient in cell density becomes less steep. As a result, at late times, spreading is primarily due to chemotaxis, which dominates over active diffusion, as shown in the lower panel of Fig 3C. This behavior is also reflected by the bottom set of circles and squares in Fig 3D, which represent the maximal diffusive and chemotactic fluxes across the population (exemplified by the circles and squares in Fig 3B and 3C) over time. Initially, the diffusive flux dominates over the chemotactic flux; however, as the population continues to spread and consume nutrient, the diffusive flux decreases and the chemotactic flux increases, with both eventually approaching constant values at long times.

Another key factor that hinders the initial population spreading is cellular crowding. To assess the influence of crowding, we compare the maximal diffusive and chemotactic fluxes across the population, but with or without the crowding correction factor μ_crowd(b) (corresponding to the “with crowding” and “no crowding” datasets, respectively, in Fig 3D). In both cases, the active diffusive flux dominates over the chemotactic flux initially, but chemotaxis eventually dominates as the population continues to spread and establish the nutrient gradient (e.g. top set of squares for t ≤ 0.1 h). The spreading of the population remains hindered, however; due to the high initial density of cells, crowding continues to limit the chemotactic flux of cells, only enabling a small fraction at the leading edge of the population to spread outward—as exemplified by the first two profiles in Fig 2A, the sharp decrease in both diffusive and chemotactic fluxes for x < 0.5 mm in Fig 3B, and the large difference between the two sets of squares in Fig 3D. Eventually, as this leading edge continues to spread, crowding in the forward face of the population becomes sufficiently low, enabling the coherent pulse of cells to detach from the population—as exemplified by the t = 3.7 h profile in Fig 2A and the “kink” in the top set of squares at t ≈ 3 h in Fig 3D. Hindrance due to crowding continues to decrease over time, as shown by the diminishing difference between the two sets of squares in Fig 3D for t > 3 h, and eventually approaches a constant value.

Hence, population spreading is initially slow due to the time required for cellular consumption to establish a sufficiently strong nutrient gradient to drive chemotactic spreading. Cellular crowding near the initial inoculum then continues to hinder spreading until enough of the forward face of the population has spread outward—enabling cells to detach as a coherent pulse that continues to move outward, eventually approaching a constant speed.

3.1.3 Influence of initial cell density.

Our analysis thus far considered a dense initial inoculum, for which cellular crowding hinders the formation and detachment of a pulse; at much longer times, this less-crowded pulse no longer resembles the initial inoculum, but instead is shaped by the interplay of chemotaxis and growth (Fig 3). We therefore expect that for a lower-density inoculum, a similar pulse also emerges at long times, but with initial dynamics that are limited instead by the time required for cellular consumption to establish a sufficiently strong nutrient gradient. To test this expectation, we repeat the simulation shown in Fig 2, but using an initial peak number density of cells that is 10⁵ times smaller (b₀ = 10⁻⁵b_max)—shown in S2 Video.

In the previously-considered case of a dense inoculum, the cells deplete nutrient rapidly via consumption, and the population subsequently spreads from its leading edge as a growing jammed front (first two curves in Fig 2). By contrast, with a more dilute inoculum, nutrient depletion takes much longer. Instead, the population continually grows and spreads as a whole (red to blue curves in Fig 4A inset), not just at its leading edge, without appreciably depleting nutrient. It eventually reaches a maximal density b′ = b₀e^γt′ for which the time scale of subsequent nutrient depletion t_dep ∼ c_∞/(κb′) is comparable to the time scale of subsequent growth t_g ∼ γ⁻¹; equating these time scales yields . Therefore, we expect that nutrient is fully depleted at the initial inoculum after . The simulation results are consistent with this estimate, which neglects spatial variation in nutrient availability through the entire population and thus serves as a lower bound, showing that nutrient is fully depleted at the initial inoculum after ∼ 14 h (red to blue curves in Fig 4B inset). The nutrient gradient again extends over a large distance ahead of the population, as expected from our calculation of the diffusion parameter δ₀ ≫ 1.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Results from a numerical simulation of population spreading in intermediate confinement starting from a more dilute inoculum.

The simulation incorporates both motility and growth. (A) shows the dynamics of the cells while (B) shows the corresponding dynamics of the nutrient, quantified by the normalized density b/b_max and concentration c/c_∞, respectively. Different colors indicate different times as listed. The dilute inoculum (jamming density shown by the dashed grey line in A) initially centered about the origin first grows exponentially and spreads diffusively until nutrient is locally depleted (upper inset shows the same data, but zoomed in to the vertical axis); only then does a coherent pulse detach and propagate continually via chemotaxis in response to the nutrient gradient (lower inset shows the same data but with both axes zoomed out). Even though the short-time behavior is different from the case of a more dense inoculum shown in Fig 2, the long-time behavior of this pulse is identical. To facilitate comparison with Fig 2, in B, the three dashed grey lines again show the characteristic concentrations of sensing c₊ and c₋ and the characteristic Monod concentration c_char; the corresponding positions are shown by the pluses, diamonds, and triangles, respectively, in A-B. An animated form of this Figure is shown in S2 Video.

https://doi.org/10.1371/journal.pcbi.1010063.g004

Unlike the case of a dense inoculum, the population does not subsequently spread as a jammed front. Instead, once the nutrient gradient is sufficiently strong, a coherent pulse of cells again detaches without the prior formation of a jammed front, continues to propagate the nutrient gradient with it, and continues to spread outward (t > 15 h in Fig 4). Consistent with our expectation, this pulse is noticeably similar to that which arises in the dense inoculum case: it has a nearly-identical shape and also appears to move without an appreciable change in shape, eventually reaching approximately a similar constant speed v_pulse ≈ 0.1 mm/h (compare late-time profiles in Figs 2 and 4). Evaluating the speed by instead tracking the position of the peak, instead of the leading edge, also yields a comparable value of v_pulse ≈ 0.16 mm/h.

To further characterize the population spreading dynamics, we again plot the leading edge position x_l as a function of time t. As in the case of a dense inoculum, x_l ∼ t^ν with ν ≪ 1 at early times, transitioning to ν ≈ 1 at later times (Fig 5A); however, these seemingly similar dynamics reflect fundamentally different underlying processes at early times. With a more dilute inoculum, slower nutrient depletion causes the diffusive flux to initially dominate over chemotaxis (Fig 5B) without any influence of cellular crowding—indicated by the overlap of the early-time points with/without the crowding correction in Fig 5D. As cells continue to grow and consume nutrient, they eventually establish a sufficiently strong gradient and spread as a coherent pulse via chemotaxis—as indicated by the dominant role of the chemotactic flux at long times (Fig 5C and 5D). At these later times, the different contributions to the bacterial flux are nearly identical to those that drive pulse propagation in the case of a dense inoculum (compare Figs 5C and 3C). Indeed, the fractions of the overall pulse speed attributable to chemotaxis and growth, as quantified by Eq 9, are ≈ 58% and 42%, respectively—nearly identical to the case of a dense inoculum. Hence, while the initial dynamics of population spreading are sensitive to the initial cell density—consistent with experiments [28]—the properties of the pulse that forms and continues to drive spreading at long times are not, instead being set solely by the interplay between chemotaxis and growth.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Population dynamics, morphology, and fluxes driving spreading for the simulation of bacteria in intermediate confinement, but from a more dilute initial inoculum (Fig 4).

The simulation incorporates both motility and growth. (A) Increase in the position of the leading edge of the population is initially hindered (red), but approaches constant-speed motion indicated by the triangle at long times (purple). The corresponding speed v_pulse is shown in the inset. (B) At a short time corresponding to the red point in A, the population grows exponentially (top panel) and spreads primarily through growth and diffusion (lower panel). (C) At a long time corresponding to the purple point in A, the population spreads as a coherent pulse (top panel). Lower panel shows that chemotaxis is the primary contributor to pulse propagation. Positions of the maximal diffusive and chemotactic fluxes are indicated by the circles and squares, respectively, in B-C; note the slight upward kink in the diffusive flux in C indicated by the circle. (D) Variation of the maximal diffusive and chemotactic fluxes, indicated by the circles and squares in B-C, over time. The initial population dynamics are dominated by cellular diffusion (circles), while at longer times chemotaxis dominates (squares). To illustrate the role played by cellular collisions, we show the same data with (black) and without (grey) the crowding correction μ_crowd in the datasets indicated by the red lines; the data are identical except at long times, when crowding slightly hinders population spreading.

https://doi.org/10.1371/journal.pcbi.1010063.g005

3.2 Influence of confinement

For the case of intermediate confinement explored thus far, we have established that chemotaxis and growth both drive population spreading at long times. How does this behavior change with confinement? As quantified in Figs 3D and 5, confinement-induced crowding limits the chemotactic flux; therefore, we expect that with reduced or increased confinement, chemotaxis or growth plays a more dominant role in driving spreading, respectively. To test this expectation, we perform the same simulation with a dense inoculum as in Figs 2 and 3, but with different values of the confinement-dependent parameters as summarized in Table 1. In particular, our simulations explore Λ = 0.95, 0.037, and 3.6 × 10⁻⁴, representing weak, intermediate, and strong confinement (top, middle, and bottom rows in Figs 6 and 7), respectively—also shown in S3 Video.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Increasing confinement causes a transition from fast chemotactic pulse propagation to slower jammed growth expansion.

Panels show results from numerical simulations of population spreading from the same dense inoculum initially centered about the origin in weak, intermediate, and strong confinement, shown by top, middle, and bottom rows respectively. First column shows the results of the full model, while second and third columns show the same simulations with growth or chemotaxis omitted, respectively. We only show the normalized cellular density b/b_max for clarity. As noted in §2.4, the initial inoculum is composed purely of dense-packed cells with liquid in between them (ϕ = 1), surrounded by the obstacle-filled medium (ϕ < 1); hence, the initial inoculum has b = b_max, which is larger than b_jammed, the jamming density of cells in confinement. Different colors indicate different times as listed in the color scale. In weak confinement, a coherent pulse rapidly detaches and continually propagates; this pulse is driven primarily by chemotaxis, and thus, omitting growth barely changes the dynamics while omitting chemotaxis abolishes the propagation altogether. Conversely, in strong confinement, the population spreads slowly as a jammed front, driven primarily by growth. In intermediate confinement, both growth and chemotaxis drive population spreading. The dashed grey line shows the jamming density, which varies depending on confinement (and is larger than the vertical scale in the top row). An animated form of this Figure, along with the nutrient profiles, is shown in S3 Video.

https://doi.org/10.1371/journal.pcbi.1010063.g006

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Population dynamics, morphology, and fluxes driving spreading for simulations of bacteria in weak, intermediate, and strong confinement (Fig 6) as shown by the top, middle, and bottom rows, respectively.

(A-B) Increase in the position of the leading edge of the population is initially hindered, but approaches constant-speed motion indicated by the triangle at long times. In strong confinement (C), however, the long-time behavior approaches diffusive-like scaling instead. (D-E) At long times, the population spreads as a coherent pulse (solid line) driven primarily by chemotaxis in weak confinement, and by both chemotaxis and growth in intermediate confinement. (F) In strong confinement, however, the population spreads slowly as a jammed front, driven primarily by growth.

https://doi.org/10.1371/journal.pcbi.1010063.g007

In all cases, the cells first rapidly deplete nutrient locally via consumption, generating a nutrient gradient that again extends over a large distance and drives subsequent spreading at the leading edge of the population. However, consistent with our expectation, and with experimental observations [28], the nature of this spreading is strongly confinement-dependent.