Role of hydrodynamic interactions in chemotaxis of bacterial populations

We briefly review the classic velocity jump model used in previous works such as [42]. Bacteria chemotax by changing their frequencies of tumbling in response to external signal gradients. Biologically, this can be thought of as searching for local surroundings with more favorable conditions. This is achieved through complicated intracellular signaling events triggered by receptor-ligand binding and unbinding and involves internal signal excitation, adaptation and amplification as well as modification of flagellar rotation. In past works, a tumbling rate dependent on the swimming direction has been implemented to enforce chemotatic movement (e.g. [43]). In this paper, we employ a simplified approach that assumes the tumbling rate $\lambda$ depends on the external signal gradient directly [25, 44, 45]. Specifically, we assume

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eqn:lambda} \lambda^i = \lambda_0\left(1 - \frac{\kappa{\nabla S({\bf x}^i,t) \cdot {\bf d}^i}}{\eta + | \nabla S({\bf x}^i,t)|}\right). \nonumber \end{align} \tag{ 1 }$$

Here $S({\bf x}^i,t)$ is the external signal concentration at the bacterium location ${\bf x}^i$. The sensitivity to the chemical gradient is captured in the dimensionless constant $\kappa \in [0,1]$ and the noise a bacterium may experience in trying to identify the correct chemical gradient is encoded into the constant $\newcommand{\e}{{\rm e}} \eta >0$. Both parameters help quantify the dependence of bacterial tumbling on the signal gradient. In particular, the noise captured by $\newcommand{\e}{{\rm e}} \eta$ could represent thermal noise in the signal passing through the fluid or a bacterium's inability to completely identify the signal gradient from local chemical measurements. A recent work studying E. coli shows that these sources of noise occur at the interface of the chemical ligands and chemoreceptors can interact with one another resulting in an affect on the chemotactic motility of the bacterium [46]. This work is the first to use the simple form given in (1) to model the modified tumbling rate and it has not been verified for rod-shaped cells prior to this work; however, it captures the essential response of a bacterium and is compatible with the individual-based modeling framework employed here.

The natural, unbiased tumbling rate is given by $\lambda_0$. To build intuition for the change in tumbling rate as a function of the chemical signal gradient, we consider the limiting cases. Observe that in the limit of zero noise $\newcommand{\e}{{\rm e}} \eta \to 0$, (1) reduces to $\lambda = \lambda_0\left(1 - \kappa {\bf u}_S \cdot {\bf d}^i\right)$ where the dot product gives the cosine of the angle between the unit vector in the direction of the chemical gradient ${\bf u}_S$ and the bacterium's current orientation ${\bf d}^i$. Thus, if the bacterium is swimming in the direction of the chemical gradient the tumbling rate reduces to near zero, $\lambda_0(1-\kappa)$, but if a bacterium is swimming orthogonal to the chemical gradient, then $\lambda \approx \lambda_0$ and the average time between tumbles increases to the natural, unbiased rate.

Classic works studying the tumbling of E. coli, experimentally observed that the tumbling rate should be proportional to the chemical gradient [47]. Furthermore, this work makes the stronger statement that once a bacterium is swimming in the direction of the chemical gradient the tumbling rate should be decreasing monotonically. The importance of the swimming direction and the magnitude of the chemical gradient are crucial for the individual-based modeling framework and are captured in (1). This is accounted for by observing that in the limit as $|\nabla S| \to \infty$ the tumbling rate converges to $\lambda_0(1 - \kappa {\bf u}_S \cdot {\bf d}^i)$ which enforces the boundedness of the response to the chemical gradient even if very large. Also, as $|\nabla S| \to 0$ the tumbling rate returns to a value close to the typical isolated constant rate $\newcommand{\e}{{\rm e}} \lambda_0(1 - \frac{\kappa}{\eta} \nabla S \cdot {\bf d}^i)$ where the tumbling rate is linearly dependent on the chemical gradient with stiffness $\newcommand{\e}{{\rm e}} \kappa / \eta$ (see [44, 48] for alternative modeling approaches where $\lambda$ is directly proportional to the gradient). This behavior illustrating the stiffness and boundedness of the chemotactic response directly relates to the same response in the chemotaxis flux in the continuum description (e.g. the flux-limited Keller–Segel model). In contrast with the classic Keller–Segel model, the flux-limited Keller–Segel model has solutions which do not blow up in finite or infinite time [49]. This is also related to a linear instability of the kinetic (continuum) chemotaxis model based on the velocity jump process. In relation to bacterial suspensions, the length of a run is decided by the stiff response to temporal sensing of chemical cues [50].

