Spatial modulation of individual behaviors enables an ordered structure of diverse phenotypes during bacterial group migration
Abstract
Coordination of diverse individuals often requires sophisticated communications and highorder computational abilities. Microbial populations can exhibit diverse individualistic behaviors, and yet can engage in collective migratory patterns with a spatially sorted arrangement of phenotypes. However, it is unclear how such spatially sorted patterns emerge from diverse individuals without complex computational abilities. Here, by investigating the singlecell trajectories during group migration, we discovered that, despite the constant migrating speed of a group, the drift velocities of individual bacteria decrease from the back to the front. With a Langevintype modeling framework, we showed that this decreasing profile of drift velocities implies the spatial modulation of individual runandtumble random motions, and enables the bacterial population to migrate as a pushed wave front. Theoretical analysis and stochastic simulations further predicted that the pushed wave front can help a diverse population to stay in a tight group, while diverse individuals perform the same type of mean reverting processes around centers orderly aligned by their chemotactic abilities. This mechanism about the emergence of orderly collective migration from diverse individuals is experimentally demonstrated by titration of bacterial chemoreceptor abundance. These results reveal a simple computational principle for emergent ordered behaviors from heterogeneous individuals.
eLife digest
Organisms living in large groups often have to move together in order to navigate, forage for food, and increase their roaming range. Such groups are often made up of distinct individuals that must integrate their different behaviors in order to migrate in the same direction at a similar pace. For instance, for the bacteria Escherichia coli to travel as a condensed group, they must coordinate their response to a set of chemical signals called chemoattractants that tell them where to go.
The chemoattractants surrounding the bacteria are unequally distributed so that there is more of them at the front than the back of the group. During migration, each bacterium moves towards this concentration gradient in a distinct way, spontaneously rotating its direction in a ‘runandtumble’ motion that guides it towards areas where there are high levels of these chemical signals. In addition to this variability, how well individual bacteria are able to swim up the gradient also differs within the population. Bacteria that are better at sensing the chemoattractant gradient are placed at the front of the group, while those that are worst are shifted towards the back. This spatial arrangement is thought to help the bacteria migrate together. But how E. coli organize themselves in to this pattern is unclear, especially as they cannot communicate directly with one another and display such diverse, randomized behaviors.
To help answer this question, Bai, He et al. discovered a general principle that describes how single bacterial cells move within a group. The results showed that E. coli alter their runandtumble motion depending on where they reside within the population: individuals at the rear drift faster so they can catch up with the group, while those leading the group drift slower to draw themselves back. This ‘reversion behavior’ allows the migrating bacteria to travel at a constant speed around a mean position relative to the group.
A cell’s drifting speed is determined by how well it moves towards the chemoattractant and its response to the concentration gradient. As a result, the mean position around which the bacterium accelerates or deaccelerates will vary depending on how sensitive it is to the chemoattractant gradient. The E. coli therefore spatially arrange themselves so that the more sensitive bacteria are located at the front of the group where the gradient is shallower; and cells that are less sensitive are located towards the back where the gradient is steeper.
These findings suggest a general principle for how bacteria form ordered patterns whilst migrating as a collective group. This behavior could also apply to other populations of distinct individuals, such as ants following a trail or flocks of birds migrating in between seasons.
Introduction
Collective group migration, as an important class of coordinated behaviors, is ubiquitous in living systems, such as navigation, foraging, and range expansion (Krause et al., 2002; Partridge, 1982; Sumpter, 2010). In the presence of individual heterogeneity, the migrating group often exhibits spatially ordered arrangements of phenotypes (Krause et al., 2002; Parrish and EdelsteinKeshet, 1999; Partridge, 1982; Sumpter, 2010). In animal group migration, individual behavioral abilities (e.g., directional sensitivity) would result in social hierarchy, which further drives the spatial arrangement in a coordinated group (Couzin et al., 2005). At the same time, spatial arrangements can lead to different costs and benefits for the individuals participating in the group migration (Krause, 1994; Parrish and EdelsteinKeshet, 1999; Partridge, 1982). Participating individuals must follow disciplinary rules to organize themselves into coordinated patterns while on the move, which requires complex computational abilities to interact with the group and the environment (Couzin and Krause, 2003; Couzin et al., 2002; Vicsek and Zafeiris, 2012). Therefore, understanding how individuals of different phenotypes determine their location in the group is an essential prerequisite to uncover the organization principles of collective populations.
The chemotactic microbe, Escherichia coli, provides a simple model to address the emergence of collective decisionmaking among diverse population, as it can exhibit both individualistic behaviors (Dufour et al., 2016; Frankel et al., 2014; Kussell and Leibler, 2005; Waite et al., 2016; Waite et al., 2018) and collective migratory patterns (Adler, 1966a; Fu et al., 2018; Keller and Segel, 1971a, Wolfe and Berg, 1989). Individual cells perform runandtumble random motions by spontaneously switching the rotating direction of flagella (Berg, 2004; Berg and Brown, 1972). These cells can facilitate the chemotaxis pathway to control the switching frequency to bias their motions toward their favorable direction along the chemoattractant gradient, where the efficiency to climb the gradient is defined as the chemotactic ability ($\chi $) (Celani and Vergassola, 2010; Dufour et al., 2014; de Gennes, 2004; Si et al., 2012). In addition, the chemotactic abilities of individual cells exhibit substantial phenotypic heterogeneity even for the clonal bacterial population, which diversifies the chemotactic response to identical signals (Spudich and Koshland, 1976; Waite et al., 2016; Waite et al., 2018). Despite the stochastic solitary behavior and variations in phenotypic ability, the E. coli population can migrate as a coherent group by following a selfgenerated attractant gradient, via consumption of the whole population (Adler, 1966a; Saragosti et al., 2011; Wolfe and Berg, 1989). The migratory group form a stable pattern of phenotypes sorted by their chemotactic abilities (Figure 1A), which is believed to maintain phenotypic diversity in the group (Fu et al., 2018; Waite et al., 2018). Although a previous study showed that behavior modulation helps migrating bacteria to maintain a consistent group (Saragosti et al., 2011), how individuals with phenotypic variations manage to determine their relative positions in the group remains to be determined.
Here, we analyzed singlecell trajectories of bacterial runandtumble motions in the chemotactic migration group (see Materials and methods). We found that the expected drift velocity of individual cells decreased from the back to the front. Such a spatial profile modulates cells to behave as mean reverting processes relative to the entire group, that is, cells effectively tend to revert their direction of runs toward the mean position of the group. Using an OrnsteinUhlenbeck (OU)type model, we demonstrated that the mean reversion behavior is a result of a pushed wave, where the driving force decreases from the back to the front of the group. Cells of different phenotypes are imposed to the same type of driving force, of which the strength is coupled with their chemotactic abilities. As a result, the pushed wave front, driven by the spatially structured force, can maintain more diverse individuals in the migratory group. By theoretical analysis and stochastic simulations, we also discovered that the balanced locations of diverse phenotypes are spatially ordered by their chemotaxis abilities. Further simulations and experiments with cells of titrated chemoreceptor abundance demonstrated that this spatial modulation of individual behavior enables the ordering of bacteria with diverse chemotaxis abilities. Therefore, although the highorder computational abilities are not available to simple organisms, the spatial modulation of stochastic behaviors at the individual level reveals novel decisionmaking capabilities at the population level.
Results
The drift velocity of individual cells during group migration exhibits a spatially structured profile
To directly investigate the ordering mechanism in a coherent migration group, we quantified the stochastic behaviors of bacterial runandtumble motions relative to the stable migrating group. To achieve this, we employed a microfluidic device that generated a stable propagating band of bacteria, as previously reported (Fu et al., 2018; Saragosti et al., 2011). Using aspartate (Asp) as the only chemoattractant to drive the migration of E. coli, we tracked a small fraction of fluorescent bacteria (JCY1) as representatives of the nonfluorescent wildtype cells (RP437) (see Materials and methods, Appendix 1—figure 1, Video 1). Because the group velocity ${V}_{G}=\langle \Delta {x}_{i}\left(t\right)/\Delta t\rangle $ was constant over time, $V}_{G}\sim 3.0\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}/\mathrm{s$ (Appendix 1—figure 2), we were able to map the tracks to a moving coordinate $z=x\text{}\text{}{V}_{G}t$.
With the shifted tracks, we calculated the key statistics of the singlecell behaviors. We first checked that the average instantaneous velocity, ${V}_{I}(z)=\u27e8\mathrm{\Delta}{x}_{i}(z)/\mathrm{\Delta}t\u27e9$, was constant along the density profile $\rho \left(z\right)$ (Figure 1B). This result confirmed that the group migrates coherently. Then, we identified all the run states and tumble states of individual trajectories using a previously described computer assistant program (Dufour et al., 2016; Waite et al., 2016), to ensure that the tracks are separated into successive tumbleandrun events (see Materials and methods). Comparing the sample events initiated from the back (b), middle (m), and front (f) of the migration group, we observed that the runs in the front were longer but distributed more uniformly in terms of the directionality, whereas the runs in the back were shorter but more oriented toward the group migration direction (Figure 1C, Appendix 1—figure 3A, B). Quantitatively, the statistics of run length (and duration) displayed exponential distributions with the means for the direction of group migration longer than those for the opposite direction (Figure 1D, Appendix 1—figure 3C). The angular distribution of the run length confirmed the difference in directionality between cells in different spatial locations (Figure 1E, Appendix 1—figure 3D). We also confirmed that cells in the back showed greater directional persistence toward the migration direction (Appendix 1—figure 4A–C), as (Saragosti et al., 2011) reported previously. All these results suggested that the bacteria in the back run more effectively toward the direction of group migration than those in the front.
To quantify the spatial extent of the drift efficiency, we defined the expected drift velocity $V}_{D}\equiv \frac{\u27e8{l}_{R}\cdot \mathrm{cos}{\theta}_{R}\text{}\u27e9}{\u27e8{\tau}_{R}+{\tau}_{T}\u27e9$ , by the projection of the average run length along the migration direction over the average duration of runs and tumbles (Dufour et al., 2014). This quantity reflects the effective run speed of runandtumble events that start running on a given location relative to the group. The drift velocity was found to decrease from the back to the front, crossing the group velocity ${V}_{G}$ in the middle of the group (Figure 1F). This particular trend of ${V}_{D}\left(z\right)$ suggests a mean reverting behavior of bacteria in the group: the cells at the back drift faster than the group ($V}_{D}>{V}_{G$), enabling the cells to catch up with the group; at the same time, the cells in the front drift slower than the group ($V}_{D}<{V}_{G$), causing them to slow down and fall back (Figure 1G). Such mean reverting process also results in subdiffusion of individuals relative to the group (Appendix 1—figure 3F), for which the mean square displacement (MSD) is constrained over time (Appendix 1—figure 3G). The slope of the spatial extent of ${V}_{D}\left(z\right)$, $r\text{}=0.05\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}\mathrm{i}{\mathrm{n}}^{1}$ , quantifies the speed at which individuals revert to its center.
We noted that the spatial profile of the expected drift velocity ${V}_{D}\left(z\right)$ was different from the instantaneous velocity ${V}_{I}\left(z\right)$ . This is because the instantaneous velocity ${V}_{I}\left(z\right)$ defines the average speed of cells in a given time interval $dt$, which reflects the positional shift of a group of cells at a given position. However, the expected drift velocity ${V}_{D}\left(z\right)$ defines the average speed of runandtumble events that start tumbling at a given position. Since the run duration is explicitly modulated by the gradient of chemoattractant $g\left(z\right)$ and is dependent on the chemotactic ability $\chi $ , ${V}_{D}\left(z\right)=\chi g\left(z\right)$ represents the drift velocity driven by the external stimuli (Celani and Vergassola, 2010; de Gennes, 2004; Dufour et al., 2014; Si et al., 2012).
The decreasing drift velocity implies a pushed wave front of group migration
To understand how the spatially structured profile of the drift velocity ${V}_{D}$ impacts on the group migration, we first adopted a Langevintype model that describes bacterial motions as an active Brownian particle driven by the expected drift velocity ${V}_{D}$ and a random force: $dx={V}_{D}dt+\u03f5dW$ (Berg, 2004). In this model, the runandtumble random motions are considered as a Gaussian random force $\u03f5dW$, while the cell motions are imposed to a deterministic force ${V}_{D}$ (Rosen, 1973; Rosen, 1974). The strength of Gaussian noise can be estimated by the effective diffusion coefficient $\u03f5=\sqrt{2D}$ , while the drift velocity is determined by the product of the perceived gradient $g\left(z\right)$ and the chemotactic ability $\chi $, ${V}_{D}=\chi g\left(z\right)$ . Such a stochastic description of bacterial motions has been proven equivalent to the classic KellerSegel (KS) model that described the population dynamics of bacterial chemotactic group migration (Keller and Segel, 1971a, Rosen, 1973). In the moving coordinate, $z=x\text{}\text{}{V}_{G}t$ , this Langevintype model specifies that $dz={V}_{D}\left(z\right)dt{V}_{G}dt+\u03f5dW$. Thereby, the cell motions relative to the migrating group are modulated by two effective forces: one generated by ${V}_{D}\left(z\right)$, which pushes the cells to catch up with the wave; and another generated by ${V}_{G}$ , which drags the cells to fall behind the wave. These two ‘forces’ constrain the random motions of individuals in an effective potential well $U\left(z\right)$ (Figure 2A).
To estimate the spatial profile of the effective driven force ${V}_{D}\left(z\right)$, we analyzed the KS model with moving ansatz (supplementary text). Assuming the density profile $\rho \left(z\right)$ directly measured from experiments (Figure 1B), we can deduce the chemoattractant concentration profile $S\left(z\right)$ (Appendix 1—figure 4D), as well as the perceived gradient $g\left(z\right)$ (Figure 2B). Since the perceived gradient $g\left(z\right)$ is almost linear in the main part of the wave profile, we approximated it by $g\left(z\right)\approx {g}_{0}+{g}_{1}z$ (Figure 2B, dashed line), which allows us to transform the stochastic model to an OUtype equation:
By this simplified model, we obtained a clear picture how individual behaviors are regulated relative to the group: cells are imposed to a driving force linearly dependent on their relative positions in the group, $F\left(z\right)=\chi {g}_{1}z+\left(\chi {g}_{0}{V}_{G}\right)$ . This suggests that the random motions of bacteria are constrained in a parabolic moving potential well, $U\left(z\right)=\frac{1}{2}\chi \left{g}_{1}\right{z}^{2}+\left(\chi {g}_{0}{V}_{G}\right)z+{U}_{0}$ (Figure 2A), where ${U}_{0}$ is set by $U\left(\infty \right)=0$. The position that minimizes $U\left(z\right)$ is also the balanced position, ${z}_{0}=\frac{{g}_{0}}{{g}_{1}}\frac{{V}_{G}}{\chi \left{g}_{1}\right}$ , where the driving force is null $F\left({z}_{0}\right)=0$. Behind the balanced position, cells experience a pushed force to catch up the group. Cells would start to fall back once they exceed the balanced position, where the driving force becomes negative. Therefore, the decreasing profile of the driving force enable cells to perform mean reverting behaviors around the balanced position. In addition, the expected rate that cells tend to revert back to the balanced position is defined by the slope of the spatially dependent driving force, $r=\chi \left{g}_{1}\right$ (Figure 2B).
Given the knowledge of individual behaviors, we studied the spatial distribution of population on the group migration. The OU model (Equation 1) which describes the singlecell stochastic motions has an equivalent form, known as the FokkerPlanck equation (Equation S19) which describes the spatialtemporal dynamics of cell density distribution. This population model provides a traveling wave solution with the mean position around ${z}_{0}$ and standard deviation $\sigma =\frac{\u03f5}{\sqrt{2\chi \left{g}_{1}\right}}$ .
From the solution, we noted that the effective driving force, as well as the drift velocity, has a negative slope (${g}_{1}<0$). The decreasing profile of the drift velocity suggests that the leading front of the group migrates as a pushed wave front (Gandhi et al., 2016; van Saarloos, 2003). As a key feature of the pushed wave, the leading front of the group drops parabolically, which is much faster than that of diffusion front. This further leads to a tight density profile of group migration for a single phenotype population.
The pushed wave front results ordered pattern of phenotypes
To address whether this pushed wave front still holds in presence of phenotypic diversity, we further examined the above analysis with cells of diverse chemotactic abilities ${\chi}_{i}$ imposed to the same moving perceived gradient $g\left(z\right)$ (Figure 2B). Given a monotonically decreasing profile of perceived gradient $g\left(z\right)$ , the driving force that each phenotypic individual experiences exhibits the same spatial dependency with the slopes depend on the intrinsic phenotypic properties of each phenotype ${r}_{i}={\chi}_{i}\left{g}_{1}\right$ . This monotonic dependency means that the balanced positions ${z}_{0}$ of the diverse phenotypes are orderly arranged. By the stochastic Langevintype model with phenotypic diversity in chemotactic ability, we first confirmed that each phenotypic population migrates at a constant speed ${V}_{G}$ , following the moving gradient $g\left(z\right)$ (see supplementary text). The density profiles of cells with different ${\chi}_{i}$ follow the same shape but are spatially orderly aligned (Figure 2C). Under the same moving gradient $g\left(z\right)$, the driving force ${\chi}_{i}g\left(z\right)$ is phenotypedependent, so that the bottom position of the potential well, ${z}_{0,i}=\frac{{g}_{0}}{{g}_{1}}\frac{{V}_{G}}{{\chi}_{i}\left{g}_{1}\right}$, is also spatially arranged according to ${\chi}_{i}$ (Figure 2D). As predicted by the OU model, the balanced positions ${z}_{0,i}$ of different phenotypes increase with their chemotactic abilities ${\chi}_{i}$ (Figure 2E, black line), while the standard deviation (${\sigma}_{i}$) of the density profiles decreases with ${\chi}_{i}$ (Figure 2E, blue line). We also confirmed the ordered structure of phenotypes by a particlebased model of the Langevintype dynamics coupled with chemoattractant consumption (Appendix 1—figure 5). Therefore, under the spatially decreasing driving force, cells with phenotypic diversity perform the same type of mean reverting processes with spatially ordered mean positions.
The spatial order of phenotypes does not directly promise a compact group migration with phenotypic diversity. By close examination of the density profile, we found that each phenotypic subpopulation propagates as a pushed wave front. We further calculated the total density profile of the entire migratory group with Gaussian distribution of chemotactic abilities under a decreasing linear gradient. Simulation reveals a pushed wave front for the combination of these subpopulations with different chemotactic abilities (Figure 2F, insert). In addition, we checked that the width of the entire group maintains in a converged width over long time, suggesting that the pushed wave profile enables the migratory population with diversity to keep in a tight shape (Figure 2F).
We examined the migration profiles under other forms of perceived gradient $g\left(z\right)$ : a constant perceived gradient and a spatially increasing perceived gradient. In the first case, individual bacteria of identical phenotype follow the diffusion process relative to the group, where the standard deviation of the population increasing with time by $\sigma =2\sqrt{Dt}$ (Figure 2G, dashed line). Each phenotype subpopulation is expected to have a constant drift velocity over space. However, as the drift velocity would vary by the chemotactic ability, each subpopulation migrates in different group speeds, suggesting a compact group migration of diverse population cannot be maintained in a constant perceived gradient (Figure 2G). In the latter case, when imposed to a spatially increasing driving force (a pulled wave, by definition), individual cells display superdiffusion processes. The density profile easily diverges in this case (Figure 2H). Therefore, we concluded that the pushed wave front, driven by a decreasing shape of driving force, enables a spatially ordered and compact pattern of phenotypes while on the move.
Spatially ordered individual behaviors predicted by agentbased simulation
Although the simplified OUtype model (Equation1) represents a key aspect of the ordering mechanism of phenotypes, it does not detail the signaling processes of bacterial chemotaxis, such as receptor amplification, adaptation, and motor responses (Sumpter, 2010), and it cannot predict bacterial runandtumble behavior. To consolidate the proposed mechanism underlying the emergence of spatial orders from the individual random motions, we further performed agentbased simulations integrated with the chemotactic pathway, multiflagella competition, and boundary effect in three dimensions (3D) (Dufour et al., 2014; Jiang et al., 2010; Sneddon et al., 2012). In the agentbased model, the attractant dynamics governed by diffusion and bacterial consumption is described by a reactiondiffusion equation (Equation S2) as previous multiscale models (Erban and Othmer, 2005; Xue and Othmer, 2009). We constructed a population with multiple phenotypes defined by different chemotactic abilities ${\chi}_{i}$ , where ${\chi}_{i}$ was varied by changing the receptor gain $N$ (for details, see supplementary text). Since the receptor gain $N$ only affects the amplification factor by which a cell responds to the gradient, the variation in bacterial motility $\u03f5$ is unchanged. As a result, a dense band of migrating cells that follow a selfgenerated moving chemoattractant gradient via consumption were recaptured as experiments (Appendix 1—figure 6A). The phenotypes were spatially ordered as $\chi $ varies (Appendix 1—figure 6B), and the velocity profile of each phenotype decreases from the back to front (Appendix 1—figure 6C) as predicted by the OUtype model.
To better analyze the simulations, we simplified the simulation with a nonconsumable attractant profile $S\left(z\right)$ moving with constant speed ${V}_{G}$ (Appendix 1—figure 4D). Using this simplified model, we first checked that the mean positions of the density profiles of cells with different receptor gains $N$, as well as their peaks, were orderly aligned with respect to chemotactic ability ${\chi}_{i}$ .
As an important advantage of the agentbased simulations, the model allowed us to analyze singlecell behavior during the ordered group migrations. For each phenotype $i$, the expected drift velocity ${V}_{D,i}\left(z\right)$ decreased along the density profile (Figure 3A). Consistent with the ordered structure of the density profiles, the intersection between ${V}_{D,i}\left(z\right)$ and ${V}_{G}$ exhibited the same order of chemotactic ability ${\chi}_{i}$ (Appendix 1—figure 7). As the reversion rate ${r}_{i}=\left\frac{d{V}_{D,i}\left(z\right)}{dz}\right$ showed a positive correlation to ${\chi}_{i}$ , cells with lower receptor gain $N$ (resulting in a smaller $\chi $) experienced a weaker reverting force toward the center (Figure 3A insert). Thus, the effective moving potential, ${U}_{i}\left(z\right)$ , which constrains the cells around the mean positions sorted by their chemotactic abilities, becomes flat for cells with lower chemotaxis ability $\chi $ (Figure 3B; Long, 2019). As a result, cells of each phenotype perform subdiffusion, whereby the MSD along the migration coordinate relative to the group was bound at a level negatively correlated to $\chi $ (Figure 3C). The width of the density distribution, as an effect of the reversion force, decreased with the reversion rate, as an approximate linear function. Using this agentbased model, we further obtained similar results for populations of different ${\chi}_{i}$ through adaptation time $\tau $, or basal CheY protein level ${Y}_{p0}$ , which determined the basic tumble bias $T{B}_{0}$ (Dufour et al., 2014; Jiang et al., 2010; Sneddon et al., 2012; Appendix 1—figure 8).
Experiments with titrated phenotypes confirm ordered behavior modulation
To verify the model predictions on the individual behaviors of different phenotypes, we experimentally measured the trajectories of cells with different chemotactic abilities during the group migration. Specifically, we altered the chemotaxis abilities of cells by titrating the expression level of Tar, which is under the control of a small molecule inducer aTc (Sourjik and Berg, 2004; Zheng et al., 2016) (see Materials and methods and Figure 4A). The variations in the expression of Tar would lead to different receptor gains in response to the Asp gradient (Adler, 1966b; Adler, 1969; Sumpter, 2010), but the tumble bias and growth rate would not change. Using the migration speed of the bacterial range expansion to quantify the chemotaxis ability of the titrated strains, we found that the chemotaxis ability increases with aTc concentration (see Materials and methods and Appendix 1—figure 9).
The Tartitrated cells labeled with yellow fluorescent protein (YFP; strain JCY20) were added to the wildtype population at a ratio of 1 in 400. Within the wildtype population, 1 in 50 cells was labeled with red fluorescent protein (RFP) (strain JCY2). As the Tartitrated strain constituted a small portion of the premixed population, we considered that the density profile of the population would be invariant to different levels of induction of Tar. The premixed population could generate collective group migration, similar to the wildtype population (Appendix 1—figure 9). The trajectories of the YFPlabeled cells were tracked to represent the behavior of cells with different chemotactic abilities, while the profile of wildtype cells with RFP was also measured to characterize the density distribution of the entire migratory population.
By comparing the statistics of cells with different Tar expression levels, we found that the expected drift velocity ${V}_{D,i}\left(z\right)$ followed the same decreasing pattern from back to front (Figure 4B). More importantly, as the Tar expression level (chemotactic ability) increased, the slope of the decreasing pattern increased, which was consistent with the model prediction shown in Figure 3A. The intersections between ${V}_{D,i}\left(z\right)$ and ${V}_{G}$ , as well as the peak positions and mean positions of each Tartitrated density profile (Figure 4C), shifted toward the front as the chemotactic ability increased (as measured by the migration rate on agar plates; Cremer et al., 2019; Liu et al., 2019). The ${V}_{D}$ cross point was always behind the peak position and the mean position (Figure 4C), suggesting that cells were leaking behind. Moreover, the width of each Tartitrated density profile (defined by $2{\sigma}_{i}$) decreased as the reversion rate ${r}_{i}$ increased (Figure 4D), consistent with the model results in Figure 3C. Thus, as the OUtype model predicts, the width of the density profile is controlled by the reversion rate determined by the chemotactic ability ${\chi}_{i}$ .
Discussion
In summary, coordinated behaviors with ordered spatial arrangements of phenotypes are abundant in a wide range of biological and humanengineered systems, and are believed to involve elaborate control mechanisms. For animal migrations, it is challenging to characterize simultaneously the computational strategy and behavior at the individual level so as to avoid averaging out phenotypic diversity, and the emergent behavior at the population level (Couzin et al., 2005; Couzin et al., 2002; Vicsek and Zafeiris, 2012). For bacterial chemotactic migration, cells with different phenotypes are spatially aligned based on their chemotactic abilities. This observation was explained as a selfconsistent result with the decreasing profile of attractant predicted by KS model (Fu et al., 2018). In this study, we demonstrated that the collective consumption of attractant by bacterial group generates a spatial structure of individual drift velocity along the migrating group profile. Such a spatial profile of drift velocity results a pushed wave front on population level, and maintains diverse phenotypes in a compact migration group. Moreover, this pushed wave front enables spatial modulation of individuals to perform mean reverting random motions around centers sequentially aligned by their chemotactic abilities, thereby giving rise to a spatially ordered pattern. Therefore, we demonstrated that the population order could emerge among diverse individuals that following the same rule of behavioral modulation.
This strategy of selforganization does not require sophisticated communications (Curatolo et al., 2020; Karig et al., 2018; Liu et al., 2011; Payne et al., 2013) nor other hydrodynamic interactions (Chen et al., 2017; Drescher et al., 2011; Zhang et al., 2010) among individuals. Our observation of the decreasing drift velocity can imply the effective perceived gradient that cells experience. We believe that this spatial profile is mainly contributed by the consumption of the chemoattractant. By using a Taronly strain (UU1624) (Gosink et al., 2006), we demonstrated that the mutant could also generate a stable group migration in our experimental condition similar to the wildtype strain (Appendix 1—figure 9G). This further suggests that the secretion of selfattractants is unlikely a necessary condition of collective group migration (Cremer et al., 2019; Fu et al., 2018; Keller and Segel, 1971a), although there are doubts about the existence of selfattractants in high density (Budrene and Berg, 1995).
The spatially dependent drift velocity provides a structured driving force of a migration group, resulting a pushed wave front. Pushed and pulled waves are determined by the spatial distribution of the spreading velocity of a propagation front (van Saarloos, 2003; Figure 2F–H). The properties of pushed and pulled waves have been discussed in growthdriven range expansion (Birzu et al., 2018; Erm and Phillips, 2020; Gandhi et al., 2016), where the wave type is determined by densitydependent growth rates. A prominent example of pulled wave is known as the Fisher wave (Fisher, 1937), where the population expansion is driven by constant diffusion and logistic growth of individuals. However, in such biological systems, the spatial dependence of front speeds is hardly quantified in experiments. Here, by direct measurement of drift velocity on singlecell level, we identified the chemotactic migration group of bacteria as a pushed wave. This chemotaxis system would provide a unique multiscale model to study further details of pushed wave.
The spatially decreasing profile of driving force does not only cause an ordered pattern of phenotypes, but also results a pushed wave front that enables a negative feedback control on the propagating speed. This further allows a compact density profile for heterogenous population to migrate at the collective level. The advantageous to keep diversity in the pushed wave front was also reported in the growthdriven range expansion system (Birzu et al., 2018). Our study revealed that, other than spatial regulation of fitness, the direct modulation of individual drift velocity in space could also maintain diversity in range expansion.
Detailed analysis of the spatial structured driving force could also provide the limits of phenotype that is allowed in the group. In the migratory group, the same rule of behavioral modulation applied to cells with different phenotypes, such that the random motions of cells were bound by moving potential wells whose basins were sequentially aligned. However, it is noteworthy that cells could skip the potential wells from the back because the ‘driving force’ decreased again at the far back of the group (Long, 2019). This results in leakage of cells in the migratory group (Holz and Chen, 1978; Novickcohen and Segel, 1984; Scribner et al., 1974). Phenotypes with weaker chemotactic abilities were located at the back of the group, where the effective potential well was shallower (Figure 3C), allowing for more chances to skip. Thus, such collective migration selects bacteria with higher chemotactic abilities (Liu et al., 2021; Liu et al., 2019).
The simple computational principle of behavioral modulation to allocate different phenotypes in the collective group is likely not limited to sensing the selfgenerated signal by consumption of attractant. A prominent example of trailfollowing migration (Couzin and Krause, 2003; Helbing et al., 1997) and a typical class of collective behavior is represented by a modified Langevintype model, where individuals tracing the accumulated signal secreted by all participants (Equation S20) can reproduce similar spatialtemporal dynamics of behavioral modulation, as well as ordered arrangements of phenotypes in the migratory group (Appendix 1—figure 10). Thus, this mechanism of matching individual abilities by the signal strength might provide an explanation of how other higher organisms organize ordered structures during group migration.
Materials and methods
Strains
The wildtype strain E. coli (RP437) and its mutants were used in this study, where all plasmids were kindly provided by Dr Chenli Liu. Specifically, the Tartitratable strain was constructed by recombineering according to previous research (Zheng et al., 2016). Specifically, the DNA cassette of the PtettetRtar feedback loop was amplified and inserted into the chromosomal attB site by recombineering with the aid of plasmid pSim5. The tar gene at the native locus was seamlessly replaced with the aph gene by using the same recombineering protocol. To colorcode the strains, we use plasmids with chloramphenicolresistant gene carrying yfp under constitutive promoter (for JCY1 strain) and pLambdadriven mRFP1 plasmids maintained by kanamycin (for JCY2). To colorcode Tartitratable strain (JCY20), a plasmid carrying yfp chloramphenicolresistant gene was transformed into constructed Tartitratable strain. The taronly strain (UU1624) was modified from RP734 and was kindly provided by Prof. Johan Sandy Parkinson.
Media and growth conditions
Request a detailed protocolFor bacterial culture, the M9 supplemented medium was used. The preparation of the M9 supplemented medium follows the recipe in previous study (Fu et al., 2018): 1×M9 salts, supplemented with 0.4% (v/v) glycerol, 0.1% (w/v) casamino acids, 1.0 mM magnesium sulfate, and 0.05% (w/v) polyvinylpyrrolidone40. 1×M9 salts were prepared to be 5×M9 salts stock solution: 33.9 g/L Na_{2}HPO_{4}, 15 g/L KH_{2}PO_{4}, 2.5 g/L NaCl, 5.0 g/L NH_{4}Cl.
For migration experiments in the microchannel, the M9 motility buffer was used. The recipe was:1×M9 salts, supplemented with 0.4% (v/v) glycerol, 1.0 mM magnesium sulfate, and 0.05% (w/v) polyvinylpyrrolidone40, 0.1 mM EDTA, 0.01 mM methionine, and supplemented with 200 µM aspartic acid.
For the migration rate measurements, the M9 amino acid medium with 0.2% (w/v) agar was used to prepare swim plate (Liu et al., 2019). The recipe was: 1× M9 salts, supplemented with 0.4% (v/v) glycerol, 1× amino acid, 200 µM aspartic acid, 1.0 mM magnesium sulfate, and 0.05% (w/v) polyvinylpyrrolidone40. 1× amino acid were prepared to be 5× amino acid stock solution: 4 mM alanine, 26 mM arginine (HCl), 0.5 mM cysteine (HCl·H_{2}O), 3.3 mM glutamic acid (K salt), 3 mM glutamine, 4 mM glycine, 1 mM histidine (HCl·H_{2}O), 2 mM isoleucine, 4 mM leucine, lysine, 1 mM methionine, 2 mM phenylalanine, 2 mM proline, threonine, 0.5 mM tryptophane, 1 mM tyrosine, 3 mM valine. All experiments were carried out at 30°C. Plasmids were maintained by 50 µg/mL kanamycin or 25 µg/mL chloramphenicol.
Sample preparation
Request a detailed protocolThe bacteria from frozen stock were streaked onto the standard LuriaBertani agar plate with 2% (w/v) agar and cultured at 37°C overnight. Three to five separate colonies were picked and inoculated in 2 mL M9 supplemented medium for overnight culture with corresponding antibiotics to maintain plasmids. The overnight culture was diluted by 1:100 into 2 mL M9 supplemented medium the next morning. For Tar titration strains, related aTc were added in this step. When the culture OD600 reaches 0.2–0.25, it was then diluted into prewarmed 15 mL M9 supplemented medium so that the final OD600 was about 0.05 (Liu et al., 2019; Zheng et al., 2020; Zheng et al., 2016).
Bacteria were washed with the M9 motility buffer and were resuspended in fresh M9 motility buffer to concentrate cell density at OD600 about 1.0. Then, the wildtype strain and fluorescent strain were mixed with ratio of 400:1 before loaded in the microfluidic chamber (Fu et al., 2018; Saragosti et al., 2011). For Tar titration experiments, the wildtype strain (RP437) was mixed with two fluorescent strains (JCY2 and JCY20) by 400:8:1.
Microfabrication
Request a detailed protocolThe microfluidic devices were fabricated with the same protocol and the same design as previous research (Bai et al., 2018; Fu et al., 2018), except that the capillary channel was designed longer than that of previous ones. The size of the main channel was 20 mm×0.6 mm× 0.02 mm and only one gate at the end of the channel was kept (Appendix 1—figure 1A).
Band formation
Request a detailed protocolFluorescent cells were mixed with nonfluorescent cells by 1:400 for cell tracking in the dense band. Sample of mixed cells with density OD600 ≈1.0 was gently loaded into the microfluidic device and then the device was spun for 15 min at 3000 rpm in a 30°C environmental room so that almost 100,000150,000 cells were placed to the end of the channel. After spinning, the microfluidic device was placed on an inverted microscope (Nikon TiE) equipped with a custom environmental chamber set to 50% humidity and 30°C.
Imaging
Request a detailed protocolThe microscope and its automated stage were controlled by a custom MATLAB script via the μManager interface (Edelstein et al., 2014; Fu et al., 2018). About 30 min after the sample loading, a 4× objective (Nikon CFI Plan Fluor DL4X F, NA 0.13, WD 16.4 mm , PhL) was placed in the wave front for imaging. The fluorescent bacteria, seen as randomly picked samples of the migrating group, were captured continuously with frame rate of 9 fps in 10 min until they leave the view. Typical tracks are longer than 300 s. Timelapsed images were acquired by a ZYLA 4.2MP Plus CL10 camera (2048 × 2048 arrays of 6.5 µm×6.5 µm pixels) at 9 frames/s (fps) . An LED illuminator (0034RXCite 110LED) and an EYFP block (Chroma 49003; Ex: ET500/X 20, Em: ET535/30 m) compose the lightening system.
For the Tar titration experiments, the channel was first scanned with 10× objective (CFI Plan Fluor DL 10× A, NA 0.30, WD 15.2 mm, PH1) enlighten by an LED illuminator (0034RXCite 110LED) through the RFP block (Chroma 49005, Ex: ET545/X 30, Em: ET620/60 m) and EYFP block channels for seven neighbored views around the migration group. These images were further combined to two large pictures of the RFP strains and YFP strains. The channel was scanned twice, respectively before and after the 10 min tracking of fluorescent Tartitrated cells.
Track extraction and state assignment
Request a detailed protocolThe acquired movie was first analyzed with the Utrack software package to identify bacteria and to get their trajectories (Jaqaman et al., 2008). Then the tracks were labeled by run state and tumble state by a custom MATLAB package (Waite et al., 2016) using a previously described clustering algorithm (Dufour et al., 2016).
Track analysis
Request a detailed protocolThe group velocity ${V}_{G}$ was calculated by averaging the frametoframe velocity ($dt\approx 0.11\mathrm{s}$) over all tracks and all time. The cell number for the first frame over a spatial bin of $\mathrm{\Delta}x=60\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$ and a channel section a=12,000 µm^{2} were calculated to get the density profile $\rho \left(x,t=0\right)=\frac{\sum i\left(x,t\right)}{a\cdot dx}$ . The peak position of the first frame (${x}_{peak}\left(t=0\right)$) was then determined by the maximum of $\rho (x,t=0)$. The position of each bacterium (${x}_{i}\left(t\right)$) was transformed to moving coordinate position ${z}_{i}$ by the group velocity ${V}_{G}$ and origin of the axis on the density peak by ${z}_{i}={x}_{i}\left(t\right){V}_{G}t{x}_{peak}\left(t=0\right)$ . Given the relative position of each cell, we recalculated the density profile in moving coordinate $\rho \left(z\right)=\frac{\sum i\left(z\right)}{a\cdot dx}$ . The width of the density profile was defined by two times the standard deviation of relative positions $2\sigma =2\sqrt{\frac{1}{n1}{\sum}_{i=1}^{n}{\left({z}_{i}\langle z\rangle \right)}^{2}}$ . The spatial distribution of the instantaneous velocities $\langle {V}_{I}\left(z\right)\rangle $ was calculated by averaging the velocity in spatial bin of $\mathrm{\Delta}z=240\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$.
A tumblerun event is the minimal element of bacterial behavior. The typical spatial scale of a tumblerun event is about 20 μm, which is much smaller than the spatial bin size chosen in this study (240 µm). The spatial distributions of run time $\langle {\tau}_{R}\left(z\right)\rangle $, tumble time $\langle {\tau}_{T}\left(z\right)\rangle $, and run length $\langle {l}_{R}\left(z\right)\rangle $ were calculated by averaging the related values of all the events with tumbling position (${z}_{T}$) located in each spatial bin ($z$). As the displacement of tumble is small, the tumbling position is approximately the starting position of runs. For each tumblerun event, we have the vector linking starting position and end position of the run. The running angle ${\theta}_{R}$ is then defined by the angle between run direction and the group migration direction. One can easily deduce all the other quantities with the formulations in Table 1.
Growth rate and migration rate measurement
Request a detailed protocolGrowth rates of Tartitrated strains were calculated from exponential fitting (${R}^{2}>0.99$) over measured curves of cell density (OD600) with respect to time. A 250 mL flask with 20 mL M9 supplement medium were used. All measurements were performed in a vibrator of rotation rate of 150 rpm at 30 °C. OD600 was measured by a spectrophotometer reader every 25 min. Each strain has been measured for at least three times.
The semisolid agar plate was illuminated from bottom by a circular white LED array with a light box as described previously (Liu et al., 2011; Liu et al., 2019; Wolfe and Berg, 1989) and was imaged at each 2 hr by a camera located on the top. As bacteria swimming in the plate forms ‘Adler ring’, we used the first maximal cell density from the edge to define the moving edge of bacterial chemotaxis. The migration rate was then calculated from a linear fit over the data of edge positions in respect to time (${R}^{2}>0.99$).
Models and simulations
Request a detailed protocolDetails of the theoretical models and numerical simulations were presented in the appendix notes. In which, the Langevin equation was deduced and solved numerically with a particlebased simulation; the approximated OUtype equation and its traveling wave solution were deduced; an agentbased simulation of bacterial with chemotaxis pathway was performed following previous works (Dufour et al., 2014; Jiang et al., 2010; Sneddon et al., 2012).
Appendix 1
Theory and numerical simulations
Summary
In this note, we first deduced the Langevintype equation from the classic KS model that describes bacteria group migration in a capillary. The Langevintype equation describes the dynamics of singlecell particle motion as active particles (Section 1). Together with the chemoattractant consumption equation, we investigated the Langevintype equation through particlebased simulations for a group of particles of single phenotype or multiphenotypes (Section 2). Then, we transformed the equations into moving coordinate, so that the stable density profiles as well as the chemoattractant profiles were observed (Section 3). In this moving coordinate, the chemoattractant consumption can be decoupled from particle dynamics, which allows us to simulate particle motions over a nonconsumable field of chemoattractant on the move (Section 4). Furthermore, the Langevin dynamics describes that the particle motion was simplified to an OUtype process, of which the particle dynamics can be solved analytically (Section 5). To gain more insights of the chemotactic response, we described an agentbased model which integrated the bacterial chemotaxis pathway in Section 6. And specifically, we deduced the chemotactic ability $\chi $ from the pathway parameters (Section 7). The agentbased simulations were performed under a moving gradient, from which ordered structures of phenotypes were recaptured (Section 8). The model was further extended to system with signal secretion and shows similar results (Section 9).
1. The Langevin equation with chemoattractant consumption
The bacterial chemotactic migration in a capillary was first modeled by Keller and Segel in 1971 (Keller and Segel, 1971a). In this model, the bacterial cell division was neglected as the nutrition in the capillary was poor. In 1D, the dynamics of cell density $\rho \left(x,t\right)$ is a combination of diffusion and chemotaxis:
where $D$ and $\chi $ respectively represent the diffusion coefficient and chemotactic ability of bacteria, and $g\left(x,t\right)$ represents the perceived gradient of chemoattractant.
The dynamics of the chemoattractant $S(x,t)$ follows a reactiondiffusiontype equation with a diffusion coefficient ${D}_{s}$ and consumption with a constant rate $k$:
From the chemoattractant profile $S(x,t)$, we can get $g\left(x,t\right)$ by:
with ${K}_{on},{K}_{off}$ were dissociation constants for active and inactive conformation of the receptor, respectively.
Following Rosen, 1973, one can deduce the FokkerPlanck equation associated to the KS model that describes the probability $P=P(x,y;t)$ to observe of a cell which initiates from position $y$ at position $x$ and time $t$. The related equation writes:
Equation S4 is equivalent to the Langevin dynamics as:
where ${x}_{i}$ represents the position of the ith bacteria and $g\left({x}_{i},t\right)$ is the perceived gradient of each bacteria at time $t$. $dW$ is a Wiener process that describes a random walk with noise level $\u03f5=\sqrt{2D}$ .
As previously proofed, the expected drift velocity equals the production of chemotactic ability $\chi $ and the perceived gradient $g$, ${V}_{D}=\chi g$ (Waite et al., 2016). So that, the Langevintype Equation S5 can be reformed to:
The initial condition of the chemoattractant was $S\left(x,0\right)={S}_{0}=200\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{M}$. In the stochastic model, the density profile $\rho \left(x,t\right)$ was calculated by the sum of cell number ($\sum i(x,t)$) in certain location $[x,x+\Delta x]$ and time $t$ over volume $a\cdot \Delta x$, where $a$ is the capillary section area and $dx$ is the spatial bin size:
2. Particlebased simulation with consumption
The Langevintype model in the absolute coordinate $x$ (Equation S5) was simulated with a customized code using MATLAB. In each case, 100,000 particles were simulated in a 1D space where $x\in \left[0,\text{}20\right]\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}\mathrm{m}$. The absorption boundaries were applied for $x=0\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}\mathrm{m}$ and $x=20\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}\mathrm{m}$, so that particles stay at the boundary until they leave them actively ($x\left(x>0\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}\mathrm{m}\right)=0$ and $x\left(x>20\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}\mathrm{m}\right)=20\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}\mathrm{m}$). Initially all particles were placed at $x=0\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}\mathrm{m}$ and their positions were calculated according to Equation S7 at each time step of $\Delta t=0.05s$ for 60 min. The random noise follows a Gaussian distribution generated by MATLAB function ‘randn.m’ with amplitude multiplied by $\sqrt{2}\u03f5\sqrt{\Delta t}$ .
Unlike the particles that change their position continuously, the spatial profile of cell density $\rho \left(x,t\right)$ was first calculated on a space mesh of grid $\mathrm{\Delta}x=10\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$, and then linearly interpolated to a finer mesh of grid $\mathrm{\Delta}x=0.1\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$. The chemoattractant $S\left(x,t\right)$ and perceived gradient $g\left(x,t\right)$ were deduced accordingly. At each time point, following Equation S7, cells located on grid $\left[x,x+\Delta x\right]$ were summed and were divided over $a\Delta x$ to get the density profile. In order to avoid fast consumption near $S=0$, a saturation effect of consumption rate was added to Equation S3 using $S/(S+{K}_{s})$.
where $K}_{s}=0.1\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{M$, of which the value was small and does not affect the linear relation in most range of $S>{K}_{s}$ .
Then the chemoattractant profile $S\left(x,t\right)$ was calculated from $\rho \left(x,t\right)$ according to Equation S8, using a fourthorder RungeKutta method from an initial condition of $S\left(x,0\right)={S}_{0}$ . The perceived gradient $g\left(x,t\right)$ was further calculated from $S\left(x,t\right)$ on Equation S3. Having the discrete $g\left(x,t\right)$ over grid of $\mathrm{\Delta}x=0.1\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$, the perceived gradient of each particle $g\left({x}_{i}\right)$ was assigned by detecting the particle position over the spatial mesh grid.
Following the above instructions, density profile and peak positions with bacteria of $\chi =0.11\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}{\mathrm{m}}^{2}\mathrm{m}\mathrm{i}{\mathrm{n}}^{1}$ and $D=\chi /22$ were plotted in Appendix 1—figure 5A. The group migration reached steady state after 20 min simulation. The peak velocity was then defined by linear fitting of peak positions for $t\in \left[20,\text{}60\right]\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}\mathrm{i}\mathrm{n}$, which was found to be $V}_{p}=3.03\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}/\mathrm{s$.
To simulate a bacterial group of multiple phenotypes, we constructed 100,000 cells with their chemotactic ability $\chi $ follows a Gaussian distribution with the average chemotactic ability $\u27e8\chi \u27e9=0.11\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}{\mathrm{m}}^{2}\mathrm{m}\mathrm{i}{\mathrm{n}}^{1}$ and s.d. of 0.03 mm^{2}min^{1} (Appendix 1—figure 5B) while the diffusion coefficient $D$ was fixed to $\langle \chi \rangle /22$ for all cells. The consumption rate of each bacterium was also set constant of $4.6\cdot {10}^{9}\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{M}/\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}/\mathrm{s}$. The group of multiphenotypes migrates slower than single phenotype with peak velocity $V}_{p}=2.90\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}/\mathrm{s$ (Appendix 1—figure 5A). The bacterial positions after 20 min simulation were further transformed to the moving coordinate ${z}_{i}={x}_{i}{V}_{p}t$. They were subgrouped according to their chemotactic ability $\chi $ to six subgroups. And the corresponding average density profiles of each subgroups were plotted in Appendix 1—figure 5D with colors defined in Appendix 1—figure 5C. The density profile of subgroups was found ordered according to their average chemotactic ability $\chi $. The timeshifted profiles, $S\left(z\right)$ and $g\left(z\right)$, were shown in Appendix 1—figure 5E,F.
The simulations of bacteria with consumption confirmed that coherent group migration and ordered spatial structures can emerge from a system of active particles modulated by gradient of consumable chemoattractant.
3. Traveling wave solution
As the density profile was found stable during the experimental time scale (30 min), we are allowed to use a moving coordinate $z=x{V}_{G}t$ in searching for a traveling wave solution. Equation S2, Equation S5, Equation S7 become:
Given a stable density profile measured experimentally, we can integrate Equation S10 to get a stable profile of chemoattractant $S\left(z\right)$. As the chemoattractant was unconsumed at $z=\infty $ and it is completely consumed as the wave passes $z=\infty $, the boundary condition of $S\left(z\right)$ can be written as $S\left(\infty \right)={S}_{0}$ and $S\left(\infty \right)=0$. We get:
where $M\left(z\right)={\int}_{z}^{\infty}\rho \left(z\right)dz$ is the accumulated number of cells up to position $z$. We then used the experimentally measured cell density profile together with the group speed ${V}_{G}$ to calculate the moving attractant profile $S\left(z\right)$ (Appendix 1—figure 5A). The perceived gradient $g\left(z\right)$ in the moving coordinate was then calculated by:
In a stable migration group, the $g\left(z\right)$ deceases from back to front over the density profile. Suppose the gradient does not steepen toward the back, then cells with lower chemotactic abilities will fall behind and see shallower gradient, drifting even slower and causing the wave to spread out. Since the gradient is formed by local consumption, the front of the wave will have less cells and become shallower, while the back of the wave will accumulate more cells and become steeper. This will automatically readjust the gradient steepness until the back becomes steep enough to allow cells will lower chemotactic abilities to catch up, stabilizing the gradient and the traveling wave.
4. Particlebased model with preset moving gradient
The chemoattractant profile $S\left(z\right)$ can be deduced by an experimentally measured stable density profile $\rho \left(z\right)$ in the moving coordinate. With the experimental mesh grid $\mathrm{\Delta}x=240\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$, we plot the average $S\left(z\right)$ over three biological replicates in Appendix 1—figure 4D (circles). In order to apply this profile into the particlebased simulation, we first expanded the space range of $S\left(z\right)$ to $x\in \left[4800,\text{}4800\right]\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$ by setting $S\left(z<2400\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}\right)=0$ and $S\left(z>2400\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}\right)={S}_{0}$ . Next, we used a cubic spline data interpolation to construct a profile of $S\left(z\right)$ with much finer spatial mesh grid $\mathrm{\Delta}x=0.001\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$ (Appendix 1—figure 4D, line). Therefore, a $g\left(z\right)$ profile on finer meshes was calculated accordingly (Figure 2A). The $g\left(z\right)$ profile was then translated to an absolute coordinate by $g\left(x,t\right)=g\left(z\right)+{V}_{G}t$, with $V}_{G}=3.0\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}/\mathrm{s$.
Using the preset moving gradient $g\left(x,t\right)$ , we can simulate the particle motions without calculation of $\rho (x,t)$ and $S\left(x,t\right)$ . We used the same method as described in Section 2 to simulate particle motions, and at each time point $t$, we assigned $g\left({x}_{i}\right)$ by indexing bacterial position on $g\left(x,t\right)$ , with spatial precision of $\mathrm{\Delta}x=0.001\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$.
For bacteria of four different $\chi $ ranging from 0.15 to 0.3 mm^{2}min^{1} , we simulated 10,000 particles for each phenotype with the same noise level $\u03f5=15\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}\mathrm{m}\cdot \mathrm{m}\mathrm{i}{\mathrm{n}}^{0.5}$ . During simulation of 60 min, bacteria of all phenotypes move as a group with peak positions linearly increasing with time. The fitted slope of all peaks equals the preset velocity of the moving gradient ${V}_{G}$ . In the moving coordinate $z$, the density profiles of all phenotypes were ordered (Figure 2A), and they were clearly ordered by $\chi $. The peak positions and mean positions of the density profile increase with $\chi $ (Figure 2C), while the width (defined by $\sigma $) decreases with $\chi $ (Figure 2D).
5. The type model
The perceived gradient $g\left(z\right)$ deduced from experimental data shows linear decreasing trend in the major range of density profile $0.4\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}\mathrm{m}<\mathrm{z}<0.4\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}\mathrm{m}$. It can be then linearly fitted by $g\left(z\right)\approx {g}_{1}z+{g}_{0}$ with $g}_{1}\approx 2.1\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}{\mathrm{m}}^{2$ and $g}_{0}\approx 1.2\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}{\mathrm{m}}^{1$ . Subscribe this linear approximation in Equation S7, we get:
Now, we compare Equation S10 to the standard OU model given by $dx=\lambda xdt+\u03f5dW$. Equation S9 can be taken as an OUtype process with $\chi {g}_{1}$ as the reversion rate $\lambda $ and the mean position deviated from origin to ${z}_{0}$:
where the mean position ${z}_{0}$ is the steadystate position of its ensemble average $\langle z\rangle $:
The original OU model describes an active particle motion under a linear force field around $x=0$. Since Equation S10 can be taken as a generalized OU model, we can thus consider our bacteria in the migration group as an active particle under an effective force field of $\chi {g}_{1}z+\left(\chi {g}_{0}{V}_{G}\right)$ . Integrating this force field over $z$, we get an effective potential well $U\left(z\right)$ :
The ${U}_{0}$ is an integral constant controlled by boundary condition, which was set by $U\left(\infty \right)=0$ in this article. The effective potential well $U\left(z\right)$ generated by $U\left(z\right)=\int \left(\chi g\left(z\right){V}_{G}\right)dz$ from the particlebased simulation and Equation S17 are compared in Figure 2B.
Using a variable substitution of $\stackrel{~}{x}=z(\chi {g}_{0}{V}_{G})/{g}_{1}$ , and applying the general solution of the standard OU model of $\stackrel{~}{x}\left(t\right)={\stackrel{~}{x}}_{0}{e}^{\lambda t}+{\int}_{0}^{t}\u03f5dW\cdot {e}^{\lambda \left(t{t}^{`}\right)}d{t}^{`}$ (Cherstvy et al., 2018), we get the mean square displacement: $MSD\equiv \langle {\left(z\left(t\right){z}_{0}\right)}^{2}\rangle ={\langle {z}_{0}\rangle}^{2}{\left(1{e}^{\chi {g}_{1}t}\right)}^{2}+\frac{{\u03f5}^{2}}{2\chi {g}_{1}}\left(1{e}^{2\chi {g}_{1}t}\right)$ . The standard distribution $\sigma $ of the density distribution is the limit of MSD:
The trend of both mean positions (Figure 2C) and width (Figure 2D) was successfully predicted by analytical solutions of the OUtype model Equation S15 and Equation S18.
Moreover, the OUtype model Equation S14 has an equivalent FokkerPlanck equation that describes the probability density function of particles on $z$, $P(z,t)$:
At $t\to \infty $, the probability density function reaches the steady stat $P\left(z\right)$, where the density profile follows $0=\frac{d}{dz}\left(\left(\chi {g}_{1}{z}_{i}+\left(\chi {g}_{0}{V}_{G}\right)\right)P\left(z\right)\right)+D\frac{{d}^{2}P\left(z\right)}{d{z}^{2}}$ . Solving this equation with boundary conditions of $P\left(z=\pm \mathrm{\infty}\right)=0$ and $\frac{\mathrm{\partial}P\left(z=\pm \mathrm{\infty}\right)}{\mathrm{\partial}z}=0$ , we get the steady stat density profile $P\left(z\right)={e}^{\frac{1}{D}(\frac{1}{2}\chi {g}_{1}{z}^{2}+\left(\chi {g}_{0}{V}_{G}\right)z+{C}_{0})}$ with ${C}_{0}$ is the integral constant that allows ${\int}_{\infty}^{\infty}P\left(z\right)dz=1$. The exponent of this density profile is positiondependent that drops in parabolic form as $z\to \infty $. Compare to the front of Fisher wave, where the exponent of the density profile decreases linearly with $z$ ($P\left(z\right)={e}^{Cz}$), we know that the density profile of a chemotaxis migration wave is more compact than a Fisher wave.
6. Agentbased simulation with chemotaxis pathway
In order to simulate more details of bacterial group migration in a capillary, we use a previously established chemotactic pathway integrated agentbased simulation. This model integrated bacterial chemotactic signaling pathway to determine the switch rate of motors, the conformation changes of flagella, as well as the competition of multiple flagella.
We first simulated 80,000 cells evenly distributed in four phenotypes of different receptor gain $N$ with consumption of chemoattractant. The cells were simulated as described in a previous paper (Saragosti et al., 2011), while the chemoattractant field S was calculated at each time step ($\mathrm{\Delta}t=0.01\phantom{\rule{thinmathspace}{0ex}}\mathrm{s}$) according to Equation S8. In a 3D channel of 0.2 mm×0.01 mm× 10 mm, space was discretized to grid of 2.5 µm×2.5 µm× 2.5 µm . The evolution of the attractant during each time step was calculated by the sum of the consumption of every cell in the neighbor grids weighted by their effective diffusion distance. The consumption rate was fitted for reasonable migration speed. Detailed methods were described in a previous paper (Xue and Othmer, 2009). The perceived attractant of each cell was then linearly interpolated from the closest grids of the attractant field.
Simulation starts by setting all cells in one end of the channel ($x=0$), randomized in y and z directions and setting uniform chemoattractant concentration $S}_{0}=200\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{M$ within the channel. Flip boundary condition was used for all three directions. After 20 min simulation, the bacterial population from stable migration groups (Appendix 1—figure 6A) tracks with cellular runtumble stat was recorded from 30 to 32 min with time step of 0.1 s. Using the group velocity ${V}_{G}$ of this period as the speed of moving coordinate $z=x{V}_{G}t$, density distribution of subgroups of different receptor gain $N$ was plotted to show spatial ordering of subgroups (Appendix 1—figure 6B). The tracks were analyzed with exactly the same method as the experimental tracks (see Materials and methods). The expected drift velocity ${V}_{D}\left(z\right)$ of all subgroups decreases from back to front and was also spatially ordered (Appendix 1—figure 6C).
The particlebased simulations (Section 4) have already shown that the chemoattractant consumption can be decoupled from bacterial motions. To simplify the analysis, we used the same nonconsumable moving gradient $g\left(x,t\right)$ described in Figure 2A. Using the preset gradient, we simulated 100,000 bacteria of four evenly distribution phenotypes of different receptor gain $N$ or adaptation time $\tau $ or basal CheYp level ${Y}_{p0}$ . The simulation was performed in 3D with reflection boundary on $x=0\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}\mathrm{m}\text{}\mathrm{}\text{}\mathrm{x}=20\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}\mathrm{m}$ and with free boundaries in the other two directions ($y,z$). Initially all cells were placed on $x=2\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}\mathrm{m},\text{}\mathrm{y}=\mathrm{z}=0\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}\mathrm{m}$. In each time point with step $\Delta t=0.01s$, the internal states of every bacterium were updated according to their perceived gradient indexed by $g\left({x}_{i}\right)$ with spatial precision $\mathrm{\Delta}x=0.001\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$. And then the bacterial moving stat (run or tumble) was updated accordingly, which results in movement to the next time point. Most parameters used in the simulation follow previous paper (Saragosti et al., 2011), while the changed parameters were the upper limit of chemoattractant sensing $K}_{on}=1000\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{M$ (Fu et al., 2018); the lower limit of chemoattractant sensing $K}_{off}=1\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{M$ (Si et al., 2012); the diffusion coefficient of chemoattractant $D}_{s}=1000\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mu {\mathrm{m}}^{2}/\mathrm{s$ (Fu et al., 2018); the run speed of bacteria $v}_{0}=30\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}/\mathrm{s$ (Keller and Segel, 1971b); the adaptation time ${\tau}_{M}=1s$ (Waite et al., 2016); the motor number was set 5 as in Sneddon et al., 2012; while the basic CheYP level $Y{p}_{0}$ = $2.77\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{M}$, which corresponds to $T{B}_{0}\approx 0.27$. After 40 min of simulation, bacterial tracks were recorded with 10 fps for 20 min.
7. Estimation of $\chi $ from chemotaxis pathway
In order to compare the agentbased simulation results to the analytical OUtype model, we need to estimate the chemotactic ability $\chi $ from the parameters of chemotaxis pathway. Considering an simple case where perceived gradient of chemoattractant $g$ is a fixed value, the averaged drift velocity of a group of bacteria ${V}_{D}$ can be estimated by the parameters of bacterial chemotaxis pathway (Dufour et al., 2014): ${V}_{D}\approx \frac{{\tau}_{R0}^{`}}{1+{\tau}_{R0}/{\tau}_{M}}\frac{\left(1T{B}_{0}\right){v}_{0}^{2}Ng}{d}$ . In this equation, ${\tau}_{R0}$ is the run time in homogenous environment, ${\tau}_{R0}^{`}=d{\tau}_{R0}/df$ with $f$ is the free energy of the receptor, ${\tau}_{M}$ is the adaptation time, $d$ is the dimension (e.g. d is 3 in this case), $T{B}_{0}$ is the tumble bias in homogenous environment, and $N$ is the receptor gain.
On the other hand, given the expression of ${V}_{D}$ , comparing the two expressions of expected drift velocity in fixed gradient, we can deduce:
Using the expression Equation S20, we can deduce the required receptor gain $N$ from $\chi $, which is the case shown in the main text in Figure 3.
8. Agentbased simulation for different parameters
To investigate the bacterial behavior in moving gradient, we have simulated bacteria with $\chi $ ranging from 0.15 to 0.30 corresponding to receptor gain $N$ ranging from 14.5 to 25.4. Related results were shown in Figure 3 and discussed in the main text.
Receptor gain is not the only parameter that affects the chemotactic ability. We have also varied the adaptation time ${\tau}_{M}$ and basal CheYp level ${Y}_{p0}$ to investigate the effect of these parameters on bacterial behaviors with a moving chemoattractant profile of larger slope of perceived gradient, $g}_{1}=2.5\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}{\mathrm{m}}^{2$ . For both cases, the density profile $\rho \left(z\right)$ spatially ordered according to the variant parameter, and the expected velocities for call parameters ${V}_{d}\left(z\right)$ decrease from back to front. These results show that bacteria perform mean reversion flow for all these parameters and were spatially ordered according to phenotypes.
For the case of basic CheYp level ${Y}_{p0}$ ranging from 2.70 µM to 2.85 µM, the basic tumble bias $T{B}_{0}$ varies from 0.07–0.19, and the chemotactic ability $\chi $ decreases with ${Y}_{p0}$ from 0.16 to 0.25 mm^{2}min^{1} . The mean positions of subgroups increase with $\chi $ so that it decreases with ${Y}_{p0}$ . Unlike the receptor gain $N$, the basic tumble bias also increases the diffusion coefficient $D$. Analogized to the OUtype model, varying ${Y}_{p0}$ changes the reversion rate $\chi g$ and the noise level $\u03f5=\sqrt{2D}$ at the same time. Thus, the reversion rate decreases with ${Y}_{p0}$ and the width of the density profile $\sigma $ increases with the reversion rate $r$ (Appendix 1—figure 8A–D).
For the case of adaptation time ${\tau}_{M}$ ranging from 1 to 10 s, mean positions and width of the density profile show nonmonotonic dependence on ${\tau}_{M}$ . This phenomenon shows that adaptation time has a tradeoff effect on bacterial chemotaxis. As Frankel et al., 2014 have discussed, the increasing adaptation time enforces chemotactic ability, while it weakens the sensitivity of bacteria to the chemoattractant gradient changes. Such dual effect was also reflected on the reversion rate defined by the slope of ${V}_{D}$ (Appendix 1—figure 8E–H).
9. Modified Langevintype model for accumulated signals
As discussed in the main text, the reversion flow of individuals in the migrating group relies on a decreasing driving force from back to front. In the case of bacterial chemotaxis, such decreasing force comes from the gradient of perceived signal of chemoattractant, which is consumed by all the agents. Such mechanism can be generalized to the systems of trailfollowing migration (e.g. ants), where the signal is secreted by agents instead of being consumed. In this case, the signal concentration ${S}_{s}(x,t)$ decreases from back to front, because more agents have passed the back than front, thus have secreted more signal on the back. Such dynamics can be modeled by similar Langevin equations and signal accumulations as discussed in Section 1:
where $\gamma =3\ast {10}^{13}\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{M}\cdot {\mathrm{s}}^{1}\cdot \mathrm{c}\mathrm{e}\mathrm{l}{\mathrm{l}}^{1}$ is the secretion rate of the signal ${S}_{s}$ and $\beta =0.001{s}^{1}$ is the decay rate of it. Other terms and parameters are identical to equations in Section 1.
Then, we solve Equation S21 numerically by applying 10^{6} particles with normally distributed response ability ${\chi}_{i}$ of mean value of 1 and s.d. of 0.3 (Appendix 1—figure 10C). Results show that a coherent migrating group is formed spontaneously with constant migration speed ∼6.7 µm/s (Appendix 1—figure 10B). In the moving coordinate, the ‘effective global force’ ${S}_{s}\left(z\right)$ decreases from the back to front as predicted. The particles with different response ability ${\chi}_{i}$ are ordered arranged (Appendix 1—figure 10D), while particles with larger χ locate more in the front.
Data availability
All data generated or analysed during this study are included in the manuscript and supporting files.
References

Effect of Amino Acids and Oxygen on Chemotaxis in Escherichia ColiJournal of Bacteriology 92:121–129.https://doi.org/10.1128/jb.92.1.121129.1966

Applications of Microfluidics in Quantitative BiologyBiotechnology Journal 13:e1700170.https://doi.org/10.1002/biot.201700170

Collective Memory and Spatial Sorting in Animal GroupsJournal of Theoretical Biology 218:1–11.https://doi.org/10.1006/jtbi.2002.3065

Selforganization and collective behavior in vertebratesAdvances in the Study of Behavior 32:1–75.https://doi.org/10.1016/S00653454(03)010015

Chemotaxis: the role of internal delaysEuropean Biophysics Journal 33:691–693.https://doi.org/10.1007/s002490040426z

Limits of Feedback Control in Bacterial ChemotaxisPLOS Computational Biology 10:e1003694.https://doi.org/10.1371/journal.pcbi.1003694

Direct Correlation between Motile Behavior and Protein Abundance in Single CellsPLOS Computational Biology 12:e1005041.https://doi.org/10.1371/journal.pcbi.1005041

Advanced methods of microscope control using muManager softwareJournal of Biological Methods 1:e10.https://doi.org/10.14440/jbm.2014.36

From individual to collective behavior in bacterial chemotaxisSIAM Journal on Applied Mathematics 65:361–391.https://doi.org/10.1137/S0036139903433232

Evolution Transforms Pushed Waves into Pulled WavesThe American Naturalist 195:E87–E99.https://doi.org/10.1086/707324

Model for chemotaxisJournal of Theoretical Biology 30:225–234.https://doi.org/10.1016/00225193(71)900506

Traveling bands of chemotactic bacteria: A theoretical analysisJournal of Theoretical Biology 30:235–248.https://doi.org/10.1016/00225193(71)900518

Differential fitness returns in relation to spatial position in groupsBiological Reviews of the Cambridge Philosophical Society 69:187–206.https://doi.org/10.1111/j.1469185x.1994.tb01505.x

The spatial organization of microbial communities during range expansionCurrent Opinion in Microbiology 63:109–116.https://doi.org/10.1016/j.mib.2021.07.005

BookFrom Individual to Collective Behavior: The Role of Memory and Diversity in Bacterial NavigationYale University.

A Gradually Slowing Traveling Band of Chemotactic BacteriaJournal of Mathematical Biology 19:125–132.https://doi.org/10.1007/BF00275935

The structure and function of fish schoolsScientific American 246:114–123.https://doi.org/10.1038/scientificamerican0682114

Fundamental theoretical aspects of bacterial chemotaxisJournal of Theoretical Biology 41:201–208.https://doi.org/10.1016/00225193(73)901136

On the propagation theory for bands of chemotactic bacteriaMathematical Biosciences 20:185–189.https://doi.org/10.1016/00255564(74)900789

Numerical Study of Formation and Propagation of Travelling Bands of Chemotactic BacteriaJournal of Theoretical Biology 46:189–219.https://doi.org/10.1016/00225193(74)901477

Pathwaybased meanfield model for Escherichia coli chemotaxisPhysical Review Letters 109:0–4.https://doi.org/10.1103/PhysRevLett.109.048101

Front propagation into unstable statesPhysics Reports 386:29–222.https://doi.org/10.1016/j.physrep.2003.08.001

Non‐genetic diversity modulates population performanceMolecular Systems Biology 12:895.https://doi.org/10.15252/msb.20167044

Behavioral Variability and Phenotypic Diversity in Bacterial ChemotaxisAnnual Review of Biophysics 47:595–616.https://doi.org/10.1146/annurevbiophys062215010954

Multiscale Models of TaxisDriven Patterning in Bacterial PopulationsSIAM Journal on Applied Mathematics 70:133–169.https://doi.org/10.1137/070711505
Article and author information
Author details
Funding
National Key Research and Development Program of China (2018YFA0903400)
 Xiongfei Fu
National Natural Science Foundation of China (32071417)
 Xiongfei Fu
CAS Interdisciplinary Innovation Team (JCTD201916)
 Xiongfei Fu
Guangdong Provincial Key Laboratory of Synthetic Genomics (2019B030301006)
 Xiongfei Fu
Strategic Priority Research Program of Chinese Academy of Sciences (XDPB1803)
 Xiongfei Fu
National Natural Science Foundation of China (11804355)
 Yang Bai
National Natural Science Foundation of China (31770111)
 Yang Bai
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
The authors acknowledge C Liu, for sharing E. coli strains and plasmids; S Huang for help with the microfluidics; T Emonet for help with singlecell tracking and agentbased simulations; T Hwa, F Jin, L Luo, J Long for insightful discussions. This work was financially supported by the National Key R&D Program of China (Grant No. 2018YFA0903400), National Natural Science Foundation of China (Grant No. 32071417), CAS Interdisciplinary Innovation Team (Grant No. JCTD201916), Guangdong Provincial Key Laboratory of Synthetic Genomics (Grant No. 2019B030301006), Strategic Priority Research Program of Chinese Academy of Sciences (XDPB1803) to XF, National Natural Science Foundation of China (Grant Nos. 11804355 and 31770111) to YB.
Copyright
© 2021, Bai et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 1,278
 views

 307
 downloads

 10
 citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Immunology and Inflammation
 Microbiology and Infectious Disease
Mortality and morbidity from tuberculous meningitis (TBM) are common, primarily due to inflammatory response to Mycobacterium tuberculosis infection, yet the underlying mechanisms remain poorly understood. We aimed to uncover genes and pathways associated with TBM pathogenesis and mortality, and determine the best predictors of death, utilizing wholeblood RNA sequencing from 281 Vietnamese adults with TBM, 295 pulmonary tuberculosis (PTB), and 30 healthy controls. Through weighted gene coexpression network analysis, we identified hub genes and pathways linked to TBM severity and mortality, with a consensus analysis revealing distinct patterns between HIVpositive and HIVnegative individuals. We employed multivariate elasticnet Cox regression to select candidate predictors of death, then logistic regression and internal bootstrap validation to choose best predictors. Increased neutrophil activation and decreased T and B cell activation pathways were associated with TBM mortality. Among HIVpositive individuals, mortality associated with increased angiogenesis, while HIVnegative individuals exhibited elevated TNF signaling and impaired extracellular matrix organization. Four hub genes—MCEMP1, NELL2, ZNF354C, and CD4—were strong TBM mortality predictors. These findings indicate that TBM induces a systemic inflammatory response similar to PTB, highlighting critical genes and pathways related to death, offering insights for potential therapeutic targets alongside a novel fourgene biomarker for predicting outcomes.

 Microbiology and Infectious Disease
Patients present a wide range of clinical severities in response severe acute respiratory syndrome coronavirus 2 infection, but the underlying molecular and cellular reasons why clinical outcomes vary so greatly within the population remains unknown. Here, we report that negative clinical outcomes in severely ill patients were associated with divergent RNA transcriptome profiles in peripheral immune cells compared with mild cases during the first weeks after disease onset. Protein–protein interaction analysis indicated that earlyresponding cytotoxic natural killer cells were associated with an effective clearance of the virus and a less severe outcome. This innate immune response was associated with the activation of select cytokine–cytokine receptor pathways and robust Th1/Th2 cell differentiation profiles. In contrast, severely ill patients exhibited a dysregulation between innate and adaptive responses affiliated with divergent Th1/Th2 profiles and negative outcomes. This knowledge forms the basis of clinical triage that may be used to preemptively detect highrisk patients before lifethreatening outcomes ensue.