Fluid dynamics of bacterial turbulence

Self-sustained turbulent structures have been observed in a wide range of living fluids, yet no quantitative theory exists to explain their properties. We report experiments on active turbulence in highly concentrated 3D suspensions of Bacillus subtilis and compare them with a minimal fourth-order vector-field theory for incompressible bacterial dynamics. Velocimetry of bacteria and surrounding fluid, determined by imaging cells and tracking colloidal tracers, yields consistent results for velocity statistics and correlations over two orders of magnitude in kinetic energy, revealing a decrease of fluid memory with increasing swimming activity and linear scaling between energy and enstrophy. The best-fit model parameters allow for quantitative agreement with experimental data.

Self-sustained turbulent structures have been observed in a wide range of living fluids, yet no quantitative theory exists to explain their properties. We report experiments on active turbulence in highly concentrated 3D suspensions of Bacillus subtilis and compare them with a minimal fourth-order vector-field theory for incompressible bacterial dynamics. Velocimetry of bacteria and surrounding fluid, determined by imaging cells and tracking colloidal tracers, yields consistent results for velocity statistics and correlations over two orders of magnitude in kinetic energy, revealing a decrease of fluid memory with increasing swimming activity and linear scaling between energy and enstrophy. The best-fit model parameters allow for quantitative agreement with experimental data.  [1][2][3][4][5][6][7][8][9][10] has shed light on generic ordering principles that appear to govern collective dynamics of living matter [11][12][13][14][15], from large-scale animal swarming [1,2] to meso-scale turbulence in microbial suspensions [3][4][5][6][7][8] and micro-scale selforganization in motility assays [9,10]. Although very different in size and composition, these systems are often jointly termed 'active' fluids, for which there is now a range of continuum theories [12,[14][15][16][17][18][19][20][21][22]. From these have come important qualitative insights into instability mechanisms [13-16, 20, 23] driving dynamical pattern formation, but a quantitative picture remains inchoate; even for the simplest active (e.g., bacterial or algal) suspensions uncertainty remains about which hydrodynamic equations and transport coefficients [24,25] provide an adequate minimal description, due in large part to the inability of existing data to constrain the manifold parameters in these models. One approach to remedy this problem is to characterize collective dynamics as in high Reynolds number fluid turbulence, in terms of kinetic energy, enstrophy and spatio-temporal correlation functions, and to compare with an appropriate longwavelength theory (i.e. Navier-Stokes-type equations). We present such an analysis here, measuring collective behavior in dense suspensions of the bacterium Bacillus subtilis in comparison to predictions of a (fourth-order) continuum model for bacterial flow [7,26].
Previous experimental studies of bacterial suspensions in open droplets [3,4,27,28], freestanding films [5,8,25,29], on surfaces [6,30,31], or quasi-2D microfluidic chambers [7] focused separately on the bacterial and fluid components, leaving uncertain how accurately passive tracers [32,33] reflect collective bacterial dynamics. The experiments reported here, performed in closed 3D microfluidic chambers, allowed near-simultaneous mea-surements of cell and tracer motion, and exploit a natural reduction in bacterial swimming activity due to oxygen depletion [8,27,34] to obtain data spanning two orders of magnitude in fluid kinetic energy. Combined with extensive 3D numerical simulations of the model, this data allows robust parameter estimates. Quantitative agreement between experiment and theory suggests that this model presents a viable generalization of the Navier-Stokes equations to incompressible active fluids.
Wild-type strain 168 of B. subtilis has cigar-shaped cell bodies on average 0.8 µm in diameter and 5 µm long [7]. It was streaked on LB medium agar plates from frozen stocks. Colonies from these plates were used to inoculate overnight cultures in Terrific Broth (TB; Sigma), which were back-diluted 1:100 into 100 ml of TB and grown to mid-logarithmic phase on an orbital shaker at 37 • C. These cultures were then concentrated 400× at 4, 000 g (final volume fraction ∼ 50%), and fluorescent microspheres (diameter 1 µm, F-8816, Invitrogen) were added at a final concentration of ∼10 9 beads/ml. The resulting suspensions were loaded into polydimethylsiloxane (PDMS) microfluidic devices, consisting of a series of cylindrical chambers (radius 750 µm, height 80 µm), connected by thin channels [7,36]. The inlet and outlet of the device were sealed with vacuum grease, and images were acquired in the (xy)-midplane of the chambers, ≈40 µm above the bottom, using a Zeiss 40× (NA 1.3) oil immersion objective and a high-speed camera at 40 fps (Fastcam SA-3, Photron). Movies were recorded in pairs for each field of view (768 × 800 pix; 1 pix = 0.36 × 0.36 µm 2 ), one with bright-field illumination and one with fluorescence excitation by a 633 nm laser (B&W Tek) at ∼20 mW. These movies were taken immediately after each other with a ∼3 min time lag between subsequent pairs. During the ∼10 min imaging period for each  Table I). Scale bars 70 µm.
device, the motility of B. subtilis cells decreased markedly due to oxygen depletion [27]. The experimental setup yields 2D projected velocities of 3D suspension motion ( Fig. 1). Data were analyzed under the assumption that the flow structures are isotropic, as verified by test measurements at different distances from the chamber bottom. Commercial particle tracking velocimetry (PIV) software (Dantec Flow Manager) was used to determine the bacterial flow velocity (v x , v y ) from bright-field images ( Fig. 1a,b), corrected for systematic pixel-locking errors [28]. Data shown in Figs. 2 and 3 are based on 7 movie segments (40 fps, each 50 s long) corresponding to 7 different activity levels.
Global bacterial flows were quantified by the in-plane kinetic energy E xy (t) = (v 2 x + v 2 y )/2 and in-plane en- is the vertical component of vorticity and · is a spatial average. While E xy and Ω z fluctuate, their time averages (E xy , Ω z ) are approximately constant during the 50 s time interval used in the data analysis (Fig. 2b,c). Over two orders of magnitude in energy (Fig. 2d) we observe the linear scaling Ω z = E xy /Λ 2 , with Λ ≈ 24 µm being roughly one half of the typical vortex radius.
Probability distribution functions (PDFs) of the in-plane bacterial velocity are approximately Gaussian, with a slight broadening due to collective swimming (Fig. 2a). The negative values of the equal-time spatial velocity correlation function (VCF; Fig. 3a) indicate the existence of vortices [4] (Fig. 1). The VCF is remarkably robust with respect to changes in the bacterial activity; in particular, the typical vortex radius R v ∼ 40 µm, estimated from the first zero of the VCF, depends only weakly on the kinetic energy. This result is consistent with recent findings by Sokolov and Aranson [8] for free-standing films. The vortex size in 3D is roughly five times larger than for quasi-2D turbulence in thin microfluidic chambers [7], where bacterial swimming and hydrodynamic interactions are suppressed by the nearby no-slip boundaries [36,37]. Unlike the spatial VCF, the two-time velocity auto-correlation function (VACF) varies systematically with energy or vorticity (Fig. 3b), but they collapse when plotted as functions of the dimensionless lag-parameter τ Ω 1/2 z (inset of Fig. 3b), implying that the higher the activity the shorter the memory of the bacterial fluid. Generally, the statistics of 3D bacterial turbulence differ strongly from conventional 3D Navier-Stokes turbulence [38,39], as bacteria inject energy on the smallest scales, inducing an 'upward' energy cascade towards larger length scales.
We infer the flow of the solvent medium from particle tracking velocimetry (PTV) analysis of the fluorescence images, which only show the tracer particles, assuming that they are passively advected. Data shown in Figs. 2 and 3 are based on 7 movies (40 fps, length 100 s) at different activities. Trajectories of individual tracer particles were found with a custom algorithm which, depending on seeding density and tracer dynamics, was able to identify up to 10 4 in-plane tracks, the longest typically lasting 5 − 8 s. The effective sample size was insufficient to determine reliably the tracer VACFs, but did yield global flow properties, velocity histograms and equaltime VCFs. The velocity PDFs, calculated directly from individual tracer velocities, are approximately Gaussian with a peak at small velocities from tracer accumulation near the vortex centers (Fig. 2a).
Estimates from PTV for the medium VCF and enstrophy were obtained by interpolating tracer velocities on a 450 × 450 pix subwindow in the center of the imaging plane using MATLAB's Delaunay triangulation with a lattice spacing ∆ = 90 pix/N f , where N f is the mean number of tracers detected per frame. The accuracy of this reconstruction procedure is controlled by the tracer concentration, which was kept low to limit effects on the bacteria motion and to avoid tracking ambiguities (typically N f ∈ [47,144] for data shown in Figs. 2 and 3). As a result, the uncertainties for the PTV data are considerably larger than for PIV data (see Fig. 2d). The interpolated tracer flow fields were used to estimate the kinetic energy E xy , enstrophy Ω z , and spatial correlation functions of the in-plane medium flow components. In agreement with the PIV results for the bacterial flow, we find again a linear enstrophy-energy relation (Fig. 2d) and comparable vortex radii, using the first zero of VCF as an estimate (Fig. 3a). We may therefore conclude that, at our very high bacterial concentrations, solvent and bacterial flow statistics become tightly linked.
We now examine how these data compare to predictions of a theory of active fluids introduced recently [7,26]. This minimal continuum model assumes that, at high concentrations, the bacterial flow due to swimming and advection can be described by a single velocity field v(t, x) and a pressure p(t, x). They obey the incompressibility condition ∇ · v = 0 and Equation (1) extends the incompressible Toner-Tu theory [14,18,40] with a fourth order term as in the Swift-Hohenberg equation [41]. The parameter λ 0 describes ad-vection and nematic interactions, and λ 1 an active pressure contribution [26]. For pusher swimmers [36] like B. subtilis, general considerations of hydrodynamic [42] and nematic stresses [26,43] suggest that λ 0 ≥ 1 and λ 1 (λ 0 − 1)/3 ≥ 0 in 3D. The (β, v 0 )-terms correspond to a quartic Landau-type velocity potential [14,18,40] and are physically motivated by the observation of extended jet-like streaming regions in B. subtilis suspensions at intermediate concentrations [28]. The parameter v 0 defines the collective speed that would be achieved if all bacteria were to move in the same direction. When β = 0 the model does not conserve momentum or energy, as it describes exclusively the bacterial flow component, which may exchange energy and momentum with the solvent. The nonlocal (Γ 0 , Γ 2 )-terms encode passive and active stresses due to hydrodynamic and steric interactions. For λ 0 = 1, λ 1 = β = Γ 2 = 0 and Γ 0 > 0, the model reduces to the incompressible Navier-Stokes equation. A detailed stability analysis [26] shows that when λ 0 = 0, β > 0, v 0 > 0, Γ 2 > 0 and Γ 0 < 0 this is one of the simplest vector models to describe phenomenologically the formation of jets and turbulent vortices in quasiincompressible active suspensions. Very recently, the 2D version of Eq. (1) has been shown to provide a quantitative mean field description of bacterial meso-scale turbulence in quasi-2D suspensions [7]. Its applicability to the physically more relevant 3D case is first explored here. We simulated Eq. (1) in 3D with periodic boundary conditions using a pseudospectral operator-splitting algorithm [44,45] and a pressure correction subroutine to ensure incompressibility [7,26]. Simulation grids ranged from 128 3 lattice points for parameter pre-screening to 256 3 for statistical analysis. Numerical stability of the solver was verified for a wide range of parameters and (color online) Correlation functions for solvent (PTV) and bacterial (PIV) flow at different energies and bestfit continuum theory (see Table I), using the same colors as in . Since Ωz ∝ Exy (Fig. 2d), this implies that the higher the bacterial activity the shorter the flow memory. space-time discretizations. All simulations were initiated with randomly chosen velocities. Figure 4 shows structure-formation in a typical simulation domain.
Since in 3D we have λ 1 (λ 0 − 1)/3 [26], Eq. (1) has essentially five free parameters (λ 0 , β, v 0 , Γ 0 , Γ 2 ). Two of those can be eliminated by choice of appropriate length and time units. We adopt a natural unit system such that the vortex wave-length scale Λ Γ = 2π Γ 2 /(−Γ 0 ) = 2π and v 0 = 1. In our simulations, the box length is fixed as L = 12Λ Γ , corresponding to approximately twice the experimental field of view, and the time step as ∆t = 0.05Λ Γ /(2πv 0 ). To estimate the three remaining parameters (λ 0 , β, Γ 0 ), we note that Γ 0 and Γ 2 define a typical vortex speed V Γ = −Γ 3 0 /Γ 2 . In the turbulent regime, it is plausible that V Γ is smaller than but close to v 0 , i.e. V Γ = ζv 0 where ζ 1. Furthermore, for pushers, the dimensionless parameter λ 0 should be larger than 1, but smaller than for quasi-2D suspensions [7], since nematic (steric) stresses can be more easily avoided in 3D; we infer λ 0 ∼ 2. Finally, the acceleration time scale τ 0 = (βv 2 0 ) −1 should be of the order of the vortex time-scale Λ Γ /V Γ . Using these estimates as initial values in a systematic parameter scan, and by comparing with the bacterial PIV data, we obtained the best-fit parameters in Table I. Generally, the VCFs and VACFs respond sensitively to parameter variations in the simulations, suggesting that the estimates in Table I are accurate within 10-15% for quasi-incompressible B. subtilis suspensions. As an independent cross-check, we computed Λ = (E xy /Ω z ) 1/2 from the best-fit simulation using Λ Γ ∼ 50 µm and found Λ ∼ 29 µm which compares well with the experimental PIV value in Fig. 2d. We stress that the conserved form of the bacterial velocity PDFs (Fig. 2a), VCFs, and VACFs (Fig. 3) implies that all our experiments can be fitted by a single set of rescaled parameters (λ 0 , β, Γ 0 ), as it suffices to adjust the physical  I: Parameters of the best-fit continuum model. To match a specific experiment, one must merely adjust the physical value of v0 by equating Exy = 0.54v 2 0 to the corresponding kinetic energy value in Fig. 2d. values of v 0 and Λ Γ to match the kinetic energy and vortex length at a given bacterial activity level. As evident from the flow patterns in Fig. 1 and from the solid curves in Figs. 2a and 3, the best-fit parameters yield good qualitative and quantitative agreement with the experiments.
For incompressible 'passive' fluids, that are governed by the Navier-Stokes equations transport parameters have of course been measured for a wide range of materials [46]. In contrast, quantitative theories of even the simplest active fluids have been lacking. We have shown here that the minimal fourth-order vector model [7,26] in Eq. (1) reproduces the main statistical features of selfsustained 3D bulk turbulence in concentrated bacterial suspensions, suggesting that this theory is a viable candidate for the quantitative description of incompressible active fluids. Due to the close correlation between bacterial and medium (tracer) flow observed in our experiments, we expect that this generic model will be useful in a wide range of future applications, in particular for predicting the effects of confining geometries on collective microbial dynamics [47,48] and for understanding the anomalous viscosities of active fluids [24,25].