For the simulations presented throughout this work, we take $\lambda_0 = 1.0$ tumb s⁻¹, $\newcommand{\e}{{\rm e}} \kappa = .9, \eta = 5 \times 10^{-5}~\mu$M mm⁻³ for comparison with prior results of a classic velocity jump model [25]. To incorporate tumbling in the dynamic equations, it is assumed that during a tumbling event a bacterium chooses a random direction ${\bf d}'$ instantaneously after stopping with no directional persistence, according to the constant turning kernel

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eqn:tumb} T({\bf d}, {\bf d}') = 1/|\mathcal{S}^2| \nonumber \end{align} \tag{ 2 }$$

where $\mathcal{S}^2$ is the unit sphere and $|\mathcal{S}^2| = 4\pi$ as in [43]. While the classic VJP model (1) and (2) only allows bacteria to perform straight periods of runs followed by a reorientation through tumbling, the approach here will allow for more complex fluid-mediated movement. In particular, bacteria will propel themselves forward, be advected by the local fluid flow generated by other bacteria, and attempt to follow the local chemical signal gradient which is also advected by the fluid.

2.1.1. Bacterial motion

We summarize the remaining components of the thin film quasi-2D model used here and originally developed in [26, 40, 41]. Denote the position of the center of mass and the orientation direction of the ith bacterium by ${\bf x}^i \in R^2$ and ${\bf d}^i \in \mathcal{S}^1$. In addition, it is assumed that a bacterium propels itself in the direction ${\bf d}^i$ with relative swimming speed $V_0$. We let ${\bf v}({\bf x}^i)$ be the fluid velocity on the location of the ith bacterium, which is modeled as a Stokesian fluid defined in section 2.1.3. Moreover, we introduce a short-range, volume exclusion force ${\bf F}^i = \sum\nolimits_{j\neq i} {\bf F}_j^i$ acting on bacterium i by all other bacteria. The form of ${\bf F}_j^i$ is of a truncated Lennard–Jones type (for greater detail refer to [26]). In the present model, volume exclusion only occurs when two bacteria are very close to each other to prevent them from occupying the same region. An alternative approach accounting for similar forces using squirming spherical swimmers has been considered in [51]. For simplicity we assume any excluded-volume interaction acts in the direction of the centers of mass of each bacterium resulting in no effect on bacterial orientation. By neglecting inertia, we arrive at the following system describing bacterial motion during a run

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle d{\bf x}^i =\Big \{V_0{\bf d}^i + {\bf v}({\bf x}^i) + {\bf F}^i/ \zeta \Big \} dt, \qquad 1\leqslant i\leqslant N, \label{eqn:xdot} \nonumber \end{align} \tag{ 3 }$$

where $\zeta$ is the drag coefficient for the bacterium, representing an effective Stokes' Law relating a force to a velocity at low Reynolds number, and N is the total bacterium number.

The fluid flow can also change the direction of the bacterial movement. We assume that the bacteria are ellipsoids with major and minor axes denoted as a and b. The governing equation for ${\bf d}^i$ during a 'run' is given by Jeffrey's equation [52]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle d{\bf d}^i = \Big \{-{\bf d}^i \times \boldsymbol{\omega}({\bf x}^i) - {\bf d}^i \times \left[B{\bf d}^i \times \left({\bf E}({\bf x}^i)\right)\cdot {\bf d}^i\right]\Big \} dt.\displaystyle\label{eqn:ddot} \nonumber \end{align} \tag{ 4 }$$

Here $\boldsymbol{\omega}({\bf x}^i) := \nabla \times {\bf v}({\bf x}^i)$ and ${\bf E}^{\,j}({\bf x}^i) := \{\nabla_{\bf x} {\bf v}({\bf x}^i) +$$[\nabla_{\bf x} {\bf v}({\bf x}^i)]^T\}/2$ are the vorticity and rate of strain respectively of the flow induced by all other bacteria on the location of the ith bacterium, and $B = (a^2-b^2)/(a^2+b^2)$ is the Bretherton constant accounting for the bacterial shape (with B  =  0 for spheres and B  =  1 for needles). The derivation of equation (4) and further details of the physical implications can be found in [26, 29, 52, 53].

Incorporating the directional jumps due to tumbling into (4), we formally obtain

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle d{\bf d}^i &= \Big\{-{\bf d}^i \times \boldsymbol{\omega}({\bf x}^i) - {\bf d}^i \times \left[B{\bf d}^i \times \left({\bf E}({\bf x}^i)\right)\cdot {\bf d}^i\right] \Big\} dt\nonumber \\ &\quad+ (\xi^i - {\bf d}^i) dN^i(t; \lambda^i), \nonumber \end{align*}$$

where $\xi^i$ are independent random variables uniformly distributed on $\mathcal{S}^2$ and $N^i(t; \lambda^i)$ are independent inhomogeneous Poisson processes with intensity $\lambda^i$. In this work, the only random component of motion comes from the tumbling mechanism. In principle, one could consider random fluctuations, but these effects appear minimal compared with the swimming motion and tumbling.

2.1.2. Signal dynamics

The tumbling frequency of a bacterium is determined by the external signal $S({\bf x},t)$. In this paper, we consider the situation in which there is no external signal, but rather the signal is secreted by the bacteria and causes the bacteria to aggregate as demonstrated experimentally in [54]. The chemical signal is then advected by the fluid itself. Thus, the governing equation is a reaction–diffusion–advection equation for $S({\bf x},t)$ for ${\bf x} \in \Omega, t >0$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eqn:chem} \displaystyle\partial_t S + \nabla \cdot \left({\bf v} S\right) = \alpha\Delta S + \sum_{j=1}^N f_{sec}\delta({\bf x}-{\bf x}^{\,j}(t)) - \gamma S. \nonumber \end{align} \tag{ 5 }$$

