Chaotic dynamics and oxygen transport in thin films of aerotactic bacteria

A kinetic model and three-dimensional numerical simulations are applied to study the dynamics in suspensions of run-and-tumble aerotactic bacteria confined in freestanding liquid films surrounded by air. In thin films, oxygen and bacterial concentration profiles approach steady states. In thicker films, a transition to chaotic dynamics is shown to occur and is characterized by unsteady correlated motions, the formation of bacterial plumes, and enhanced oxygen transport and consumption. This transition, also observed in previous experiments, arises as a result of the coupling between the aerotactic response of the bacteria and the flow fields they generate via hydrodynamic interactions. C © 2012 American Institute of Physics. [http://dx.doi.org/10.1063/1.4752764]

The ability of motile bacteria to modulate their motions in response to changes in concentration of chemical cues or nutrients plays a central role in the development and growth of micro-organismal colonies, as well as in cell-cell communication processes such as quorum sensing. 1 The directed migration of motile bacteria such as Escherichia coli and Bacillus subtilis along a chemical gradient, or "chemotaxis," is a consequence of their run-and-tumble dynamics, by which periods of motion along straight lines ("runs") are interspersed by random reorientation events ("tumbles") resulting from the unbundling and rebundling of their flagella. 2,3 y modulating the frequency of tumbling events based on experienced temporal changes in chemical concentration, bacteria are able to undergo a biased random walk in the gradient direction.In large-scale suspensions, this mode of locomotion often results in complex spatial patterns, with accumulation near sources of oxygen and nutrients, 4 or the formation of stable multicellular aggregates. 5s they locomote through fluid environments, bacteria not only interact with chemical substances but can also interact with each other in sufficiently concentrated suspensions via fluidmediated hydrodynamic interactions.A neutrally buoyant self-propelled particle swimming through a viscous medium exerts a net force dipole on the fluid, 6 thereby driving a flow that can affect the motions of its neighbors.In a suspension of many swimmers, the resulting fluid flows can be quite strong and overcome the individual swimming motions.0][11] While these effects have been studied both theoretically [12][13][14][15][16][17] and numerically, [18][19][20] their impact on nutrient transport and bacterial chemotactic response has not yet been addressed.
In recent experiments, Sokolov et al. 21studied suspensions of the bacterium Bacillus subtilis confined in free-standing liquid films.As a result of aerotaxis (or directed response to an oxygen stimulus), the bacteria were found to accumulate near the interfaces where the oxygen concentration was the highest.In thin films ( 100 μm), a steady concentration profile was reached with a spatially homogeneous bacterial distribution in the plane of the film.As the thickness was increased (100-500 μm), a transition to inhomogeneous three-dimensional chaotic motions took place, with unsteady plumes of bacteria penetrating the bulk of the film and enhanced oxygen transport.In yet thicker films ( 500 μm), a stagnant bacterial layer at the center of the film was also reported.As argued by Sokolov et al., these dynamics are likely a consequence of the coupling between the aerotactic response of the bacteria and the flow fields set-up via hydrodynamic interactions.Recent attempts at elucidating this coupling have focused on linear stability analyses in uniform gradients, 22,23 and have shown that the convective motion driven by bacterial active stresses may be responsible for these observations.However, these models did not account for oxygen transport and were not able to capture the strongly nonlinear regimes of the experiments.
In this Letter, we address this problem by extending a previous kinetic model for active suspensions [12][13][14] to account for run-and-tumble aerotaxis in thin liquid films, and explicitly model the three-way coupling between bacterial dynamics, oxygen transport, and fluid flow.Using threedimensional numerical simulations, we study the emergence of nonlinear dynamics and pattern formation in these films, and are able to quantitatively capture and explain the three regimes reported in experiments 21 and to correlate them with the prediction from the linear stability analysis of Kasyap and Koch. 23he configuration of a suspension of motile bacteria is modeled by a distribution function (x, p, t) of center-of-mass position x and orientation p, normalized as 1 V V dx dp (x, p, t) = n, where n is the mean number density in the domain of volume V .The distribution function satisfies a conservation equation ][26] Specifically, the first term −λ describes bacteria that tumble away from orientation p (with tumbling rate λ), whereas the integral term accounts for tumbling from orientation p to p.The flux velocities ẋ and ṗ in Eq. ( 1) are modeled as 13 The center-of-mass velocity in Eq. ( 2) results from the constant self-propulsion velocity U 0 p (in the direction of p), the local fluid velocity u(x, t), and translational diffusion with diffusivity D. Similarly, Eq. ( 3) accounts for the rotation of the particle in the local flow field via Jeffery's equation, 27 where E = (∇u + ∇u T )/2 and W = (∇u − ∇u T )/2 are the rate-of-strain and vorticity tensors, respectively, and γ ∈ [ − 1, 1] is a shape parameter.For a spheroidal particle, γ = (r 2 − 1)/(r 2 + 1), where r is the aspect ratio, and γ ≈ 1 for a slender particle such as most types of bacteria.We also account for rotary diffusion with diffusivity d.
The fluid velocity u appearing in Eqs. ( 2) and ( 3) is driven by the active stresses generated by the swimming bacteria, and satisfies the steady Stokes equations where q is the pressure and η is the viscosity of the medium.The active particle stress tensor p (x, t) captures the effect of the force dipoles due to self-propulsion and can be modeled as 28, 29 where the stress magnitude σ 0 is a negative constant for rear-actuated swimmers, also known as pushers, such as most motile bacteria.We model oxygen transport using an advection-diffusion-reaction equation for the oxygen concentration field s(x, t) (normalized by its constant value at the film boundaries, taken to be the oxygen saturation concentration) In Eq. ( 6), oxygen is advected by the local flow and diffuses with diffusivity D 0 .Consumption by the bacteria is modeled as a second-order reaction, where κ is the consumption rate in oxygen-rich environments, c(x, t) = dp (x, p, t) is the bacterial concentration, and the function f(s) = max {0, [(s − 0.1)/0.9] 0.2 } accounts for the experimental observation that oxygen consumption diminishes drastically when its concentration falls below ≈10% of its saturation value. 30inally, the tumbling bias, which couples bacterial dynamics to oxygen concentration and leads to aerotaxis, is modeled in Eq. ( 1) by letting the tumbling rate depend on s as 3,24 λ = λ 0 exp(−ξ D t s), where ξ is the aerotactic response strength, and D t s is a Lagrangian derivative capturing the rate of change of s as experienced by a bacterium moving with velocity U 0 p + u: Bacteria for which D t s > 0 therefore have a lower tumbling rate, causing them to migrate towards oxygen-rich regions.Within this framework, the tumbling rate only depends on the instantaneous rate of change of s along bacterial trajectories rather than on the time history of the concentration field experienced by the swimmers.Equations ( 1)-( 6) form a closed system for the coupled dynamics of bacteria, oxygen, and fluid flow, and are the basis of this study.A similar system was also recently proposed by Lushi et al. 31 to model suspensions of chemotactic bacteria where the chemoattractant is secreted by the swimmers.
We solve this system numerically in a doubly periodic geometry (periodic in the x-y plane, bounded in the z-direction) mimicking the liquid films of the experiments of Sokolov et al. 21We denote by L the film thickness, with z = 0, L corresponding to the free surfaces.The boundary condition on is expressed as ẑ • ẋ dp = 0 at z = 0, L, and ensures that the mean concentration of bacteria is conserved.For the fluid flow, a free-shear-stress boundary condition is used to model air-liquid interfaces: u z = 0 and ∂ z u x = ∂ z u y = 0 at z = 0, L. Both boundary conditions are easily implemented using a reflection condition for the distribution function.The normalized oxygen concentration has a prescribed value of 1 at both interfaces.Spectral solutions of Eqs. ( 4) and ( 6) are used, coupled to a finite-difference solution (second order in time, space, and orientation) of the conservation equation for the distribution function. 32The initial condition for (x, p, 0) is uniform isotropic ( 0 = n/4π ) with a weak random perturbation, while the oxygen concentration inside the film is initially zero.
Parameter values are chosen for direct comparison with the experiments of Sokolov et al. 21The bacteria are Bacillus subtilis, with length = 4 μm, U 0 = 20 μm/s, σ 0 = −2.7 × 10 −18 kg m 2 /s, 33 and a mean number density of n = 2 × 10 16 cell/m 3 . 21Based on the experiments, we estimate the translational diffusivity to be D = 1.3 × 10 −10 m 2 /s, and an estimate for the rotational diffusivity can then be obtained from the generalized Taylor dispersion 19 relation d = U 2 0 /6D ≈ 0.5 s −1 .The oxygen diffusivity and consumption rate are set to D 0 = 2 × 10 −9 m 2 /s and κ = 6.7 × 10 −18 m 3 /cell/s, respectively. 21There is some uncertainty as to the actual values of the base tumbling rate and chemotactic strength in bacterial suspensions; we choose the values λ 0 = 1.9 s −1 and ξ = 3.1 s, which are close to previous estimates from experiments. 2 Finally, the film thickness is varied over the range L = 50-600 μm, which is also the range considered in the experiments of Sokolov et al. 21s the simulations start, oxygen rapidly diffuses into the film on a time scale L 2 /8D 0 of the order of a few seconds, and consumption by the bacteria begins.Both effects drive the formation of oxygen gradients inside the film, to which the bacteria respond by modifying their run-and-tumble dynamics so as to migrate preferentially towards the oxygen-rich regions near the boundaries.This aerotactic response leads to anisotropic orientation distributions, thereby inducing active stresses that, in some cases, can drive three-dimensional flows.
Simulations for various film thicknesses are illustrated in Fig. 1 (and accompanying video), showing typical bacterial and oxygen concentration fields at an arbitrary time past the initial transient.Based on the film thickness L, we distinguish three regimes with qualitatively different dynamics.In regime I, which is shown in Figs.1(a) and 1(b) and corresponds to thin films (L 200 μm), bacterial and oxygen concentrations quickly reach steady-state profiles that are uniform in the x-y plane.At steady state, diffusion across the film leads to a nearly uniform oxygen profile, and by consequence to a weak bacterial migration.This weak aerotaxis, along with the strong confinement which constrains the growth of instabilities and limits flows in the z-direction, explains the lack of unsteady dynamics in this case.As L is increased to ≈200-400 μm, regime II Phys. is observed and is characterized by the emergence of unsteady dynamics on the scale of the system [Figs.1(c) and 1(d)].After an initial transient, a statistical steady state is reached when the effects of active input power, oxygen diffusion, and consumption by the bacteria are balanced.The bacteria are found to concentrate significantly near the boundaries, but their distribution constantly fluctuates in time with the formation of dense unsteady plumes penetrating the bulk of the film and thereby resulting in fluid mixing.The oxygen concentration away from the boundaries is diminished, which further enhances aerotaxis.Above L 400 μm [Figs.1(e) and 1(f)], regime III is also characterized by a strong bacterial migration towards the boundaries and unsteady flow dynamics, but also exhibits a strong oxygen depletion layer and an additional accumulation of bacteria near the center of the film (see Fig. 2).The bacterial plumes are found to intensify, and fluid velocities are reduced near the centerline.The existence and characteristics of all three regimes are consistent with the experimental observations of Sokolov et al., 21 as are the critical film thicknesses above which they arise.
Figure 2 shows vertical profiles of various fields averaged over time (after the initial transient) and in the plane of the film.We compare these profiles to a "base-state" equilibrium solution of the kinetic equations corresponding to a steady distribution function that only depends on z and θ = cos −1 p z , with no fluid flow.Bacterial concentration profiles are plotted in Fig. 2(a), and exhibit strong peaks near the boundaries owing to aerotaxis.The profiles become sharper as L increases, yet they are not as sharp as in the base state, a consequence of the large-scale flows that emerge in thick films and facilitate transport of bacteria into the bulk.Interestingly, we find that the oxygen profiles in Fig. 2(b) are only weakly affected by interactions and remain close to the base state for all values of L; an estimate of the Péclet number for oxygen transport, Pe o = LU 0 /D 0 , yields values of order O(1), suggesting that this is likely a result of the strong oxygen diffusion.While s(z) is nearly uniform in thin films (hence the weak aerotaxis in regime I), strong gradients form in thicker films.In regime III (L = 600 μm), an oxygen depletion layer (defined as s(z) ≤ 0.1) forms in the center of the film, where consumption nearly ceases.The appearance of this layer is accompanied by a slight accumulation of bacteria at the centerline [see Fig. 2(a)].This accumulation should really be interpreted as a lack of depletion, as bacteria in this region no longer experience oxygen gradients and therefore become incapable of migrating away by aerotaxis.A similar dense immobile layer was reported in 530 μm films by Sokolov et al. 21he oxygen consumption K(z) = cf[s] (where • denotes averaging over x, y, and t) is plotted in Fig. 2(c).Consumption occurs primarily near the boundaries as expected.For all film thicknesses, hydrodynamic interactions are found to increase consumption slightly near the center of the film with respect to the base state, as a result of the enhanced transport by the flow.In very thick films (regime III), however, consumption ceases completely near the centerline (where s ≤ 0.1), leading to bacterial accumulation.
The emergence of large-scale flows with increasing L is illustrated in Fig. 2(d), showing the rms u z (z) of the vertical velocity across the film.There is no flow in regime I, which coincides with the base state.In regime II, significant velocities [∼O(U 0 )] are seen across most of the film, with the effect of enhancing transport from the boundaries into the bulk (bacterial plumes).However, the profile shape changes in regime III as the depletion layer forms, with high velocities now localized near the boundaries and low velocities in the center.In both regimes II and III, horizontal velocities (not shown) are also significant and can reach values of up to ≈5U 0 .
The unsteady dynamics arising in thick films are directly linked to the coupling between the bacterial aerotactic response and hydrodynamic interactions, and cannot be explained solely based on either effect: in fact, setting either ξ = 0 (unbiased tumbling) or u = 0 (no hydrodynamic interactions) in the simulations has the effect of stabilizing the suspensions for all values of L considered here.A mechanism for such taxis-driven instabilities was recently proposed by Kasyap and Koch 23 for bacteria swimming in a uniform gradient, where they argued that the anisotropic orientation distributions resulting from chemotaxis yield active stresses that drive flows tending to reinforce density fluctuations in directions normal to the chemical gradient.Based on a linear analysis, they showed that the onset of instability is governed by two dimensionless groups:  3(c) and compared to β crit .The marginal stability curve is found to predict the onset of unsteady dynamics in the simulations remarkably well (transition from regime I to II).We find that this transition is primarily associated with an increase in Pe b (or aerotactic drift), whereas the transition from II to III is accompanied by a sharp increase in β (or active stresses).
Using a kinetic model and numerical simulations, we have studied the transition to chaotic dynamics occurring in liquid films of aerotactic bacteria as film thickness is increased.Our results, which are consistent with experiments 21 and a linear stability theory, 23 emphasize the importance of active stresses in these suspensions, which can drive unsteady flows and enhance both bacterial and oxygen transport across the film.A similar conclusion was also reached by Lushi et al., 31 who studied the effects of self-generated flows on chemotactic aggregation when the chemoattractant is secreted by the swimmers: they find that, in suspensions of pushers such as bacteria, self-generated flows have a strong impact on the aggregation patterns, which become unsteady as a result of strongly mixing flows advecting both swimmers and chemoattractant.These observations all suggest that mixing and transport by hydrodynamic interactions likely play a significant role in facilitating the access of motile bacteria to oxygen and other nutrients in large-scale suspensions.

FIG. 2 .
FIG. 2. Vertical profiles for different film thicknesses L: (a) bacterial concentration c(z)/n, (b) oxygen concentration s(z), (c) oxygen consumption K(z)/n, and (d) rms of the vertical fluid velocity u z (z).Simulation results (symbols) are compared to the base state (lines).
FIG. 3. Base-state profiles of (a) active stress component p zz /nσ 0 , and (b) aerotaxic drift U d /U 0 .(c) Base-state values of β and Pe b and corresponding regimes in simulations, compared to the marginal stability curve β crit of Kasyap and Koch.23The data points correspond to values of L in the range 50-600 μm.