Authors: Henry H. Mattingly, Thierry Emonet
Categories: Biological Sciences, bacterial chemotaxis, bet-hedging, collective behavior, nongenetic diversity, nongenetic inheritance, 435
Source: Proceedings of the National Academy of Sciences of the United States of America
Natural selection enriches a population with the best-adapted phenotypes. But collective behaviors can also shape a population’s phenotype composition, even without selection. We study this in the context of collective migration of bacteria, in which the spatial arrangement of individuals by their chemotaxis abilities determines which individuals keep up with the migrating group. Since this spatial organization is environment-dependent, we find that a slow loss of low-performing phenotypes enables an isogenic population to nongenetically adapt its phenotype composition to migrate in changing environments. An important part of this adaptation strategy is the time scale on which individuals with new phenotypes are produced. Nongenetic inheritance provides a way to tune this time scale and may be widespread among microbes.
Keywords: collective behavior, nongenetic diversity, nongenetic inheritance, bet-hedging, bacterial chemotaxis
Collective behaviors enable groups to perform tasks that an individual cannot, but they require coordination (1–4). Therefore, even when collectives can accommodate individual heterogeneity, they often reduce it by removing some phenotypes from the population (1, 5). At the same time, isogenic populations with nongenetic diversity produce new phenotypes by reproducing (6–8), which can replace lost phenotypes. How do these processes of removal and production of phenotypes balance, and how do they affect collective performance of the population?
Collective migration of chemotactic bacteria is a powerful model system for studying the interaction of collective behaviors and growth. Groups of Escherichia coli can migrate collectively by consuming an attractant cue in their environment, creating a traveling gradient that they chase (5, 9–27). Experiments have recently shown that this process leads to directed range expansion of the population and can increase the overall population’s growth rate, relative to undirected expansion, by maintaining a larger fraction of the population in exponential growth (21, 28).
Even isogenic populations of E. coli exhibit a range of behavioral phenotypes (5, 29)—quantified, for example, by the fraction of time a cell spends tumbling, or tumble bias (30, 31)—which climb gradients at different speeds (32). We showed that these diverse phenotypes can migrate together by spatially sorting along the traveling gradient they create, matching individual chemotaxis ability to the local gradient steepness so that all individuals migrate at approximately the same speed (5). However, this compensatory mechanism is The concentration of attractant at the back of the group drops below the detection limit of the cells’ receptors (21, 33–35), causing low-performing phenotypes located there to slowly fall behind (5). Thus, collective migration gradually removes these low performers.
Those experiments were performed in liquid, where cells with low tumble bias climb static attractant gradients the fastest (32, 36). However, bacteria encounter other natural environments (37), such as soil and animal hosts. In porous media, in particular, experiments indicate that cells with intermediate tumble bias navigate best (12, 16, 18, 19). Theory (38) and agent-based simulations (39) suggest that this is because cells that tumble infrequently (low tumble bias) get stuck in traps, whereas cells that tumble too frequently (high tumble bias) move very little during each run. Cells with intermediate tumble frequency balance these two competing effects to diffuse and perform chemotaxis effectively. Thus, individual and collective chemotaxis in environments of different porosity require different swimming behaviors.
If collective migration reduces diversity over time, how can populations leverage their multiple behavioral phenotypes to migrate effectively in both liquid and porous environments? One solution may be the generation of individuals with new phenotypes. Recent experiments have shown that a substantial portion of behavioral variability in E. coli populations is generated by imperfect inheritance of phenotype upon cell division (8). While daughter-cell phenotypes are correlated with those of their mothers, they are not identical. Thus, cell divisions provide a way to replenish phenotypes lost during collective behaviors.
Here, we use a generalized Keller–Segel model (10) to quantitatively study how collective migration and growth balance to determine the population’s phenotype composition and performance in multiple environments. By specifically removing low performers in each environment, we predict that collective behaviors interact with growth to dynamically adapt the population composition to the collective task. We find that a key parameter in this balance is the level of phenotypic inheritance. While higher inheritance enriches the population with high performers for the current environment, which increases migration speed, it comes at the cost of slower responsiveness to new environments. Finally, this trade-off leads to an optimal level of inheritance that maximizes average migration speed in varying environments and enables a diverse population to outperform a nondiverse one. Thus, the tension between collective behavior and individuality may generally balance with growth to enable diverse populations to dynamically and nongenetically shape their own compositions of phenotypes to match the demands of collective tasks.
To quantitatively study how the conflict between collective behavior and individuality shapes migrating bacterial populations, we extended the classic Keller–Segel model (10) to include multiple behavioral phenotypes, finite sensitivity for attractants, and cell growth (Methods). In this model, cells of different tumble bias, TB, climb gradients with drift speeds that depend on their chemotactic performance χ(TB) (36), and they move randomly with diffusivity μ(TB).
For liquid environments, we previously showed that, despite differences in chemotactic abilities, cells of different phenotypes, TB, are able to migrate together by spatially organizing themselves by performance, χ(TB) (5) (Fig. 1A; reproduced from ref. 5). However, the mapping from phenotype TB to performance χ(TB) depends on the porosity of the environment in which the cells swim. The logical, but previously unknown, consequence of this is that phenotypes should spatially organize in an environment-dependent manner during group migration. Furthermore, by placing low performers at the back, we expect that collective migration removes whichever phenotypes perform poorly in the current environment.
Fig. 1. Collective migration differentially removes behavioral phenotypes in an environment-dependent manner. (A) Experimental data from ref. 5. Two isogenic populations with different distributions of phenotypes (TB) travel together in liquid by spontaneously sorting themselves within the migrating group. High-performing (low-TB; blue) cells are in front, where the traveling gradient is shallow, and low-performing (high-TB; red) cells are in back, where the traveling gradient is steep, equalizing chemotactic drift speed within the group. Cells located at the back are at a higher risk of falling behind. Reprinted from ref. 5, which is licensed under CC BY 4.0. (B) The chemotactic performance of each phenotype, quantified by χ, depends on the physical properties of the environment in which it swims. The highest-performing phenotypes in liquid (green) and porous (orange) environments are different. (C) Schematic: Since cells in the group sort by performance χ, the spatial organization of phenotypes TB is different in different environments. (D) Since χ(TB) and spatial organization are environment-dependent, so is the differential loss of phenotypes. Phenotype composition P(TB) of the migrating group in two simulations, in liquid (green; D, Left) and porous (orange; D, Right) environments. Black lines indicate the initial phenotype compositions. Each colored line corresponds to a snapshot in time, with lighter color corresponding to later time. Group migration removes different phenotypes in each environment.
To test these predictions, we modeled chemotaxis in porous media using a biophysical model for cell diffusion μ(TB) developed by Licata et al. (38). We took χ(TB) to be closely related to μ(TB) (Fig. 1B and SI Appendix) because both arise from cells’ run-and-tumble motility (21, 22, 35). To isolate the effect of collective migration on the dynamics of the migrating group’s phenotype composition, we first conducted simulations in which the cells did not grow.
Simulations confirmed our predictions (schematic in Fig. 1C and simulation results in SI Appendix, Figs. S1 and S2): Phenotypes in the migrating groups were organized differently in liquid and porous environments, but they always sorted by performance, χ(TB). As a result of this different spatial organization, the phenotypes that fell behind the migrating group differed in each case, with high-TB cells falling behind in liquid and low-TB cells falling behind in porous media. To quantify this environment-dependent loss of cells based on chemotactic performance, we derive in SI Appendix a semianalytical expression for how fast cells fall behind, or “leak” out of, the migrating
Here, N(t) is the number of cells traveling in the group at time t, c(t) is the time-dependent migration speed, ϵ is the cells’ sensitivity for the attractant, and 〈χ−μ|z0〉 is the average chemotactic performance among phenotypes located at the (moving) back edge of the group, z0. Individual phenotypes fall behind at rates set by their relative abundance at z0 (SI Appendix). Due to the spatial organization of phenotypes by performance, this expression reveals that the total flux with which cells fall behind the migrating group is set by the lowest-performing phenotypes traveling.
Over time, this differential loss of phenotypes shifts the composition of the group from its initial condition toward individuals that are high-performing in each environment (Fig. 1D). These composition changes do not occur due to differences in growth or death. Instead, collective migration leaves behind phenotypes from the traveling group that perform poorly in the current environment.
To investigate how growth and loss of cells balance during collective migration, we needed a model of cell growth. Modeling growth with diversity requires specifying which phenotypes are produced when a cell divides—that is, the extent to which daughter cells nongenetically inherit the phenotypes of their mothers (8). Imperfect inheritance may arise due to stochasticity in partitioning chemotaxis-related proteins and organelles between the two daughter cells at the time of cell division (6, 7). To make minimal assumptions, we modeled production of phenotypes as a Gaussian, auto-regressive process (Methods and SI Appendix), similar to recent work studying nongenetic inheritance of growth rate (40). Since tumble bias TB lies between zero and one, we defined TB=1/(1+eF) and, instead, modeled production of the transformed phenotype F=log (1−TBTB).
In this model, when the population is grown in batch culture—i.e., when it is not migrating—its distribution of phenotypes F approaches a Gaussian with population mean 〈F〉 and variance σ2. This description captures experimental measurements of the TB distribution in a wild-type strain of E. coli cells (SI Appendix, Fig. S3). Upon division, mother cells with phenotype F′ produce daughter cells with phenotypes F given
where η is a Gaussian random variable with mean 〈η〉=0 and variance 〈η2〉=σ2 (1−ϕ2), and ϕ is the correlation between mother and daughter phenotypes. When ϕ=0, daughter-cell phenotypes are chosen at random from the batch-culture distribution—there is no memory. As ϕ increases, daughters’ phenotypes are more similar to their mothers’ phenotypes, increasing the characteristic time, τ, that it takes a perturbation to the population’s phenotype composition to τ−1=2 r (1−ϕ) (SI Appendix). While others have studied the role of growth rate, r, and other model parameters on group migration speed in populations containing one phenotype (35), this inheritance parameter is unique to populations that exhibit nongenetic diversity. We focus on it here because we will see that it has significant effects on group migration dynamics.
Simulating the model of collective migration with growth, we found that the migrating group not only approaches a stable migration speed and population size, but also a stable phenotype composition (Fig. 2 A and B). Leakage of cells selectively removed low-performing phenotypes at the back of the group, but growth did not specifically replace the phenotypes that were lost. Therefore, the phenotype composition of the migrating group became dynamically enriched with phenotypes that perform well in the current environment, relative to the batch-culture distribution. This composition shift depended on the phenotypic inheritance parameter ϕ: With increasing ϕ, the composition shifted more because the traveling cells produced daughters that were more similar to themselves and more dissimilar to the cells that fell behind. In the limit of ϕ→1, the migrating population would eventually purify to the single phenotype that performs best in the current environment.
Fig. 2. Growth and leakage balance to adapt the group’s composition to the environment, depending on nongenetic inheritance of phenotype. (A and B) Steady-state phenotype compositions P(TB) in groups migrating through liquid (A; green throughout) or porous (B; orange throughout) environments, for varying values of mother–daughter phenotype correlation ϕ, which increases from dark to light lines. Black solid phenotype composition in batch culture, when not migrating. (C) Larger shifts in phenotype composition away from the batch-culture composition lead to faster steady-state migration, css. (D) Theory captures dependence of steady-state migration speed css on group composition. Line with circular markers is the predicted css (Eq. 3) plotted versus css in simulations. Each marker is a simulation of a population with different value of ϕ. Solid lines without markers are the predictions of Eq. 3 if 〈χ−μ|z0〉, the mean performance at the back of the group, is replaced with 〈χ−μ〉, the average performance over the entire migrating group. Predictions in this panel used 〈χ−μ|z0〉 and 〈χ−μ〉 at steady state from simulations. Dashed black line is y=x. (E) Dynamic changes in phenotype composition circumvent the trade-offs between migrating in liquid and porous environments faced by any population with fixed composition. Black normalized css of nondiverse populations with different phenotypes (TB=0.013 at the bottom right; TB=0.4 at the top left) in a porous environment versus their css in a liquid environment. Green liquid specialist (highest css in liquid). Orange porous specialist (highest css in porous). Black generalist (highest minimum css among the two environments). Colored css of diverse populations in each environment for varying values of ϕ (shown are ϕ=0 in blue through ϕ=0.8 in red, in increments of 0.2). (F) Faster steady-state migration css of diverse populations trades off against the time it takes the population’s phenotype composition to adapt to a new environment. Inverse adaptation time (τadapt)−1 (or adaptation rate) and css from simulations were normalized to lie between zero and one. Larger values of ϕ lead to faster migration in both environments, but slower adaptation rate.
These shifts in phenotype composition were accompanied by increases in steady-state migration speed (Fig. 2C). While this seems intuitive, migration speed is set by how fast the population consumes the attractant—i.e., the size of the traveling population—not necessarily the performance of its individuals (refs. 10 and 35 and SI Appendix). Since the traveling gradient is generated by the cells’ consumption, the cells essentially make the gradient as steep as necessary to maintain the migration speed set by the pace of their consumption (SI Appendix). Therefore, for composition shifts to increase migration speed, the traveling population’s size must increase at steady state.
To quantitatively understand this, we extended an approach from a recent theoretical study of collective migration with a single phenotype (35) to include multiple phenotypes. Detailing this approach in SI Appendix, we find that the steady-state migration speed css
Eq. 3 shows that css is set by the diffusivity of the attractant, Ds, and the average chemotactic performance, 〈χ−μ|z0〉, among phenotypes located at the back edge of the group, z0. The latter quantity appears because it affects the loss of cells at the back (Eq. 1), which is exactly the net flux at z0 (SI Appendix). Higher inheritance ϕ increases 〈χ−μ|z0〉, thus reducing the rate at which all phenotypes fall behind and increasing the number of cells traveling and the migration speed at steady state. By extracting steady-state css and 〈χ−μ|z0〉 from simulations with varying values of ϕ, Fig. 2D shows that Eq. 3 accurately captures how migration speed depends on the phenotypes traveling. These effects could not be captured by instead using the mean chemotactic abilities of the entire group—i.e., 〈χ−μ〉 (Fig. 2D)—highlighting the importance of the lowest performers in the traveling group.
An additional consequence of these dynamic composition changes is that a diverse population with inheritance can avoid trade-offs faced by any nondiverse population. Fig. 2E shows the steady-state migration speeds css of a range of nondiverse populations in each environment, scaled by the maximum migration speed in each (black line). The specialist phenotypes in agar and liquid are those that migrate the fastest in each environment, but they perform less well in the alternative environment. The generalist phenotype has the highest minimum migration speed across the two environments. The convex shape of the black line suggests that no diverse population with a fixed composition of phenotypes can outperform the nondiverse populations—i.e., by being closer to the (1, 1) point on that plot. However, dynamic adaptation enables diverse populations to outperform the nondiverse populations by enriching themselves with specialists for the current environment, relative to the batch-culture composition. This is shown by comparing the highest-performing homogeneous population (black dot in Fig. 2E) to a diverse population with the same mode phenotype (colored dots in Fig. 2E) and different values of inheritance ϕ. As ϕ increases, the diverse population approaches the (1, 1) corner of the plot, accessing a region of performance space that is inaccessible to any homogeneous population (black line in Fig. 2E). Thus, by dynamically adapting its own composition, the population is able to effectively migrate through multiple environments that require different phenotypes.
Finally, while higher inheritance leads to faster steady-state migration, we found that it also slows down the population’s response to a new environment. To quantify this trade-off, we extracted the relaxation time of migration speed near steady state for varying values of ϕ (Methods). While other quantities could be used to quantify the adaptation time scale, the relaxation time near steady state is independent of the population composition before the environment change, and it is the limiting time scale over which the composition adjusts to a new environment. Note that there is also a faster time scale in this Upon encountering a new environment, the phenotypes in the population spatially reorganize according to their new values of χ, and some phenotypes immediately fall behind. This initial change in population composition occurs in the time it takes the migrating group to travel its own width.
Normalizing the adaptation rates and migration speeds to lie between zero and one and plotting them against each other, we visualized this trade-off (Fig. 2F). In both environments, adaptation rate τadapt−1 decreased with steady-state migration speed css with an approximate shape of τadapt−1∝1/css. The intuition for this trade-off can be understood by making an analogy to an Ornstein–Uhlenbeck (OU) process. In an OU process, a particle diffuses in an energy-potential well, and the relaxation time of the particle’s position depends on its mass and the drag of the medium. A larger-mass particle with a shorter relaxation time is also less responsive to external forces. Likewise, lower inheritance shortens the relaxation time of the group’s phenotype composition, but makes it less responsive to the differential loss of phenotypes resulting from collective migration.
The time it takes a migrating population to adapt its composition becomes important when it encounters changing environments. How does the time scale of the environment constrain the ability of diverse populations to adapt, and when is diversity useful? To investigate this, we simulated diverse and nondiverse migrating populations that periodically encounter liquid and porous environments after fixed amounts of time T (Fig. 3). We varied the level of inheritance ϕ of the diverse populations since it controls the relaxation time of their phenotype compositions. Then, we defined population performance as the time-averaged migration speed across the two environments, after many environment switches.
Fig. 3. Diverse populations nongenetically adapt to migrate through varying environments and can outrun a single-phenotype, generalist population in slowly-varying environments. (Top) Schematic of the simulation. Populations with varying mother–daughter phenotype correlation ϕ were simulated migrating through liquid and porous environments that switch after a fixed time T. (Middle) After many environment switches, the stationary, time-averaged migration speed across the two environments cdiverse was computed and normalized by that of a single-phenotype, generalist population cgeneralist. Migration speeds (colors) are plotted against the relaxation time of the population’s composition in batch culture, τ=1/(2r (1−ϕ)), scaled by the cell-division time 1/r, and the environment duration T, also scaled by 1/r. Note that τ increases as ϕ increases. For a plot with ϕ on the y axis, see SI Appendix, Fig. S6. White line separates two qualitatively different regimes. Insets show the phenotype composition at the end of each liquid (green) and porous (orange) environment. Black dashed lines in Insets indicate the generalist phenotype. When the environment varies rapidly (r T<∼3), the population narrows its composition around the generalist phenotype, and larger ϕ (larger τ) increases migration speed by more effectively narrowing the phenotype composition. When the environment varies slowly (r T>∼3), the diverse population adapts its composition to each environment. An optimal level of inheritance ϕ (or τ) maximizes migration speed in this regime (black dashed line, Middle).
Plotting the average migration speeds of the diverse population as a function of the environment duration, T, and the population’s composition relaxation time in batch culture, τ, we found that the average migration speed exhibited two qualitatively different regimes (Fig. 3 and Methods). In environments that changed quickly compared to the cells’ doubling time (r T<∼3), higher inheritance ϕ (longer composition-relaxation time in batch culture) monotonically increased the average migration speed. In this regime, the population’s phenotype composition narrowed around the generalist phenotype, and higher inheritance helped the population maintain the generalist phenotype when cells divided (Fig. 3, Left). Increasing diversity by increasing σ reduced average migration speed in this regime because the population was less able to narrow its phenotype distribution around the generalist phenotype (SI Appendix, Fig. S6 A and B). The nondiverse, generalist population always migrated fastest in this regime.
In slowly changing environments (r T>∼3), the population’s phenotype composition tracked the environment. An optimal level of inheritance ϕ maximized the average migration speed by balancing the trade-off between maximally adapting the composition to each environment and the need to do so before the environment switched (Fig. 3, Right). Once the period of each environment was sufficiently long, the diverse population outperformed the nondiverse one by shifting and narrowing its composition around the specialist phenotype in each environment. Increasing diversity in this regime made the “highs higher” and the “lows lower”: Populations without inheritance (ϕ=0) performed worse with more diversity, but populations with sufficient inheritance performed better with more diversity (SI Appendix, Fig. S6 A and B). Increasing diversity effectively increases the range of phenotypes to which the population has access, but inheritance is required to maintain the high-performing ones. These results suggest that nongenetic inheritance should be as large as the environment time scales allow. Thus, differential leakage paired with inheritance of phenotypes transiently and nongenetically adapts the population composition to migrate effectively through varying environments.
Here, we investigated how growth in a nongenetically diverse population balances with loss of phenotypes during collective behaviors and how this balance impacts population performance. We used collective migration of chemotactic bacteria as a model system, where the demands of migration change in environments of different porosity. We realized that collective migration causes differential loss of phenotypes in an environment-dependent Migration excludes whichever individuals perform poorly in the current environment. Our central finding is that, with growth and phenotypic inheritance, this differential loss dynamically enriches the population phenotype composition with high performers for the current environment, thus enhancing group migration speed in multiple environments. Growth of diverse phenotypes requires specifying how much daughters “remember” the phenotypes of their mothers, and we found that this nongenetic inheritance controls a general trade-off between larger composition shifts relative to the batch culture and slower responsiveness to new environments. In varying environments, this trade-off resulted in an optimal level of inheritance that maximized the average migration speed and enabled a diverse population to outrun a nondiverse one composed of the generalist phenotype.
Differential loss of phenotypes during group migration resembles the effects of natural Individuals that are “better adapted” to the current environment become enriched in the traveling group over time. However, these composition changes were not due to differences in growth or death, by assumption. Instead, differential loss of migrating individuals emerged from their spatial arrangement by chemotaxis performance. This distinguishes our results from past work on the role of nongenetic diversity in well-mixed, growing populations (41–60) and populations undergoing undirected range expansion (61–63). Spatial structure in past work has affected population growth or evolutionary success by governing which individuals access limited nutrients in dilute conditions (64, 65), which populations occupy finite habitats (22, 23), or the ability of populations to exchange metabolites (66). Our results point to another way that spatial structure can shape the composition of a by dictating which individuals participate in a collective behavior. Since this spatial structure is environment-dependent, it enables migrating populations to dynamically shape their compositions to migrate effectively in multiple environments.
The dynamic nature of these nongenetic changes in phenotype composition also differentiate our results from past work by us (65) and others (47, 53, 54, 57, 67). In Frankel et al. (65), diversity enabled a population to avoid trade-offs among multiple tasks by having phenotypes in the population that performed well at each task. These tasks were performed independently by each individual cell, without collective behavior. While mutations allowed the population to evolve its standing distribution of phenotypes, this distribution was fixed for a given genotype during performance of a task. Our results here imply that, over cell divisions, diverse populations can do Growth and selection-like effects (either natural selection or the kind described here) dynamically shift the population’s neutral distribution of phenotypes toward those that are high performing at each task, without mutations, as the task is being performed. Likewise, approaches to studying the fitness effects of diversity that geometrically analyze Pareto fronts or optimize phenotype distributions without considering how the phenotype distribution dynamically responds to selection or selection-like effects will miss this (47, 53, 54, 57, 67, 68). The problem is that the population can have a different phenotype composition in each environment, giving it more flexibility than these procedures allow. As a result, the population is not necessarily limited to the convex hull of nondiverse populations’ performance across environments. This, in turn, affects the conditions in which a diverse population performs better than a nondiverse, generalist population.
Predicting changes in phenotype composition requires knowing how new phenotypes are produced. We modeled this as resulting from imperfect inheritance of swimming phenotypes at cell division (40, 42, 48, 69–73), consistent with recent measurements (8). Stochastic switching is another known mechanism for producing new phenotypes (41, 43, 49–52, 74, 75), such as in the case of persisters (50, 51, 55, 76). Since these mechanisms produce phenotypes independently of the cells’ environment, we expect that they are both subject to the trade-off that increasing the “susceptibility” of the population’s phenotype composition to selection-like pressures comes at the cost of slower responses to new pressures. This trade-off underlies past results showing that an optimal level of inheritance or an optimal switching frequency maximizes population growth in changing environments (41, 42, 56). Despite the very different context, we found that migrating groups of bacteria exhibit the same qualitative behavior. The existence of an optimum suggests that, on longer time scales, the population can evolve the level of nongenetic inheritance to match the time scales of environmental changes.
Emergent adaptation of phenotype composition during collective behaviors could have important effects on a range of biological processes and could interact with natural selection. For example, the adaptation studied here could play a role during competition among multiple comigrating species of bacteria (22, 23, 28). More generally, group migration with spatial ordering of phenotypes is not limited to bacterial It can arise whenever individuals with diverse motility phenotypes consume or deplete a signal in the environment that they then chase. Self-generated gradients drive collective migration of eukaryotic cells during development, wound healing, and cancer metastasis (77–82), in which a leader–follower structure emerges (83–85), and, thus, adaptation mediated by spatial organization of motility phenotypes could affect these processes. A similar kind of nongenetic adaptation could potentially occur during biofilm formation, in which the collective behavior confers a selective advantage to phenotypes that produce costly biofilm-associated materials (86). Collective behaviors may generally adapt the phenotype compositions of diverse, growing populations to perform a variety of collective tasks.
Our generalized Keller–Segel model consists of the following partial differential equations (PDEs):
Here, ρ(TB)=ρ(x,t,TB) is the density of cells with tumble bias TB at location x and time t, ρ=ρ(x,t)=∫01ρ(TB′) dTB′ is the total density of cells, and s=s(x,t) is the concentration of attractant. χ(TB) and μ(TB) are the phenotype-dependent chemotactic (36) and diffusion coefficients characterizing cells’ chemotactic drift speed and random diffusive motion, respectively. Models for χ(TB) and μ(TB) in liquid and porous environments are described in SI Appendix. Ki and Ka are the dissociation constants of the cells’ receptors for the attractant when in the inactive and active forms (33), which roughly set the lower and upper limits of concentration that the cells can detect. We assume all phenotypes grow with rate r, and ρc is the carrying capacity. We chose ρc sufficiently large that it had negligible effects on growth in the migrating group. R(TB,TB′) encodes how divisions of a cell with phenotype TB′ contribute to production of cells with phenotype TB. Ds is the diffusivity of the attractant, and k and Km characterize consumption of the attractant by all cells. These PDEs are paired with boundary conditions stipulating that, far ahead of the wave, there are no cells, and the concentration of attractant is s∞ (i.e., there is no attractant gradient ahead of the group): ρ(x→∞,t,TB)=0, ∂xρ(x→∞,t,TB)=0, s(x→∞,t)=s∞, and ∂xs(x→∞,t)=0. Variants of this model that simulate individual agents with a single phenotype (21, 87) and multiple phenotypes (29) have been studied in the past, finding that the large number of cells in these migrating populations [typically more than 108 cells/mL (5, 21)] justifies modeling them with a continuous density ρ(TB).
In the growth term, R(TB,TB′) =2 P(TB|TB′)−δ(TB−TB′) is a kernel that maps divisions of mother cells to production of daughters. P(TB|TB′) is the distribution of offspring phenotypes TB given the mother’s phenotype TB′ that is implied by the dynamics in Eq. 2 and by the mapping from TB to F. The factor of two accounts for production of two daughter cells upon division (whose phenotypes may be correlated with each other; SI Appendix), and the Dirac delta function accounts for removal of the mother cell from the population upon division. We assumed that nutrients are abundant in the migrating group and distinct from the attractant (equivalent to large carrying capacity ρc), such that the migrating cells grow exponentially (21, 35). For simplicity, we assumed that all phenotypes divide with the same rate r, and we neglected changes in a cell’s behavioral phenotype within its lifetime (8).
The code in this work built on previous code from ref. 5. Simulations were performed by converting the PDEs into a system of ordinary differential equations by discretizing space x. Spatial derivatives were computed by using central differences. Integration forward in time was done by using an explicit, fourth-order Runge–Kutta scheme (88, 89). The initial condition was taken to be a steep sigmoidal function in ρ(x,TB,t=0)=OD0 P0(TB)1+exp (30(x−xc)xc), where OD0 is the initial density of cells, P0(TB) is the initial composition, and xc sets the extent of the initial inoculation. We used xc=〈χ−μ〉 s∞k OD0, where the average in 〈χ−μ〉 is with respect to P0(TB). This makes xc roughly the characteristic length scale (width) of the wave that would be formed by the inoculated cells. For this choice of initial condition, nearly all phenotypes in the initial inoculation join the traveling group. P0(TB) was always the batch-culture distribution. In simulations without growth, we used OD0=6; with growth in liquid, we used OD0=3; and with growth in porous media, we used OD0=0.5. The spatial grid size was dx=xc/50, using knowledge of the relevant length scale of the problem. Time steps were dt=dx2/(4 μmax), where μmax was the largest diffusion coefficient in the problem (including Ds).
Simulations were done over a moving segment of space with no-flux boundary conditions on both sides. Several conditions were used to decide when to move the window of space being simulated. A new segment was simulated the attractant concentration at the right boundary dropped below a threshold; the peak cell density (and, hence, problem length scale) changed by more than 15%; or the fraction of traveling cells that had the smallest or largest value of TB fell below a threshold (and could be removed from the simulation) or increased above a threshold (meaning that the new phenotypes would imminently be produced by growth). If any of these conditions were met, the simulation was stopped, space far behind the wave peak was discarded, and the spatial region of the simulation was extended ahead of the direction of motion. The spatial grid size and time step were updated during this step, with xc=〈χ−μ〉 s∞k ρmax, where the average was over phenotypes currently traveling, and ρmax was the peak value of total cell density at the end of the previous simulation. The number of cells of each phenotype traveling was found by integrating over space from the location where the attractant s first reached zero behind the wave to space far ahead of the group. Then, the cell density and attractant profiles were interpolated onto the new spatial grid, and the simulation was resumed.
With cell growth, phenotypes that weren’t present at the start of the simulation could be produced. In this case, if the simulation of a spatial block stopped because the fraction of edge phenotypes increased above a threshold due to growth, then the range of phenotypes in the simulation was extended. Simulations were implemented with phenotype defined as F. The grid size of F was chosen to be 1/10th of the SD of the offspring phenotype distribution in order to resolve the growth matrix R(F,F′).
For the traveling populations to reach steady state in Fig. 2, group migration was simulated for 15 h. The composition-relaxation time was determined from the time course of migration speed as it approached steady state. Once the speed was within 10% of the steady-state value, the speed profile was fit with an exponential function to extract the time scale of the slowest decaying mode, which was used to quantify the adaptation time. To compare theory to simulations, N(t,TB) in simulations was computed by integrating cell density over the wave from z=z0 to z→∞ (where s(z0)=Km Ki defines the back of the wave; SI Appendix, Steady State Migration Speed with Attractant Diffusion).
Simulations in periodically varying environments were done by first simulating migration through one environment for the fixed environment duration. Then, to make minimal assumptions about how transitions between environments occurred, we took the cell density and attractant profiles at the end of the simulation in one environment and used those profiles as the initial conditions for a simulation in the next environment. The spatial grid size was adjusted for the expected change in speed in the new environment. In order for the simulation to reach a steady state, environments were alternated 15 times this way, for a total of 30 simulations for each parameter pair (ϕ, T). The time-averaged speed in each condition was computed by dividing the total displacement over the last two environments (one liquid, one porous) by twice the duration of each environment, T. The white line in Fig. 3 that separates the regime of fast (r T<3) and slow (r T>3) environment switching was determined by finding, for each ϕ, the value of r T at which the average speed had a local minimum.
All simulations were performed on Yale’s High-Performance Computing Clusters.
We thank Jude Ong for collecting the tracking data in SI Appendix, Fig. S4, and we thank Paul Turner and Sujit Datta for helpful discussions. This work was supported by NIH Awards R01GM106189 (to H.H.M. and T.E.), R01GM138533 (to H.H.M. and T.E.), and F32GM131583 (to H.H.M.).
Simulation-code data have been deposited in GitHub (https://github.com/emonetlab/ks) (90).
Simulation-code data have been deposited in GitHub (https://github.com/emonetlab/ks) (90).