Here $\alpha$ is the diffusion coefficient of the signal, f_sec is the secretion rate by a single bacterium, and $\gamma$ is a constant degradation rate. The advection term $\nabla \cdot \left({\bf v} S\right)$ describes the dispersal of the signal due to the fluid motion and can be simplified to ${\bf v} \cdot \nabla S$ due to the incom- pressibility of the fluid. The parameters are taken as $\alpha = 4 \times 10^{-4}$ mm² s⁻¹ [55], $f_{sec} = \frac{1}{54} \times 10^{-3}~\mu$M s⁻¹ per bacterium, and $\gamma = 1/3 \times 10^{-3}$ s⁻¹ similar to the parameters leading to aggregation in the purely chemotaxis model presented in [42, 56].

2.1.3. Fluid flow

Observe that the typical size of a bacterium under consideration (e.g. E. coli) is 1–2 $\mu$m and each bacterium is capable of swimming at speeds up to $20~\mu {\rm m}~{\rm s}^{-1}$ in isolation [26, 36]. Given these values the Reynold's number of the suspension is Re $\approx$$10^{-5} \ll 1$. In addition, the time scale for the flow to reach its steady state is much smaller than the characteristic time scale for bacterial swimming. Based on these considerations, the fluid is governed by a linear Stokes equation which must account for a bacterium confined to a thin layer of fluid. An explicit analytical expression for the fluid velocity generated by a Stokeslet (point force monopole) in a thin film of thickness w satisying Stokes equation has been derived by using the method of images in [57] and simplified to an expression for the flow in the xy-plane (z  =  0) in [58]. Following a similar procedure, the dipolar solution can be obtained by taking the flow due to two oppositely oriented force monopoles and letting the distance between them go to zero as in [26]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eqn:flsolnp} {\bf u}({\bf x},{\bf d}) = K_u\nabla^3[\log(r)] : {\bf d}{\bf d}, \nonumber \end{align} \tag{ 6 }$$

where $r = |{\bf x}| = \sqrt{x^2+y^2}$ is the distance from the dipole location to the origin in the xy-plane. Here the prefactor $K_u = \frac{U_0}{6\pi~\mu w}h(z)$ has units of length⁴/time while $\nabla^3[\log(r)] : {\bf d}{\bf d}$ has units of 1/length³. This pre-factor contains the physical quantities of the dipole moment U₀, the film thickness w, the ambient fluid viscosity $\mu$. Here $h(z)$ is a function of the depth and has a parabolic profile h(z)  =  c(z)w² where c is quadratic in z [57]. At z  =  0, $h(0) = \frac{9}{16}w^2$, where $w \sim 10~\mu$m. Thus, in the computational domain (xy-plane, z  =  0) $K_u = \frac{2U_0w}{32\pi\mu}$. As opposed to a completely two-dimensional formulation which has a fundamental solution ${\bf u}({\bf x}) \sim \ln(r)$, the quasi-2d suspension decays as 1/r³ due to the confinement of the suspension. This results in a three-dimensional suspension where the z-dimension is considered to be much smaller than the x and y dimensions, $z \ll x,y$, giving it a quasi-two-dimensional appearance (see figure 1). A much more detailed explanation of the model formulation for the fluid velocity can be found in the authors previous works [26, 40, 41].

Zoom In Zoom Out Reset image size

This allows for a two-dimensional simulation in which three-dimensional features are encoded into the coefficient through K_u. This approach compares qualitatively well with the streamlines observed in recent experiment (see figure 2 and [59]). Using the linearity of Stokes equation, we approximate the true fluid flow ${\bf v}({\bf x}, t)$ from N bacteria at location ${\bf x}$ by the superposition of thin film solutions (6) for a single dipole

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eqn:flowp} {\bf v}({\bf x},t) = \sum_{j = 1}^N {\bf u}({\bf x} - {\bf x}^{\,j},{\bf d}^{\,j},t). \nonumber \end{align} \tag{ 7 }$$

Zoom In Zoom Out Reset image size

By considering pairwise interactions between the force dipolar representations of the bacteria in the semi-dilute regime, the flow generated by all the bacteria is simply the superposition of the flow generated by each cell. Semi-dilute models combined with an excluded-volume potential as implemented in this work have been successfully used to study many concentration regimes (e.g. see [26, 41]). In those recent works among others, the model starts to break down at extremely high concentration (>70% area fraction) due to the dominance of near-field forces over hydrodynamic. We save the discussion of the high concentration regime and the resulting effect on the distribution of the chemical for appendix C, since it diverts from the focus of this work on the semidilute regime.

2.1.4. Advantages of modeling approach

To handle the chemical secretion by the bacteria one would need to solve the coupled ODE/PDE system for bacterial movement and signal dynamics. To simplify the numerical approach we implement a so-called semi-discrete algorithm where we continuously solve the ODE equations of motion (3) and (4), but use a finite difference approximation to solve the chemical PDE (5). The bacteria are free to move continuously about the computational domain, but the finite difference scheme relies on the chemical concentration at discrete grid points. In order for the scheme to work well we need to account for two things: (i) the distribution of secreted chemical to neighboring grid points and (ii) an accurate representation of the chemical gradient, $\nabla S$ needed for the biased tumbling rate $\lambda$ in (1).

To answer the first point, we consider the relative amount of chemical near each grid point. Since the bacteria can move freely throughout the computational domain, we assume the signal secreted at a given time appears at the bacterium's present location. In order to capture the true chemical concentration we cannot apply a constant concentration at each of the four nearest grid points, because the gradient used in the tumbling dynamic would be far from reality. Instead, we use the location of the bacterium to compute the appropriate weights on each of the four nearest grid points so that the total amount secreted is consistent with biological considerations. Therefore, the grid points nearest the bacterium receive most of the signal, which is closer to reality, see figure 3.

Zoom In Zoom Out Reset image size

Numerically the equation for the signal (5) is solved using an ADI (alternating direction implicit) method, which is a multi-dimensional version of the Crank–Nicholson method. It is easier to solve the problem numerically using the non-dimensional equation (A.3) with full derivation in appendix A. The position of a bacterium inside the computational domain is given by ${\bf x}^i = (x^i,y^i)$. At each time step a bacterium secretes a given amount of chemical signal $\zeta$. This must be interpolated into a local chemical concentration S, deposited at the neighboring discrete grid points (see figure 3). For a bacterium inside a square subdomain with top right vertex $(i,j)$ the total amount of chemical a given bacterium senses at its location is based on weights given to the four grid points

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \tilde{S}({\bf x}) = \frac{A_4}{A}S_{i-1,j-1} + \frac{A_3}{A}S_{i,j-1} + \frac{A_2}{A}S_{i-1,j} +\frac{A_1}{A}S_{i,j}. \nonumber \end{align*}$$

Here $\tilde{S}$ is the interpolated chemical concentration, A  =  dx² is the total area occupied by the bacterium, and the sub-area of each bacterium is A_k for $k = 1, 2, 3, 4$. The quantity A_k/A is the fraction of the square occupied by subdomain A_k. The complete step-by-step procedure can be found in [42], but here we provide the idea. The additional chemical secreted by a bacterium is added to the four neighboring grid points based on its location at each time step is $f_{sec} = A_4\zeta / A^2$ at $(i-1,j-1)$, $A_3\zeta / A^2$ at $(i,j-1)$, $A_2\zeta / A^2$ at $(i-1,j)$, $A_1\zeta / A^2$ at $(i,j)$, and 0 otherwise Here $\zeta$ is the amount of chemical secreted by a bacterium and $\zeta/A$ is the concentration of chemical secreted into the square subdomain containing the bacterium. This method has also been shown to be unconditionally stable for multi-dimensonal diffusion and heat condution problems as well as being second order in time and space (e.g. [62, 63]). It has also been shown to produce results consistent with experimental observation of E. coli aggregation as noted in [25, 64].

To address the second point we need an accurate numerical representation for the gradient of the signal, $\nabla S$. The biased tumbling rate depends crucially on the chemical gradient. Thus, we must extend beyond the methods outlined in [42] to consider a weighted gradient using the discrete grid points $\nabla \tilde{S} = (\tilde{S}_x, \tilde{S}_y)$. This is done by incorporating more neighboring bacteria and approximating the partial derivatives present in the gradient by central differences.

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \tilde{S}_x({\bf x},t) &:= \frac{A_4}{A}\left[\frac{S_{i,j-1}-S_{i-2,j-1}}{2dx}\right]+\frac{A_3}{A}\left[\frac{S_{i+1,j-1}-S_{i-1,j-1}}{2dx}\right]\nonumber \\ &+\frac{A_2}{A}\left[\frac{S_{i,j}-S_{i-2,j}}{2dx}\right]+\frac{A_1}{A}\left[\frac{S_{i+1,j}-S_{i-1,j}}{2dx}\right]\nonumber \\ \tilde{S}_y({\bf x},t) &:= \frac{A_4}{A}\left[\frac{S_{i-1,j}-S_{i-1,j-2}}{2dy}\right]+\frac{A_3}{A}\left[\frac{S_{i,j}-S_{i,j-2}}{2dy}\right]\nonumber \\ &+\frac{A_2}{A}\left[\frac{S_{i-1,j+1}-S_{i-1,j-1}}{2dy}\right]+\frac{A_1}{A}\left[\frac{S_{i,j+1}-S_{i,j-1}}{2dy}\right]. \nonumber \end{align*}$$

This has been successfully implemented in a similar case with regard to bacterial chemotaxis in [42]. Just as in the case of bacterium secretion, we want to give more weight in the central difference approximation to the derivative nearest to the bacterium, but all should contribute to the total discrete gradient of S. Thus, this spatial grid choice is important, if too fine then the local gradients are only determined by the chemical contributions of individual bacteria, but if too coarse then bacteria do not have enough time to leave their local subdomain before secreting again. The spatial grid used for the finite differences in the chemical equation in the simulations is $\newcommand{\e}{{\rm e}} 10~\mu {\rm m} \approx 5\ell_0$. Observe that the amount of chemical secreted needs to be scaled with the spatial grid to ensure the desired concentration among the local grid points.

One of the advantages of the individual-based modeling approach used throughout this work when compared to continuum modeling with PDEs is the ability to investigate individual bacterium trajectories in time. Before considering the macroscopic dynamics and pattern formation for the bacterial suspension, we investigate the effect of including hydrodynamics on individual bacterium motion. It is observed that for early times when the cells only begin to secrete chemical into the environment the trajectories in each model are basically the same. Since the chemical concentration is small the advection term in (5) has only a small contribution and the magnitude of S remains small in (1). Thus, the tumbling rates in each approach remain the same producing similar trajectories given the same initial data (see figures 4(a) and (b)). The effect of hydrodynamics is not prevalent until the post-aggregation phase when large groups of bacteria secrete chemical signal in one location. The advection term has the ability to distribute the signal to larger areas of the surrounding medium leading to larger area coverage by a bacterium governed by hydrodynamic forces (see figures 4(c), (d) and appendix C). At later times, in a chemical rich environment, we observe that individual dynamics can differ in terms of area covered. Also, the effective tumbling rate is reduced on average by $20\%$ when hydrodynamic effects are incorporated (1 s⁻¹ in the VJP model to .82 s⁻¹ in the hydrodynamic model on average).

Zoom In Zoom Out Reset image size

We begin by showing that starting from a uniform bacterial distribution with no external signal in the domain initially, bacteria can self-organize into aggregates in response to the self-secreted signal, and small aggregates can merge into large aggregates gradually. This is done by comparing the two models outlined above. Throughout this section, we compare a typical realization of each model as well as the average behavior. When a single realization is depicted, the same two simulations are used in all the figures.

Both models were simulated on a 1.5 mm by 1.5 mm domain with periodic boundary conditions, assuming no initial chemical signal, $S({\bf x},0) = 0$. The simulations contain N  =  7400 circular bacteria randomly placed in the computational domain, corresponding to an average volume fraction $\phi = 1.03\%$ (see figures 6(a),(d)). Both models predict that bacteria form aggregates within minutes with size comparable to observations in experiments (5–10 bacteria lengths) [42]. Figure 5 shows the early population dynamics typically found in each model. Aggregates start to form around t  =  7 min with diameter $150~\mu {\rm m}$–$250~\mu {\rm m}$ (figures 6(b) and (e)). By t  =  13 min, almost all aggregates are well-defined (figures 6(c) and (f)). This is consistent with the observation of complex fluid flows and fragmented aggregation regions observed for 'pushers' (propulsion actuated from behind) using a continuum approach in [37].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Each aggregate produces a large amount of signal which then attracts isolated bacteria and other aggregates towards it. As time progresses, all the aggregates begin to move as a collective unit and start to merge. Interestingly, the HI model predicts a faster merging process than the VJP model. To illustrate this, we calculated the distribution of the pairwise bacterium distances d_ij as a function of time. Figure 6 plots the heat maps of these distributions and displays the corresponding bacterial configuration at five moments in time for each model (HI versus VJP). The horizontal red stripes near the t-axes represent the characteristic distance within a single aggregate. The other horizontal stripes represent characteristic distances between aggregates. The vertical lighter branches that collapse into the horizontal red strips represent the merging of two or more aggregates. From figure 6(a), we see that in the HI model realization (figure 6(b)), the initial aggregates start to merge around t  =  7 min and by t  =  48 min the whole population collapses into a single aggregate. However, in the VJP model realization, the initial aggregates begin to merge around t  =  7 min, but the whole population only merges into a single aggregate around t  =  65 min.

In order to explain the enhanced merging speed in the hydrodynamic model we must understand what the flows around each bacterium look like. When two bacteria are close to each other and aligned side-by-side, they weakly attract each other due to hydrodynamic interactions, as suggested by the fluid streamlines in figure 2. It has been shown that when the bacterium volume fraction exceeds a threshold (∼20% for bacteria), this weak attraction leads to the onset of collective bacterial motion [26, 36].

Using a simple algorithm to count the number of aggregates by identifying groups of bacteria within one cell length of each other, we can then define a merge event when the cells of one cluster are within two bacteria lengths of another cluster for a significant period of time (on the order of minutes). The results show that with hydrodynamic interactions and the advection of the chemical by the fluid incorporated, the time between the merging of any two aggregates t_merge and the time to form a single aggregate t_agg are significantly less than in the pure chemotaxis model. For the hydrodynamic model $t_{agg} = 69 \,{\rm min} \pm 28 \,{\rm min}$ and the average time between any single merge event is $\langle t_{merge} \rangle = 6.6 \,{\rm min} \pm 2.8 \,{\rm min}$. For the pure chemotaxis model $t_{agg} = 115 \,{\rm min} \pm 12 \,{\rm min}$ and the average time between a single merge event is $\langle t_{merge} \rangle = 11.27 \,{\rm min} \pm 1.2 \,{\rm min}$. The discrepancy in merge times is due to the chemical signal advection in the fluid. This allows each bacterium to find a strong chemical gradient faster leading to an increased rate of aggregate formation and merging. Whereas in the pure chemotaxis case, clusters appear to merge only if they happen to be close to each other as they move in time. This is further explained by considering recent results on enhanced mixing and fluid instabilities due to the hydrodynamics flows generated by motile cells in suspensions observed in Saintillan et al [34, 65]. These striking physical features would not be present in the pure chemotaxis model where the dynamics of the fluid cells and associated instabilities are ignored.

Next, consider the evolution of the number of aggregates in time displayed in figure 7(a) for each model. If bacteria are within two bacterium lengths then they are said to belong to the same group. To identify all groups, we start with the first bacterium, label it group one, then label all bacteria within two bacterium lengths of that bacterium group one. We then repeat this process until there is no other bacteria that qualify to be in group one. Then we move to the next bacterium that is not associated with group one, label it group two, and repeat the above process to find all groups. Groups consisting of more than $5\%$ of the bacteria population are referred to as aggregates or clusters. In figure 7(a), the lines represent the average over ten realizations of each model and the error bars indicate the standard deviation. We note that the two curves agree very well for the first 10 min, but as time progresses the HI model predicts faster merging of these aggregates. In the longtime limit in both cases there should be one large aggregate, but the intermediate behavior is different.

Zoom In Zoom Out Reset image size

To demonstrate whether collective motion occurs during aggregate merging, we compute the peak volume fraction over the entire computational domain as a function of time for each model. To calculate the peak volume (indeed, in 2D the area) fraction, we divide the whole domain into boxes of dimension $60~\mu$m $\times$ 60 $\mu$m and calculate the average volume fraction inside each box. We then take the maximum volume fraction over all the boxes and call it the peak volume fraction at that time. Figure 7(b) shows the statistics for both models suggesting that the peak volume fraction predicted by the two models agrees well for t smaller than 10 min, but the HI model predicts a larger peak volume fraction for larger time, due to the faster merging of aggregates. Intriguingly, the peak volume fraction predicted by both models eventually exceeds the collective motion threshold and the behavior of the two models start to differentiate after the threshold is achieved. This is due to the fact that the hydrodynamic interactions and the fluid advection enhance the alignment and collective swimming of nearby bacteria and this effect is not present in the pure chemotaxis model. In addition, in each case, the peak volume fraction approaches a constant value in time indicating the single aggregate roughly maintains a constant volume fraction. This is consistent with the previous work in [32, 37] where the flow instability generated from hydrodynamic interactions suppressed growth in concentration.

To further understand the role of collective motion in aggregate merging, we simulate a case with much fewer bacteria in the domain. The results with an average volume fraction ten times smaller than the simulations in figure 7 are plotted in figure 8. We see that the statistics of aggregate number and peak volume fraction predicted by the two models match well in general and for each model the threshold for collective motion was never reached. This suggests that as long as the bacterial volume fraction is below the threshold for collective motion the hydrodynamic model and the VJP model demonstrate similar behavior in aggregate formation and merging. Further, the computational complexity involved with the hydrodynamic model is not needed.

Zoom In Zoom Out Reset image size

We next investigated the population behavior after all the bacteria merge into a single aggregate using both the HI model and the VJP model. In particular, figures 9(a)–(c) shows the shape of the final aggregates predicted by the two models. The HI model leads to an amorphous aggregate that can change shape over time due to the advection of the chemical by the surrounding fluid and the tendency of the bacteria to swim toward high concentrations of chemical collectively when the local density is large. The VJP model predicts a radially symmetric aggregate in which bacterium motion is solely governed by the local chemical gradient. The agent-based model used throughout this work explains the discrepancy through the hydrodynamic interactions. Unlike in the VJP model, when hydrodynamic interactions are present we observe that local groups of swimmers may move together in a direction other than the cluster center due to the advection of the chemical by the fluid before being attracted back by the chemical gradient. This is clearly observed in figure 9(b) and is only possible with the individual-based modeling approach. This gives insight into the amorphous shape; namely, the chemical concentration distribution is also non-symmetric and moving away from the aggregate center in some directions. Thus, portions of the aggregate move toward these higher concentrations creating the non-symmetric shapes and 'appendages' that grow out from the aggregate center in time. This is not possible in the VJP model because the chemical is not advected by the fluid and simply remains primarily at the center of the cluster forming the outward symmetric shape.

Zoom In Zoom Out Reset image size

This effect can only partially be explained by observations from the continuum model in [37] which shows that aggregates leading to high volume fractions create a destabilizing active stress resulting in unsteady flow. This is consistent with figure 9 which shows the high density cluster partially splitting into arms and recombining due to chemotaxis (see also the supplemental movies). For an explanation of why this occurs we return to single bacterium dynamics as discussed in section 4.1 where the advection term in (5) in a chemical signal rich environment leads to drifting trajectories of single cells away from the aggregate. This may also contribute to the amorphous shape and arms extending out from the aggregate body in the hydrodynamic case. The size of the aggregates predicted by both models are about the same, indicating that the intrinsic aggregate size is dictated by chemotaxis or the physical size of the bacterium.

Once all the bacteria form a final aggregate, they can move as a single entity. This has also been observed in recent experiments by Saragosti et al where chemotactic bacteria form traveling bands or rings swimming collectively toward high signal concentrations [66, 67]. To track this collective movement, we recorded the center of mass and investigated whether the aggregate exhibits any type of diffusive behavior. Specifically, we computed the diffusion coefficient D and exponent $\alpha$ from mean-squared displacements $\langle |{\bf r}|^2 \rangle = 4Dt^{\alpha}$ in figure 10. For the HI model, the diffusion exponent is 1.2–1.3 and the diffusion coefficient of $2.38\times 10^{-4}~{\rm mm^2~s^{-1}} = 2.38 \times 10^{-6}~{\rm cm^2~s^{-1}}$. These values indicate super-diffusive behavior as suggested by experiments [68, 69]. In contrast, the VJP model seems to predict a much smaller diffusion exponent $\alpha = .60$ and coefficient $D = 1.55 \times 10^{-4}~{\rm mm~s^{-1}}$. This confirms that the VJP model is not capable of exhibiting collective swimming even in dense aggregates and the aggregates themselves do not display any super-diffusive behavior. Figure 11 plots the aggregate trajectory during one hour of time for one particular realization for each model. Clearly we see that the HI model predicts a more profound movement of the aggregate over time and the VJP model predicts a stationary aggregate located at the highest concentration of chemical secreted.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The baseline time step is .1sec and the spatial step is $10~\mu$m (on the order of a bacterium length). Additional realizations under the same initial conditions were run with spatial and time steps twice as large and half as large. The results show no large qualitative or quantitative difference in average peak volume fraction (range from .521–.571), transition time to cluster formation (5.2–10.1 min), and peak cluster number (8–10 initial clusters). For a thorough sensitivity analysis of the effect of the model parameters on the results please see appendix B.

The microscopic modeling approach here has provided new biological insight into the aggregation formation process and resulting dynamics. First, recent experiments have determined that chemotactic bacteria can swim collectively along a concentration gradient [66]. The work here shows that in the absence of hydrodynamic interactions the aggregates become stationary and radially symmetric. Also, we observe that hydrodynamic flows lead to instability in the aggregate shape allowing a segment of the aggregate to extend and move the cell in a given direction. Observations of the specific configuration of dense aggregates are not possible with a continuum modeling approach. The microscopic model here also provides more insight into why the peak volume fraction becomes constant in the single aggregate setting. While models with and without hydrodynamic interactions can lead to temporary spikes in bacterial concentration at locations of high chemical signal, the excluded-volume constraints cap the maximum volume fraction at the packing fraction for circles in two dimensions (bacterial peak volume fraction approaches $\Phi = .7$–.8 and the maximum packing fraction for spheres is $\Phi = .907$). The excluded-volume constraints ensure the local volume fraction does not exceed what can be physically observed. This constraint is not needed in the absence of hydrodynamic interactions because in the chemotaxis only model there is no attractive force between cells. Thus, in time the volume fraction will still approach a constant value when a single aggregate is present, but this is only a fraction of the value seen in the hydrodynamic model. Recent experiments confirm that excluded-volume effects play a large role in the properties of the collective mesoscale structures emerging with the clusters [26, 36].

In addition, the simulations show that the results seen in this paper such as the enhanced clustering and amorphous shape are possible when the fluid advection term has a greater effect relative to the diffusion term. This is most easily studied in the non-dimensional equation for the chemical concentration derived in appendix A. If the effective diffusion coefficient, $\alpha^*$, is greater than one then the chemical concentration distribution is dominated by random diffusion and the added complexity of the hydrodynamic model is not needed. However, if $\alpha^* < 1$, then the fluid advection of the chemical is the dominant interaction and the vortices generated by collective cell swimming allow the chemical to collect in local areas forming the clusters faster and the unique non-symmetric shape of the cluster. This observation is new to this work and only possible with the individual-based modeling approach introduced and analyzed in this work.

Finally, one of the novel results here is the investigation of the merging times of aggregates. The model clearly shows the hydrodynamic flows combined with the advection of the chemical signal lead to significantly faster merging times of aggregates (around twice as fast). All these results suggests the purely chemotactic aggregation, while appearing similar at a snapshot of time is fundamentally different than hydrodynamically-induced autochemotactic aggregation. Additionally, the fact that results are similar at short times suggests that the chemotactic interaction initially forms the aggregates, but the hydrodynamic interactions take over resulting in collective swimming and instabilities in the aggregate state.

Though the model behaviors are qualitatively different, this work provides testable predictions for future experiment focused on isolating the effects of the hydrodynamic interactions. In particular, observations should focus on aggregate density, time to formation, and dynamics. We hypothesize that the hydrodynamic model is closer to actual experiment because it allows for collective swimming at larger concentrations which has been observed (e.g. [36, 60]) whereas the non-hydrodynamic model only allows for cluster formation. We also note the model and results can be extended to study Janus particles in comparison with [10] where the asymmetry of the catalytic cap can lead to a different interplay between chemotactic and hydrodynamic effects. This work lies at the interface of mathematics, physics, and biology providing biophysical insight into aggregate formation with agent-based microscopic models that are amenable to mathematical analysis and simulation.