Spatial self-organization resolves conflicts between individuality and collective migration

Spatial self-organization resolves conflicts between individuality and collective migration ARTICLE DOI: 10.1038/s41467-018-04539-4 OPEN Spatial self-organization resolves conflicts between individuality and collective migration 1,6 1,7 1,2 1 3 1,8 4,5 1,2 X. Fu , S. Kato , J. Long , H.H. Mattingly ,C.He , D.C. Vural , S.W. Zucker & T. Emonet Collective behavior can spontaneously emerge when individuals follow common rules of interaction. However, the behavior of each individual differs due to existing genetic and non- genetic variation within the population. It remains unclear how this individuality is managed to achieve collective behavior. We quantify individuality in bands of clonal Escherichia coli cells that migrate collectively along a channel by following a self-generated gradient of attractant. We discover that despite substantial differences in individual chemotactic abilities, the cells are able to migrate as a coherent group by spontaneously sorting themselves within the moving band. This sorting mechanism ensures that differences between individual chemo- tactic abilities are compensated by differences in the local steepness of the traveling gradient each individual must navigate, and determines the minimum performance required to travel with the band. By resolving conflicts between individuality and collective migration, this mechanism enables populations to maintain advantageous diversity while on the move. 1 2 Department of Molecular, Cellular, and Developmental Biology, Yale University, New Haven, CT 06520, USA. Department of Physics, Yale University, New Haven, CT 06520, USA. Institute of Synthetic Biology, Shenzhen Institutes of Advanced Technology, Chinese Academy of Sciences, Shenzhen 518055, 4 5 China. Department of Computer Science, Yale University, New Haven, CT 06520, USA. Department of Biomedical Engineering, Yale University, New Haven, CT 06520, USA. Present address: Institute of Synthetic Biology, Shenzhen Institutes of Advanced Technology, Chinese Academy of Sciences, Shenzhen 518055, China. Present address: Department of Molecular Biotechnology, Graduate School of Advanced Sciences of Matter, Hiroshima University, Higashi-Hiroshima, Hiroshima 739-8530, Japan. Present address: Department of Physics, University of Notre Dame, Notre Dame, IN 46556, USA. These authors contributed equally: X. Fu and S. Kato. Correspondence and requests for materials should be addressed to T.E. (email: thierry.emonet@yale.edu) NATURE COMMUNICATIONS (2018) 9:2177 DOI: 10.1038/s41467-018-04539-4 www.nature.com/naturecommunications 1 | | | 1234567890():,; ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04539-4 ells and larger organisms exhibit collective behaviors that quantify the distribution of phenotypes in the band (Fig. 1b and 1,2 5 are often advantageous to the participating individuals . Supplementary Fig. 1). Approximately 2 × 10 clonal E. coli cells CMany such collective behaviors dynamically emerge when grown in M9 glycerol medium (M9 salts, glycerol, and casamino a large number of individuals follow the same rules to interact acids; Methods) were introduced with fresh medium into the 3,4 with each other and the environment . Prominent examples are device and concentrated at the end of the channel by cen- 4–6 bird flocks and the collective migration of bacteria along trifugation (Methods). Following centrifugation, sequential bands 7–11 8,9,12 channels and on agar plates . At the same time, pheno- of cells collectively migrated along the channel at different but typic differences among even genetically identical individuals are nearly constant speeds (Fig. 1c), presumably consuming different a ubiquitous feature of biology . Phenotypic diversity can lead to compounds within the undefined media, as demonstrated in early useful leader–follower structures within a traveling group. For studies . example, in migrating neural crest cells and in fish shoals, many E. coli cells navigate by alternating straight “runs” with organisms may follow a few more informed individuals .In “tumbles” that randomly reorient their swimming direction microbial communities, maintaining diversity in the population (Fig. 1d). By transiently suppressing tumbles whenever attractant can enable bet-hedging strategies to survive uncertain environ- signal increases, they perform a biased random walk that allows 13,15–18 20 ments and resolve trade-offs . However, heterogeneity can them to move toward higher concentrations of attractant . In the also be disruptive, as is the case in simulated swarms where non- absence of a gradient, the fraction of time a cell spends tumbling aligners tend to be purged from the swarm . This raises a —its tumble bias (TB)—remains approximately constant and dilemma: although phenotypic diversity provides advantages, it therefore can be used as a quantitative measure of the phenotype also tends to reduce coordination. of the cell. Importantly, using the same strain and microfluidic One of the simplest cases of collective behavior is exhibited by channel depth, we previously demonstrated that the TB is a bacteria: clonal populations of motile Escherichia coli cells col- strong determinant of chemotactic performance in liquid: lower lectively migrate when placed at high density at the bottom of a TB cells drift significantly faster up a static gradient than higher 8–11 22 tube filled with nutrients . This collective behavior is mediated TB cells . To quantify the distribution of phenotypes in the by the well-characterized chemotaxis system , which enables the isogenic population that was introduced in the device, a low bacteria to follow chemical gradients, in this case generated by density of cells was loaded into the microfluidic device without their consumption of attractant present in the medium (Fig. 1a). centrifugation and individual cells were tracked to determine 21,22 However, populations of E. coli exhibit substantial cell-to-cell their TB, as previously described . TB was broadly distributed 13,21 variability in their swimming phenotypes and hence che- in the population with some cells tumbling < 10% of the time (i.e., motactic abilities, even when all cells are genetically identical . TB < 0.1) and others > 50% of the time (Fig. 1e black), consistent 13,21,22 How bacterial populations manage phenotypic heterogeneity to with previous studies . Given the functional consequences still allow coordinated collective migration remains largely of this non-genetic diversity, how can the same population of cells unknown, mainly because of the difficulties in measuring cellular migrate together as a coordinated group, as shown in Fig. 1c? behavior at both the collective and the individual levels in the To answer this question, we first considered whether all same experiment . phenotypes or only a subset of them traveled in each band. We Although the migration of traveling waves or “bands” of bac- used pressure valves to capture one band of cells at a time in the teria has also served as a classic model for the theoretical study of wide chamber of our device (Fig. 1b and Supplementary Fig. 1b). 10,23 emergent phenomena and pattern formation in biology , the After trapping cells in the wide chamber, it was perfused with effect of non-genetic diversity on this process has scarcely been fresh media to homogenize the environment and dilute the cell examined. Previous studies examined how two populations may density. We verified that perfusion of the wide chamber did not 24,25 travel together ; however, it was assumed that within each affect the distribution of TBs (Supplementary Fig. 3). Dilution population all of the individuals were identical. The mechanisms enabled us to track individual cells. Homogenization ensured that by which a continuum of phenotypes can achieve coherent cells had adapted back to a uniform environment and were not migration have not been investigated. responding to an attractant gradient when we measured their Here we used a microfluidic system that enables precise TBs. The distribution of TB was shifted toward lower TB in both quantitative measurements at the individual and collective scales traveling bands compared with the original distribution (Fig. 1e), to study the interplay of diversity and collective bacterial suggesting that it was more difficult for high TB cells to migration. Our central finding is that within the traveling band, participate in collective migration. Selection against high TB cells cells spontaneously sort themselves such that their chemotactic was stronger in the faster band (Fig. 1e, red) than in the slower abilities are matched to the local gradient steepness, enabling one (Fig. 1e, yellow). Cell density and number also varied between diverse cells to travel together with the same drift speed. the two bands (Fig. 1c), suggesting that there were interdepen- Extending the classic Keller–Segel model of traveling bands to dencies between the speed of the group, its size, and the diversity account for diversity predicts this spatial sorting and qualitatively of the individuals able to migrate with the group. We periodically recapitulates the experimental results. Our second finding is a tracked cells after they were trapped and diluted in the wide novel mechanism that reduces the rate at which cells fall off the chamber and found that the original TB distribution was recov- back of the band: when attractant consumption depends on local ered after growth (Fig. 1f). Thus, selection of low TB cells by the oxygen, oxygen limitation in the center of the band increases the collective migration was not due to genetic heterogeneity. In gradient of attractant at the back, helping cells there keep up. addition, it is unlikely that cell growth affected the TB selection Together, these two mechanisms enable populations of bacteria to while cells were traveling in a band, because the duration of the maintain diversity while migrating as a group. experiment (30 min) was shorter than the cell doubling time (~ 55 min, Supplementary Fig. 4). Results Collective migration selects against high TB cells. To determine Cells of diverse chemotactic abilities migrate as a group.To the relationship between the number of cells in the band, the quantify collective behavior and diversity in the same experiment, band speed, and diversity, we switched from casamino acids to a we designed a microfluidic device consisting of a long channel to defined M9 glycerol buffer containing aspartate (Asp) as the observe the traveling band , followed by a large chamber to only limited chemoattractant (Methods). In this condition, a 2 NATURE COMMUNICATIONS (2018) 9:2177 DOI: 10.1038/s41467-018-04539-4 www.nature.com/naturecommunications | | | NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04539-4 ARTICLE ab Group migration Attractant Cell density Channel Chamber cd Individual diversity 720 Band 2 Band 1 Tumble bias 0 3.5 7 0.1 0.2 0.3 0.4 Position (mm) ef 01234 5 Time elapsed (hrs) 6 6 4 4 2 2 0 0 0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 0.5 Tumble bias Tumble bias Fig. 1 Collective migration of a phenotypically-diverse clonal population. a When concentrated at the bottom of a nutrient channel, motile E. coli cells emerge from the high cell density region and travel in bands along the channel by following gradients of attractant produced by their consumption. b Microfluidic device used to quantify the band migration and the phenotypic diversity within the band. Control gates along the channel (black vertical lines) are initially open (top) and later closed to capture different bands of cells in the observation chamber (bottom), where single cells are tracked to quantify the distribution of phenotypes within the band (Supplementary Fig. 1). c Time-lapse imaging of E. coli cells expressing the fluorescent protein mRFP1 showing the collective migration of bands in M9 glycerol medium (M9 salts, glycerol, and casamino acids; Methods). In this undefined medium, −1 −1 8 several bands emerge that travel at different speeds (red: 0.68 mm min , yellow: 0.23 mm min ) . We verified that labeling cells did not affect band speeds nor tumble bias distributions (Supplementary Fig. 2). Scale bar, 0.6 mm. d The tumble bias (color)—average probability to tumble—of individual 21, 22 cells was quantified by tracking a cell for 2 min in a uniform environment (no gradient) and detecting tumbles (black dots) as previously described . Scale bar, 200 μm. e Collective migration selects against high TB cells. Tumble bias distributions from the first (red) and second (yellow) bands (n = 3), and from the population that was introduced in the device (black; n = 3). f The tumble bias distribution of the cells in the first wave (red in e) gradually shifted back toward the original distribution (black) during growth in the perfused chamber. Fresh M9 glycerol medium was supplied every 30 min. The TB distribution was measured every 20 min single band formed (Fig. 2a) and its speed could be tuned by density of cells in the band increased (Fig. 2c, d), and the dis- changing the concentration of supplemented Asp (Fig. 2b). To tribution of TB within the band shifted toward higher TB measure band speed and density, cells expressing mRFP1 were (Fig. 2e). In general, collective migration selected against high TB mixed with unlabeled cells at ratios of 1:20, 1:50, or 1:100 (for 50, cells, with selection being stronger in faster bands (Fig. 2e, f). It is 100, and 200 μM Asp, respectively), and their positions were noteworthy, however, that diversity was not eliminated—all 11,26 detected at various time points (Fig. 2a and Methods ). Band bands still exhibited a range of TBs. Thus, although collective speed, number of cells in the band, and density profiles were migration selected against high TB cells, it was still possible for a stable over time with a slow decay due to cells falling off at the diverse group to travel together. back of the band (Fig. 2b, c, d and Supplementary Fig. 5). Although there were variations across experimental replicates due Extending the Keller–Segel model to account for diversity.To to variations in the number of cells introduced in the device, the better understand how collective migration selects TB, we relationship between Asp concentration and band parameters was extended the classic Keller–Segel mathematical model describing consistent. Specifically, as the concentration of Asp increased, the traveling bands of bacteria to include phenotypic diversity speed of the band decreased (Fig. 2b), the number and peak (Methods Eqs. 1–3). In this model, cells consume the diffusible NATURE COMMUNICATIONS (2018) 9:2177 DOI: 10.1038/s41467-018-04539-4 www.nature.com/naturecommunications 3 | | | Time (sec) PDF PDF ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04539-4 a d 9 ×10 0 3.5 7 –1.5 –1 –0.5 0 0.5 1 1.5 Position (mm) Relative coordinate, z (mm) be *p = 0.004 0.5 *p = 0.0076 0.4 0.3 n = 5 0.2 n = 4 0.1 n = 5 Simulation 0 0 50 100 150 200 0.1 0.2 0.3 0.4 0.5 Asp conc. (µM) Tumble bias c f *p = 0.0017 *p = 0.0034 1/2 10 Simulation 1/4 50 100 150 200 0.1 0.2 0.3 0.4 0.5 Asp conc. (µM) Tumble bias Fig. 2 Relationship between band speed, density, and phenotypic diversity. a Time-lapse coordinates of cells (black dots) traveling in M9 glycerol buffer (Methods) with 200 µM aspartate. Only one band forms in each experiment. Cells (1:100) were labeled with mRFP1 and their coordinates detected (Methods). Scale bar, 0.6 mm. In remaining panels, colors are aspartate concentration in the buffer: 200 µM (blue; 1:100 cells labeled), 100 µM (yellow; 1:50 cells labeled), 50 µM (red; 1:20 cells labeled) (Methods). b Band speed decreased with aspartate concentration in the buffer. Circles: experiments in which tumble bias was measured, used in e, f. Diamonds: tumble bias was not measured. Dashed: simulations (Methods and Supplementary Fig. 5). p-value: one-sided, two-sample t-test assuming unequal variances, t-value= 4.1, df = 5.4 between 50 μM and 100 μM, t-value= 4.2, df= 3.8 between 100 μM and 200 μM. c Number of cells traveling in the band increased with aspartate concentration. Circles, diamonds, and dashed: same as in b. p-value: one- sided, two-sample t-test assuming unequal variances, t-value= − 3.9, df = 6.5 between 50 μM and 100 μM, t-value=− 4.6, df = 6.4 between 100 μM and 200 μM. d Cell density profiles averaged over time and experiments (Supplementary Fig. 5). Line: average over n = 5 (red), n = 4 (yellow), n = 5 (blue) replicate experiments. Shade: SD. For each experiment, nine profiles were measured at 1.3 min intervals. Dashed: simulations. e Tumble bias distribution of the cells that traveled with the band. Lines: average over experiments; shading: SD; dashed lines: simulations. f Ratio between the distribution of tumble bias in the traveling band (Fig. 2e, colored lines) and that of the original population (Fig. 2e, black) quantifies the enrichment of tumble bias. Lower aspartate concentrations enrich more for lower tumble bias. Solid lines: mean of the measurements from experiments; shading: standard error of the mean; dashed: simulations attractant Asp, which generates a traveling gradient that the cells indeed, for TB = 0 the cell is just diffusing and χ = 0; 21,22 follow by biasing their random walk (Fig. 3a). Faster consump- Methods) . Thus, in a gradient of attractant, cells with higher tion, more cells, or less Asp all lead to a faster traveling gra- TB do not diffuse as much and climb slower than cells with lower dient . The motion of a phenotype i depends on two parameters: TB. The dependency of f on Asp concentration has been char- its effective diffusion coefficient, μ = μ(TB ), which results from acterized as well . Moreover, we conducted high-performance i i the cells’ random walk, and its chemotactic coefficient, χ = χ liquid chromtography (HPLC) experiments to verify that the (TB ), which quantifies how effective that phenotype is at biasing chemotaxis response to Asp dominated that to amino acids 27–29 its motion to follow the perceived amount of Asp, f. Theory secreted as byproducts of Asp metabolism. Experiments con- and tracking of individual E. coli cells swimming in a static gra- ducted with mutants lacking the oxygen receptor aer or both dient of α-methylaspartate (non-metabolizable analog of Asp) aer and tsr indicated that aerotaxis was not essential and Tar have shown that μ and χ are decreasing functions of the TB (note response to Asp was sufficient for band migration (Methods that for very low values of TB < ~ 0.05, χ increases with TB; and Supplementary Fig. 6). 4 NATURE COMMUNICATIONS (2018) 9:2177 DOI: 10.1038/s41467-018-04539-4 www.nature.com/naturecommunications | | | Migration speed –1 Time (min) Cell number in the band (mm min ) –1 Fold enrichment Cell density (ml ) PDF NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04539-4 ARTICLE ab Total TB 0.1~0.2 TB 0.3~0.4 TB 0.2~0.3 TB 0.4~0.5 ×10 df df C = 2 = 1 8 dz dz 6 20 Falling off Traveling –1.5 –1 –0.5 0 0.5 1 Relative coordinate, z (mm) cd –1 TB 0.1 TB 0.2 10 20 TB 0.3 TB 0.4 0.2 0.1 Individual Group –2 –1.5 –1 –0.5 0 0.5 1 –1.5 –1 –0.5 0 0.5 1 Relative coordinate, z (mm) Relative coordinate, z (mm) Fig. 3 Mathematical modeling predicts a mechanism for consistent collective migration of diverse phenotypes. a Collective migration of diverse phenotypes at the same speed is made possible by the spontaneous spatial ordering of individual phenotypes within the band such that each individual’s chemotactic df ability is matched to the local gradient steepness . The proportion of better performers (larger χ, lower TB; red) should be enriched where the gradient is dz shallower (front), whereas the proportion of weaker performers (smaller χ, higher TB; blue) should be enriched where the gradient signal is steeper (back). The position where the perceived gradient steepness is maximum (dashed border of the gray regions) determines the highest tumble bias able to travel with the band. Cells in the gray region slowly fall out of the band. b Simulated density profiles of cells migrating in 200 μM aspartate (the same simulation is shown in Fig. 2def, blue dashed) show sorting based on tumble bias. c Chemotactic coefficient χ(z) (blue) defined by the phenotype whose density df profile peaks at position z and perceived gradient steepness (black). Red symbols correspond to the location of the peak cell density for individual dz df phenotypes. d Spatial sorting enables consistent migration velocity for traveling phenotypes. The migration velocity, χðÞ z , of the phenotype whose density dz profile peaks at position z gradually decreases toward the back of the band until the gray region is reached. In the gray region the migration velocity falls off more rapidly, preventing the high TB phenotypes located there from staying in the band To metabolize Asp, E. coli consumes oxygen . Introducing a Spatial sorting as a mechanism for inclusive migration.An fluorescent oxygen sensor in the M9 glycerol buffer revealed important feature of the experiments reproduced by the simula- that oxygen availability is reduced in the center of the traveling tions is the increasing selection against high TB cells as the band where cell density is high (Supplementary Fig. 7a-f). This amount of Asp is reduced (Fig. 2f). What mechanism enables cells results in a dependency of the average consumption rate of Asp with diverse chemotactic abilities to collectively migrate together on cell density (Supplementary Fig. 7h). We modeled this effect as one band, and what controls the upper bound on the TB such that the Asp consumption rate depended linearly on oxygen among those able to migrate together? For every phenotype that concentration and constrained the related parameters by travels with the band at constant speed c, the flux of cells must be measuring oxygen and Asp consumption rates in batch cultures approximately invariant in time and equal to the chemotactic flux (Supplementary Fig. 7g). For simplicity, we ignored possible minus that due to diffusion. Focusing on the partial differential phenotypic diversity in the Asp consumption rate, as well as the equation for the cell density of phenotype i and switching to the possible dependence of the diffusion coefficient on cell density , moving reference frame z = x − ct (here x and t are absolute which was found previously to be negligible in similar position and time), two predictions emerge (Methods Eqs. 7–8). experiments . We also omitted possible contributions of First, the traveling phenotypes must distribute themselves spa- hydrodynamics and physical interactions between cells, which tially within the band so that differences in local signal strengths, df can become important when bacteria swarm over surfaces . the slope of the black line in Fig. 3a, compensate for differences dz The resulting model (Methods and Table 1) qualitatively in the chemotactic abilities of each phenotype, χ , such that reproduced the main features of our experiments, including the df c ¼ χ ðÞ z , where z is the position of peak density of phenotype i dz i dependency on Asp concentration of the band speed (Fig. 2b), cell i (Fig. 3a). Therefore, this spatial sorting places the better per- number and density (Fig. 2c, d), TB distribution (Fig. 2e), formers (higher χ , lower TB ) in the front of the band, where the i i phenotypic selection (Fig. 2f), and average Asp consumption rate gradient is shallower and more difficult to follow (Fig. 3a, red), per cell as a function of cell number in the band (Supplementary and the weaker performers (lower χ , higher TB ) at the back, i i Fig. 7h). NATURE COMMUNICATIONS (2018) 9:2177 DOI: 10.1038/s41467-018-04539-4 www.nature.com/naturecommunications 5 | | | Chemotactic coefficient 2 –1 (mm min ) Perceived signal gradient –1 df /dz (mm ) Migration speed –1 (mm min ) –1 Cell density (ml ) Perceived signal f (a.u.) ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04539-4 where the gradient is steeper (Fig. 3a, blue). Furthermore, a result of spatial sorting. Due to experimental limitations, we could second prediction is that the gradient will reach a maximal not measure the TB distributions of all four strains in each of the steepness (Fig. 3a dashed border of gray zone), determining the experiments reported in Fig. 4a-e. To better quantify the rela- weakest phenotype that can travel (lowest χ , upper bound on tionship between peak separation and difference in ⟨TB⟩ value, TB ). Thus, we see the interplay between individuals i and the we instead mixed pairs of populations using only fluorescent effect of the community on the available resource f. strains induced with different levels of anhydrotetracycline (aTc). Analysis of simulations confirmed these analytical predictions For each pair, we measured the TB distributions right before df loading cells in the device, and then measured the distance (Fig. 3b, c, d). The steepness of the perceived signal , which dz between the fluorescence peaks in the resulting traveling band emerges dynamically from the cells’ consumption, peaked at the (Fig. 4f). This confirmed that there is a monotonic relationship back of the traveling band (low z) and decayed toward the front between peak separation and difference in ⟨TB⟩ values. Therefore, (high z) (Fig. 3c black). In contrast, the position z of the peak cells of various phenotypes appear to spontaneously sort them- density of phenotype i increased with its chemotaxis coefficient χ selves along the traveling band according to their TB, enabling (Fig. 3c blue), revealing an ordering of the phenotypes within the them to migrate collectively despite phenotypic differences. migratory band (Fig. 3b, c). Multiplying the two together gave a nearly constant velocity throughout the band, thus providing an explanation of how the various phenotypes might travel together Discussion (Fig. 3d). The rightmost points in Fig. 3c (blue line) and Fig. 3d How do organisms maintain collective behavior despite the (black line) correspond to the phenotype with the maximum potential conflicts created by phenotypic diversity among indi- chemotactic coefficient in the band. Ahead of that location, there viduals? We studied this question using traveling bands of che- are no more peaks in the cell density of any phenotype; however, motactic E. coli, which collectively migrate at the same speed there are cells due to diffusion. At the back of the band, cells with despite differences in chemotactic abilities of individuals in the TB higher than the predicted upper bound rapidly fell off of the band. Our key result is that spontaneous spatial organization of band (Fig. 3a, d gray zones). Importantly, these predictions phenotypes within a traveling band helps resolve the conflicts emerged from just the dynamics of cell density (Methods Eq. 1) between phenotypic diversity and collective migration. By and therefore hold true irrespective of whether oxygen is included matching individual abilities to the local difficulty of the navi- in the model. gation task within the band, this sorting mechanism ensures Comparing simulations with and without the oxygen- consistent migration speed across the band. This process also dependent consumption rate revealed that the oxygen depen- determines the minimum chemotactic performance required to dency reduces leakage of the highest-TB cells located at the back keep up with the band, therefore explaining how diversity can of the band (Supplementary Fig. 8b). The higher concentration of become limited by collective behavior. Thus, the mechanism oxygen at the back relative to the center of the band locally reported here enables a continuum of phenotypes to migrate increases the rate of consumption, and hence the slope of the coherently. traveling gradient of Asp, helping the cells there stay longer with In the traveling band, there is always a slow leakage of cells off the band (Methods Eq. 9). Thus, oxygen dependency has a similar the back of the band because of the finite sensitivity of the cells for cohesive effect as the secretion of a self-attractant, which also 38–41 11,25,37 the attractant they are chasing . High TB cells, in particular, helps reduce the leakage of cells . which are localized at the back of the wave, are at risk of falling off. We discovered that this leakage can be reduced (but not Cells in the band are spatially sorted by TB. To experimentally eliminated) if the consumption rate of the attractant is lower in test the prediction of spatial sorting in the band by phenotype, we the center of the band than at the back, where the consumption measured the relative position of two populations of cells with rate determines the local gradient steepness and chemotactic drift. different mean TB within the traveling band. The distribution of In our case, this arises because Asp consumption depends on TB in the population was controlled by manipulating the level of oxygen, which becomes limited in the center of the band where expression of the phosphatase CheZ, which deactivates the che- cell density is high. This mechanism provides an alternative to motaxis response regulator CheY (Supplementary Fig. 9a). We other mechanisms known to reduce cell leakage, such as the 11,25,37 generated multiple populations with different TB distributions of secretion of an attractant by the traveling cells . Note that varying mean TB (⟨TB⟩) (Fig. 4a). In each population, we labeled the spontaneous sorting mechanism discussed above helps com- 1 in 50 cells of the same genetic background and induction level pensate for differences in chemotactic abilities, irrespective of the by ectopically expressing either mRFP1 or YFP. The inducer was presence of such auxiliary mechanisms (oxygen or self- washed away when cells were resuspended in buffer before attractant). starting the migration experiment. TB distributions are stable for Traveling bands of bacteria have been studied for decades since 21 8,9 more than an hour in buffer . When cells from the lowest (red) Julius Adler’s experiments in capillary tubes . Adler reported the and highest (cyan) TB distributions were mixed in equal parts formation of multiple traveling bands in complex media; we and introduced in the device, spontaneous spatial order emerged, observe the same in casamino acids (Fig. 1) and expect multiple with cells from the low (high) TB distribution located at the front bands to be able to form when multiple consumable attractants (back) of the traveling band (Fig. 4b, c and Supplementary Fig. 9b are present. Within a migrating band the cells respond to the for replicates). The distance between the peaks of the density traveling gradient, some parts of which can be fairly steep. profiles of the two populations remained nearly constant over the Therefore, we expect the instantaneous TB of an individual cell to duration of the experiment, indicating that the two populations be dynamically changing depending on its direction of motion 11,29 traveled together at the same speed (Fig. 4e). We verified that the and position within the gradient, as previously reported . Here distance between the peak densities of the two populations was we showed that phenotypic (intrinsic) differences in adapted TB stable in longer experiments (Supplementary Fig. 9d). Mixing between cells contribute significantly to spatial structure within populations with closer TB distributions caused the peaks of the the traveling band. In future studies, we will separate the con- two traveling populations to be closer (magenta and green in tributions of phenotypic and dynamic diversity to group struc- Fig. 4a, d and Supplementary Fig. 9c for replicates), suggesting ture. It will also be interesting to examine the contribution of that distance between peaks increases with difference in ⟨TB⟩ as a initial conditions, dimensionality, and growth (necessary to 6 NATURE COMMUNICATIONS (2018) 9:2177 DOI: 10.1038/s41467-018-04539-4 www.nature.com/naturecommunications | | | NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04539-4 ARTICLE a d 9 ×10 40 ng/ml aTc 6 ng/ml aTc 4 ng/ml aTc 2 ng/ml aTc –1.5 –1 –0.5 0 0.5 1 1.5 0 0.2 0.4 0.6 0.8 1 Tumble bias Relative coordinate, z (mm) b e 560 600 2 468 2 468 Position (mm) Peak position (mm) c f ×10 0.4 0.3 0.2 0.1 0 0 –1.5 –1 –0.5 0 0.5 1 1.5 0 0.1 0.2 0.3 0.4 Relative coordinate, z (mm) Difference in mean tumble bias Fig. 4 Phenotypes spontaneously order themselves along the traveling band according to tumble bias. a The distribution of tumble bias in the population can be controlled by manipulating the level of expression of the phosphatase CheZ, which deactivates the chemotactic response regulator CheY (Supplementary Fig. 9a). b Time-lapse coordinates (colored dots) of an equal mixture of low (red in a) and high (cyan in a) TB cells traveling in 200 µM aspartate M9 glycerol buffer (see Methods). Scale bar, 0.6 mm. c Corresponding density profiles (colors) together with total cell density (black). (Line: mean over n = 34 time points measured at 40 s intervals for one experiment; shading: SD; five replicates are in Supplementary Fig. 9b). d Same as in c, but for the magenta and green populations in a. Two replicates are in Supplementary Fig. 9c. e Peak positions as a function of time for the experiment in b. f The distance between the fluorescence intensity peaks of the two populations traveling together in a single band increases with the difference between the mean TB of the two populations. For each independent experiment, two populations labeled with different fluorescent proteins (mRFP1 or YFP) were induced using different aTc levels to obtain distributions with different mean tumble biases. Dots: average over n = 4 experiments; error bars are SD 41–43 maintain traveling cell density over long times ) to this the range of sensitivity was assumed to extend to vanishing process. concentrations, as in the original Keller–Segel model , which is Previous analysis of bacterial traveling bands assumed that the not biologically realistic . 10,11,37,44,45 population consisted of identical cells or at most Following depletion of local resources, the spatial self- 24,25 two phenotypes . Here we extended these studies by taking organizing mechanism described here could enable populations into account the continuum of phenotypes that is always present of bacteria to maintain diversity while traveling toward better 13,21,22 in a population . One study that considered how two environments. This diversity increases the probability that a phenotypes might travel together made the theoretical assump- phenotype well-suited to unexpected environments will be avail- tion that cells sense only the direction of the gradient, not its able if needed during travel until a destination is reached where magnitude . This assumption causes the peak densities of the growth can replenish the population. As the range of phenotypes two phenotypes in that model to coincide in space, contrary to allowed within a traveling group depends on the spatial profile of our experimental observations. In another theoretical study , the the traveling gradient, this mechanism introduces important cells did respond to the gradient magnitude, and the two phe- feedback between the environment, cellular metabolism, and notypes in the traveling solution were spatially separated. phenotypic diversity, which together generate spatial patterns of Although not discussed in the paper, the phenotype with the phenotypes according to functional capabilities. The same higher chemotactic coefficient is in the front in that solution, in mechanism might also enable different bacterial species to travel agreement with our sorting mechanism. However, in that model, together, thus enabling migration of small ecosystems. NATURE COMMUNICATIONS (2018) 9:2177 DOI: 10.1038/s41467-018-04539-4 www.nature.com/naturecommunications 7 | | | Time (sec) –1 PDF Cell density (ml ) Distance between Time (sec) peak positions (mm) –1 Cell density (ml ) ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04539-4 order Runge–Kutta routine for temporal integration (time step 0.08 s). We used Collective migration of eukaryotes resulting from traveling no-flux boundary conditions. The initial condition was ρðÞ x; 0 ¼ gradients of attractants generated by consumption or breakdown ðÞ =x of an attractant has recently been found to be more important ρ  PðÞ TB e for 0 < x <1.6 mm and ρ (x,0) = 0for x >1.6 mm. Here, ρ is i 0 0 i than previously believed because it enables cell migration over the initial cell density scale, determined from the total initial cell number ~ 2 × 10 ; x = 0.8 mm; and P(TB ) was obtained from experimental measurements 0 i much larger distances than migration along externally imposed (Fig. 2eblack). gradients . Being able to maintain diversity within the traveling Assuming near-constant wave speed, we rewrite Eq. 2 in the moving group could be important in the context of immunity and coordinate z=x–ct and integrate from − ∞ to + ∞ to obtain cancer. cA hi α ¼ : ð4Þ The types of interactions between collective behavior and N=a phenotypic diversity reported here might also be at play beyond microbiology and cell biology, in contexts where individuals in a Here, α is the average attractant consumption rate, N is the number of cells in the band, and a is the cross-sectional area of the channel. In the absence of group respond to the cumulative effect on the environment of the oxygen-dependent consumption of Asp, the average consumption rate in the individuals ahead. Whereas a bird monitors its neighbors to band would be constant across experimental conditions . However, oxygen- benefit from the collective information acquired by the flock, here dependent consumption makes the average consumption rate decreases with individuals monitor the environmental gradient to benefit from increasing cell density in the band, which is correlated with the number of cells in the band (Fig. 2c). As shown in Supplementary Fig. 7h, the value of ⟨α ⟩ the information accumulated in the environment by the band. In A calculated from experimental data using Eq. 4 decreases as the number of cells in both cases, there is a group “memory” in the form of a spatial the band increases, which is captured in simulations. structure that individuals respond to . In general, both types of To derive the result of spatial sorting analytically from this model, we first memory are probably available and utilized by groups, with their rewrite Eq. 1 in the moving coordinate z = x − t: relative importance determined by the specific biology of the ∂ρ ∂ ρ ∂ ∂f ðAÞ i i c ¼ μ  χ ρ ; ð5Þ organisms. i i i ∂z ∂z ∂z ∂z Methods Noting that for each phenotype i to be traveling with the group, its density profile dρ Mathematical model. We extended the classic Keller–Segel model to include the i must have a peak. Around the density peak z we must have ¼ 0 and dz z¼z effect of phenotypic differences in TB. The key variables in the model are the d ρ density ρ (x, t) of cells of phenotype i as a function of position x and time t, and the i <0. dz z¼z concentration of Asp A(x, t). Cells of phenotype i are characterized by their che- i Rewriting Eq. 5 as motactic coefficient χ and diffusivity μ . As Asp consumption depends on oxy- i i gen , we also model the amount of oxygen dissolved in the media O(x, t). The 2 2 ∂ ρ ∂ρ ∂ρ ∂fAðÞ ∂ f ðAÞ i i i ð6Þ μ ¼c þ χ þ χ ; parameters of the model are in Table 1. The time-dependent evolution reads: i 2 i i 2 ∂z ∂z ∂z ∂z ∂z ∂ρ ∂ ρ ∂ ∂f ðAÞ i i ð1Þ ¼ μ  χ ρ we then have at the density peaks z : i i i i ∂t ∂x ∂x ∂x 2 2 d ρ d f μ ¼ χ ρ <0 ð7Þ i 2  i i 2 dz dz 2 z¼z z¼z X i i ∂A ∂ A A ¼ μ  α ðÞ O ρ ð2Þ A A i ∂t ∂x K þ A Eq. 7 indicates that in the front of the band the perceived gradient must be shallow and become progressively steeper as z decreases toward the back (Fig. 3c black). dρ 2 Integrating Eq. 5 and using  ¼ 0, we obtain dz ∂O ∂ O O z¼z ¼ μ  α ρ þ κðÞ O  O ð3Þ O O i ex ∂t ∂x K þ O df c ¼ χ ð8Þ dz z¼z Eqs. 7, 8 together show that the cell density peaks z of each phenotype i are Eq. 1 represents the motion of cells due to diffusion and chemotaxis. fAðÞ¼ hi monotonically ordered according to their chemotactic coefficients χ . 1þA=K Mlog is the perceived signal which depends on the local Asp concentration Examining the effect of oxygen analytically, at the back of the wave the Asp 1þA=K ∂fAðÞ M dA concentration is small so that the chemotactic drift there is χ  χ . A, where K = 3.5 μM and K = 1000 μM represent the dissociation constants of 0 1 i ∂z i K dz Asp for the inactive and active conformations of the Tar receptors, and M =6is Rewriting Eq. 2 in the moving coordinate, assuming the diffusion term is 29,48 the receptor gain . The effective diffusion coefficient and the chemotaxis negligible , and assuming A ≫ K gives an expression for dA/dz. From this, the coefficient are functions of the microscopic parameters of individual cell swimming drift becomes ðÞ 1TB 22,29 v i behavior. We draw on previous work to model them as μ ¼ i 3 ð1ΘÞλ þ2D R;i rot ∂fAðÞ χ M χ  α ðÞ O ρ : k TB ð9Þ 0 i −1 i A i and χ ¼ μ (Supplementary Fig. 8a). In these expressions, v ~36 μms is ∂z K c i TB þTB i 0 0 i i the cell swimming speed, which when projected in two dimesnions (2D) corresponds to the average speed we measure in our quasi-2D device (~ 28 μms At the back of the band there are fewer cells and therefore more oxygen than in the −1 22 ); Θ = 0.16 is the directional persistence between successive runs ; D = 0.062 middle of the band. Thus, when the Asp consumption rate depends on oxygen, rot qffiffiffiffiffiffiffiffiffiffi −1 21,50 TB α (O) becomes larger at the back than the mean over the band, ⟨α ⟩.Asa i A A s is the rotational diffusion coefficient during runs ; λ ¼ ω is the R;i 1TB consequence, the drift at the back is higher than in the case without oxygen −1 rate of switching from the run state to the tumble state; and ω = 3.8 s is the dependence, slowing down the decay of the band. We verified this by running effective switching frequency . μ and χ are monotonically decreasing functions of i i simulations with and without oxygen dependence (Supplementary Fig. 8b). In the the TB and χ /μ ≈ 22 = k , apart for the range 0 ≤ TB ≲ TB = 0.05, where χ(TB ) i i i 0 i 0 i simulations without oxygen, we set ⟨α ⟩ equal to a constant value, corresponding is increasing as recently observed in experiments where individual cells were to the average consumption rate in the band α in the simulation with oxygen tracked in a static gradient of α-methylaspartate (non-metabolizable analog of dependence. This was intended to make the wave speeds in the two simulations Asp) . similar, eliminating the effect of different wave speeds on cell leakage rates. Eq. 2 represents the change in the concentration of Asp due to diffusion, with diffusivity μ , and to consumption, with half-max constant K and maximum rate A A Strains, growth conditions, and sample preparation. E. coli RP437 was used as α ðÞ O ¼ α 1  g þ g . Here, α is the base consumption rate, O is the A0 ex A A0 A A O the wild type strain for chemotaxis in this study. Cells were grown in M9 glycerol ex −1 −1 −1 external oxygen concentration, and g is the fractional reduction of Asp medium: M9 salts (6.78 g L Na HPO , 3.0 g L KH PO , 0.5 g L NaCl, 1.0 g L A 2 4 2 4 −1 −1 consumption rate at zero oxygen. Eq. 3 describes the time-dependent evolution of NH Cl), supplemented with 4 mL L glycerol, 0.1 % casamino acids, 1.0 mM magnesium sulfate, and 0.05% w/v polyvinylpyrrolidone-40 at 30 °C. Appropriate oxygen, with diffusion coefficient μ , maximum consumption rate α , half-max O O −1 −1 constant K , and supply of oxygen through the polydimethylsiloxane (PDMS) with antibiotics were supplemented (ampicillin 100 µg mL , kanamycin 50 µg mL , −1 the mass transfer rate κ. and chloramphenicol 25 µg mL ) when necessary to maintain plasmids. Eqs. (1–3) were integrated in MATLAB using second-order centered For Fig. 1, cells were collected at mid-exponential phase (approximately an differences for the spatial derivatives (mesh size 20 µm) and an explicit fourth- OD600 of 0.3) and washed twice with fresh M9 glycerol medium, then resuspended 8 NATURE COMMUNICATIONS (2018) 9:2177 DOI: 10.1038/s41467-018-04539-4 www.nature.com/naturecommunications | | | NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04539-4 ARTICLE Table 1 Model parameters Symbol Definition and value Reference M The receptor gain for aspartate, M= 6 Ref. K The dissociation constant to aspartate for the inactive conformation of the Tar receptor, K = 3.5 μM Ref. 0 0 K The dissociation constant to aspartate for the active conformation of the Tar receptor, K = 1000 μM Ref. 1 1 2 −1 55 μ The diffusion coefficient of aspartate molecules, μ = 500 μm s Ref. A A K The aspartate concentration at half-max of its consumption, K = 0.5 μM Ref. A A −12 −1 −1 α The maximum aspartate consumption rate, α = 9.3 × 10 μmol cell Supplementary Fig. 7g A0 A0 O The external oxygen level, O = 250 μM Ref. ex ex g The basal ratio of relative consumption rate at zero oxygen, g = 0.27 Supplementary Fig. 7g A A 2 −1 55 μ The diffusion coefficient of dissolved oxygen, μ = 2500 μm s Ref. O O −11 −1 −1 α The maximum oxygen consumption rate, α = 7×10 μmol min cell Supplementary Fig. 7g O O K The dissolved oxygen concentration at half-max its consumption, K = 1 μM Supplementary Fig. 7g O O −1 58 κ The oxygen transfer rate through ~ 0.5 cm of PDMS, κ = 0.02 s Ref. ρ Initial cell density, ρ ~ 2.87 × 10 cells per mL This study 0 0 x Length scale of initial cell density profile, x = 0.8 mm This study 0 0 −1 v Cell swimming speed, v ~36 μms This study Θ Directional persistence Θ = 0.16 Ref. −1 50 D Rotational diffusion coefficient during runs, D = 0.062 s Ref. rot rot k Parameter in the expression for χ /μ , k = 22 Ref. 0 i i 0 TB Parameter in the expression for χ /μ,TB = 0.05 Ref. 0 i i 0 in fresh M9 glycerol medium to concentrate cell density at an OD600 of 0.7. These The laminated layers were then cut out and the remaining ports were punched to cells were then gently loaded into the microfluidics chamber, which was make external connections with the channels. To reduce the evaporation of the maintained at 30 °C throughout the experiment. microfluidic device, the PDMS device was soaked overnight in Millipore-filtered To generate a single traveling band, experiments were conducted in M9 glycerol water at 50 °C. buffer: motility buffer (M9 salts, 0.01 mM methionine, 0.1 mM EDTA, 0.05% w/v The assembled PDMS devices were bonded to 24 × 50 mm glass coverslips −1 polyvinylpyrrolidone-40) supplemented with 4 mL L glycerol and the indicated (#1.5). The PDMS was cleaned with transparent adhesive tape (Magic Tape, amount of Asp. This buffer was used to wash and resuspend cells instead of the Scotch) followed by rinsing with (in order) isopropanol, methanol, and Millipore- complex growth medium (M9 glycerol medium) mentioned above. The same filtered water, air-drying between each rinse. The glass was rinsed the same way, M9 salts were used in the M9 glycerol medium and M9 glycerol buffer to minimize with acetone, isopropanol, methanol, and Millipore-filtered water. The PDMS was osmolality changes. RP437 is auxotroph for leucine, histidine, methionine, and tape-cleaned an additional time, and then the two pieces were placed in a plasma threonine, and therefore does not grow in M9 glycerol buffer. bonding oven (Harrick Plasma) under vacuum, gently laminated, and then baked To control the TB distribution, we used a ΔcheZ strain derivative of RP437 on an 80 °C hotplate for 15 min to establish a covalent bond. Devices were stored at containing a chromosomally-integrated copy of the phosphatase CheZ under the room temperature and used within 24 h. control of the inducible promoter and tetR (gift from Dr. Chenli Liu). CheZ dephosphorylates the response regulator CheY, which, when phosphorylated, induces the motors to switch and the cells to tumble. Thus, low CheZ results in Band formation and imaging. Washed cells were gently loaded into the device higher CheY-P and more tumbling, and vice versa. aTc was added overnight in the which was then centrifuged for 20 min at 700 g in a 30 °C environmental room to culture when indicated to release the repressor TetR from the cheZ promoter concentrate cells at the end of the chamber (Supplementary Fig. 1b). After spin- region, inducing the expression of CheZ in this strain. To color-code strains, ning, the microfluidic device was placed on an inverted microscope (Nikon Eclipse pBca1020-r0040 carrying mRFP1 (obtained from BioBrick), pLambda driving Ti-U) equipped with a custom environmental chamber (50% humidity and 30 °C). mRFP1 (gift from Dr. Chenli Liu) and plasmids carrying YFP under constitutive A custom MATLAB script was used to control the microscope and its automated promoter were transformed into RP437 and into the inducible CheZ strain by stage (Prior) via the MicroManager interface . Time-lapse images (phase‐contrast electroporation. and fluorescence: mRFP1 or YFP) of the migrating cells were acquired using a Hamamatsu ORCA-Flash4.0 V2 camera (2,048 × 2,048 array of 6.5 × 6.5 μm pixels), a × 10 phase contrast objective (Nikon CFI Plan Fluor, N.A. 0.30, W.D. Microfluidic device design and fabrication. Microfluidic devices were con- 16.0 mm) and a LED illuminator (Lumencor SOLA light engine, Beaverton, OR) structed from the biocompatible and oxygen-permeable silicone polymer PDMS on through the mCherry block (Chroma 49008, Ex: ET560/× 40, Em: ET630/75 m) or cover glass following standard soft lithography protocols for two-layer devices . EYFP block (Chroma 49003; Ex: ET500/× 20, Em: ET535/30 m). Once the band The master molds for the device consisted of two silicon wafers with features formed, starting at the origin (closed end of the channel), the motorized stage created using ultraviolet (UV) photoresist lithography. The bottom wafer had two moved along the channel and paused every ~ 1.3 mm (the width of one frame with main parts: a large chamber created using SU-8-negative resist (thickness: 10 µm, a small overlap < 0.1 mm between consecutive positions) to take images in phase SU8 3010, Microchem) and a long channel together with two inlet/outlet channels contrast and fluorescence (exposure time 122 ms for both channels). After reaching designed to be opened and closed using pressure actuated valves. A second coat of the observation chamber, acquisition started over at the origin every 40 s (every 60 SPR-positive resist (thickness: 14 µm, SPR 220-7.0, MEGAPOSIT) on the same s for Fig. 1c). wafer was used to create a rounded channel profile that can collapse fully if In Fig. 1, all cells were expressing mRFP1. In Fig. 2, unlabeled and fluorescently depressed from above (Supplementary Fig. 1a). The second, top wafer contained labeled cells were separately grown to mid-exponential phase (OD600 of about 0.3). features for the control channels that close the collapsible features in the bottom For the experiments with 50, 100, and 200 µM Asp, the fluorescently labeled cells wafer. The top wafer was created using SU-8 negative resist (thickness: 10 µm, SU8 were diluted with unlabeled cells at the following ratios 1:20, 1:50, 1:100. The mixed 3010, Microchem). The resists were then cured using UV light exposure through cells were then washed and resuspended to the same predetermined density photomasks designed in CAD software and printed by CAD/Art Services, Inc. (OD600 of 0.7). The mixed populations were loaded into the microfluidic device (Bandon, Oregon), again following photoresist manufacturer specifications. and imaged as described above. In Fig. 4b-e, a similar procedure was followed to Subsequently, wafers were baked and the uncured photoresist was dissolved. After prepare the samples for different induction conditions. A 1:1 mixture of high and curing the SPR coat, the features were baked further to produce a rounded profile. low induction cultures was mixed and loaded in the device. After both wafers were complete, a protective coat of silane was applied by vapor In Fig.4f, a 1:1 mixture of high and low aTc induction (ranging from 1 to 10 ng −1 deposition. mL ) cultures with all cells fluorescently labeled by mRFP1 or YFP was mixed and To cast and manufacture the two-layer device, the top wafer was coated with a loaded in the device. The distances between the peaks of the density profiles of the 5 mm-thick layer of degassed 10:1 PDMS-to-curing agent ratio (Sylgard 184, Dow two populations were calculated by measuring the distance between the peaks in Corning). For the bottom layer, a 20:1 mixture was prepared and spin coated to the two fluorescent intensity profiles. The mean TB of each population were create a 100 μm-thick layer. The two layers were partially baked for 45 min at 70 °C. measured by loading a sample of the population on a cover slip and tracking The top layer was then cut and separated from the wafer, holes were punched from individual cells as previously described . the feature side using a sharpened 20-gauge blunt-tip needle to make external Once a band arrived in the perfused chamber, the gate near the chamber was connections to the control valve lines, then aligned and laminated onto the bottom closed (using 10 psi pressure) to capture the band (Supplementary Fig. 1b). To layer. The stacked layers were baked together for 1.5 h at 70 °C and allowed to cool. capture the second band in a separate experiment, the gate remained open until the NATURE COMMUNICATIONS (2018) 9:2177 DOI: 10.1038/s41467-018-04539-4 www.nature.com/naturecommunications 9 | | | ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04539-4 first band passed, then closed. The first band was immediately flushed away by fluorescence intensity, the opposite effect of decreased oxygen. Experiments were flowing M9 glycerol medium at 3 psi for 5 ~ 20 s through the chamber until all cells performed in 100 μm-deep and 14 μm-deep devices. The 100 μm-deep devices were were flushed from the chamber. The gate was then reopened to let the second band straight channels, 39 mm long and 600 μm wide. The 14 μm-deep devices were the migrate into the chamber. Once a band was captured in the perfused chamber, same as those used in the other experiments in this paper. However, as the control several pulses of fresh medium/buffer at 3 psi were flown in to reduce cell density gates interfere with the fluorescence signal, the top layer of PDMS was fabricated and homogenize the environment in the chamber. Cells were left to adapt in the without a mold for the gates. In both cases, a 1:500 dilution of concentrated perfused chamber for a few minutes. micelles was added to bacteria prepared as described above, just before loading into the microfluidic device. This dilution was intended to avoid sequestration of oxygen from the cells. Two control experiments were performed in the 14 μm TB detection. Once cell density became relatively homogeneous, the swimming device, one with dye and no cells, and one with cells and no dye. Imaging and data trajectories of individual cells were recorded for 2 min at 8 fps in order to extract analysis of all three types of experiments were performed in the same way. TB. We verified that perfusion of the observation chamber did not affect the Imaging equipment were the same as described above with a few exceptions. distribution of TBs (Supplementary Fig. 3). Each TB distribution was generated by The excitation filter from an ECFP block (Chroma 31044v2, D436/× 20), the acquiring four movies, tracking the individual cells, and determining their TB as emission filter from an mCherry block (Chroma 49008, ET630/75 m), and the 21,22 described before . The resulting number of sample trajectories longer than 10 s dichroic mirror from the mCherry block (Chroma 49008, T585lpxr) were used to (shorter tracks were discarded, because they provide poor estimate of TB ) was image the dye due to its large Stokes shift. limited by the number of independent movies one can acquire in the PDMS For the 100-μm device, the × 10 objective described above was used for imaging. chamber, the size of which corresponds to four field of views, and the density of Images were taken in 2 min intervals. For the 14-μm device, imaging was cells, which must be kept low for tracking to be possible. This procedure resulted in performed using a × 40 oil objective (Nikon CFI Plan Fluor, NA 1.3, W.D. 0.24 a minimum of 2823 and 3209 trajectories per distribution giving an error of at mm). Images were taken in 3.5 min intervals. most 1.15% and 1.00% in the determination of the cumulative distribution function To analyze the fluorescence images, fluorescence intensity was first averaged (CDF) estimated by bootstrapping for Figs. 1 and 2, respectively. over the width of the channel that was visible in the image, I(x, t). At each location x, the passing wave appeared as a brief peak in intensity during the time course. To Determination of the number of cells in the band. Image analysis was conducted separate this from slow variations in signal due to photobleaching and possible in MATLAB. We detected the position of the centroid of each fluorescent cell using global changes in oxygen concentration, we smoothed the time course of the MATLAB function bwconncomp. Figures 2a, 4b, and Supplementary Figs. 6d, e fluorescence intensity at each position, I (t), using MATLAB’s smooth function report these coordinates. The number of labeled cells was multiplied with the (smoothing method: lowess; window size: 4, corresponding to a time window of 14 dilution ratio to obtain the total number of cells in the band. min) to produce the slowly varying background signal, I ðtÞ. To extract the fast- Cell density profiles in Fig. 2d, 4c d, and Supplementary Figs. 9b, c were passing wave, we divided I (t) element wise by I ðtÞ for each position x. As a result, x x measured as follows: cell density profiles were extracted for each time point within the slowly changing background was normalized to 1, while faster changes in signal were different from 1. This also eliminated differences in illumination over space. one experiment and aligned before averaging. The cell density profile at a given time point was calculated by dividing the number of cells by the volume in one Finally, the noisy normalized profiles produced by this analysis were median filtered over space using MATLAB’s medfilt1 (window size: 10). spatial bin (~ 120 µm) along the observation channel. To reduce alignment error before averaging, the position of the peak density in each profile was identified by first smoothing the profile with a moving filter with a 5-bin span (MATLAB Oxygen consumption rate in batch cultures. RP437 cells were grown in 200 mL function smooth) and then identifying the position of the peak. To avoid boundary M9 glycerol medium to mid-exponential phase and washed twice with M9 glycerol effects, only the profiles with peak position located between 3 and 8 mm from the buffer supplemented with 200 µM Asp. Cells were then resuspended in 50 mL of origin were used to calculate the average density profile. The mean and SD shown the same defined buffer to an OD600 of 0.5. The sample was placed in a beaker at in the figures were calculated with raw cell density profiles (not smoothed). 30 °C. The surface of the buffer was sealed by overlaying mineral oil. The level of dissolved oxygen was measured with a portable dissolved oxygen meter (Mel- Measurement of amino acids by HPLC. When consuming Asp, E. coli cells waukee MW600) every 30 s. The cell sample was continuously stirred at 300 r.p.m. secrete other amino acids, which could affect group migration . We used HPLC to with a magnetic stirrer. The consumption rate per cell per minute was obtained by analyze the amino acids secreted by the cells when they are suspended in M9 dividing the reduction in oxygen by the number of cells and the sampling interval glycerol buffer supplemented with Asp. RP437 cells were grown in 200 mL M9 time (Supplementary Fig. 7g). glycerol medium up to mid-exponential phase and washed twice with M9 glycerol buffer supplemented with 500 µM Asp. Cells were then resuspended in 5 mL of the same defined buffer at an OD600 of 1 and placed in a 200 mL flask. The flask was Asp consumption rate in batch cultures. RP437 cells were grown in 200 mL M9 shaken at 200 r.p.m. to maximize aeration at 30 °C. Every 15 min, a 500 µl of glycerol medium up to mid-exponential phase and washed twice with M9 glycerol culture was sampled, filtered using 0.2 µm filter (Acrodisc 13 mm Syringe Filter buffer supplemented with 500 µM Asp. Cells were resuspended in 3 mL of the same with 0.2 µm HT Tuffryn Membrane, Pall Corporation), and analyzed by HPLC via defined buffer to make cell density at an OD600 of 1. Two samples were prepared pre-column derivatization method. The resulting derivatives were separated by in two test tubes. One tube was shaken at 200 r.p.m. to maximize aeration, whereas phase chromatography using a Dionex Ultimate 3000 HPLC, with a coupled DAD- the other tube was left on the bench with mineral oil overlaid on the liquid surface 3000RS diode array detector (Dionex) and FLD detector (Dionex) using an ACE to avoid supply of oxygen from air. These two tubes were incubated at 30 °C and C18 column (3 µm, 3 × 150 mm). Amino acid standard (AAS18 Sigma) was used as sampled every 10 min. The amount of Asp in the collected sample was then reference. measured by using Aspartate Assay Kit (Abcam). The consumption rate per cell Upon uptake of Asp, cells secreted small amounts of glutamate (Glu), per minute was obtained by dividing the reduction of Asp in the sample by the asparagine (Asn), and homoserine (HS), which are attractants (Supplementary number of cells and the sampling interval time (Supplementary Fig. 7g). Fig. 6c). To quantify the relative contribution of each amino acid to the chemotactic response in our experiment the measured concentration of each amino Statistical analysis and experimental reproducibility. No statistical methods acid was plotted in units of the corresponding EC50 of the dose response of the were used to predetermine sample size. Standard error in the CDF of the TB in chemotaxis system. The EC50 of the dose response of the chemotaxis system for each replicate experiment was determined by bootstrapping (1000 bootstrap each amino acid has been quantified by the Sourjik lab using in vivo Förster samples). One-sided unpaired two-sample Student’s t-test assuming unequal var- resonance energy transfer measurements in RP437, the strain used in this study. iances was used for comparison between two groups in Fig. 2b, c. P-values < 0.05 The EC50 values of RP437 E. coli are 0.3 μM for Asp, 50 μM for Glu, 30 μM for were considered statistically significant and marked with asterisks. For t-test, t- Asn, and 3 mM for HS . Supplementary Fig. 6c reveals that the response to Asp dominates by almost two orders of magnitude over the responses to the other values and degrees of freedom are provided in the figure legends. The error bars are defined in each figure caption and are standard deviation except in Fig. 2f. Data amino acids. It also shows that of the three secreted amino acids, Glu is the one that has the second largest effect, albeit still much smaller than the response to Asp. presented in the main figures were drawn from at least three independent repli- cates, with the exception of Fig. 4d(n = 2). The number of replicates is mentioned Finally, we also checked that chemotaxis towards oxygen does not play a significant role either. Mutant strains lacking the oxygen receptor aer or both aer and tsr form in the caption of each figure. bands under the same condition as in Fig. 2a, indicating that aerotaxis is not essential and Tar response to Asp is sufficient for the band to travel Code availability. To analyze the swimming behavior of E. coli cells we used (Supplementary Fig. 6d, e). 21,22 custom MATLAB code as reported in refs. . The code to simulate the math- ematical model is described above and available from the corresponding author Measurement of oxygen in the center of the traveling band. Ruthenium upon request. complexes are toxic to E. coli; hence, the need to encapsulate them in phospholipid micelles. This was achieved using the same protocol as in ref. . Fluorescence of the ruthenium complex is quenched by oxygen binding, so higher fluorescence Data availability. Data for each figure is provided as a MATLAB .fig file from corresponds to lower oxygen concentration. It is noteworthy that the high density which the data points can be extracted https://doi.org/10.6084/m9. of cells in the band could exclude the dye; however, this should decrease figshare.6207371. 10 NATURE COMMUNICATIONS (2018) 9:2177 DOI: 10.1038/s41467-018-04539-4 www.nature.com/naturecommunications | | | NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04539-4 ARTICLE Received: 25 July 2017 Accepted: 3 May 2018 35. Drescher, K., Dunkel, J., Cisneros, L. H., Ganguly, S. & Goldstein, R. E. Fluid dynamics and noise in bacterial cell-cell and cell-surface scattering. Proc. Natl Acad. Sci. USA 108, 10940–10945 (2011). 36. Zhang, H. P., Be’er, A., Florin, E. L. & Swinney, H. L. Collective motion and density fluctuations in bacterial colonies. Proc. Natl Acad. Sci. USA 107, 13626–13630 (2010). References 37. Brenner, M. P., Levitov, L. S. & Budrene, E. O. Physical mechanisms for 1. Sumpter D. J. T. Collective Animal Behavior. Princeton Univ. Press (2010). chemotactic pattern formation by bacteria. Biophys. J. 74, 1677–1693 (1998). 2. Liu, J. et al. Metabolic co-dependence gives rise to collective oscillations within 38. Holz, M. & Chen, S. H. Quasi-elastic light scattering from migrating biofilms. Nature 523, 550–554 (2015). chemotactic bands of Escherichia coli. Biophys. J. 23,15–31 (1978). 3. Couzin, I. D. & Krause, J. Self-organization and collective behavior in 39. Scribner, T. L., Segel, L. A. & Rogers, E. H. A numerical study of the formation vertebrates. Adv. Stud. Behav. 32,1–75 (2003). and propagation of traveling bands of chemotactic bacteria. J. Theor. Biol. 46, 4. Vicsek, T. & Zafeiris, A. Collective motion. Phys. Rep. 517,71–140 (2012). 189–219 (1974). 5. Potts, W. K. The chorus-line hypothesis of maneuver coordination in avian 40. Novick-Cohen, A. & Segel, L. A. A gradually slowing travelling band of flocks. Nature 309, 344–345 (1984). chemotactic bacteria. J. Math. Biol. 19, 125–132 (1984). 6. Ballerini, M. et al. Interaction ruling animal collective behavior depends on 41. Wong-Ng, J., Celani, A. & Vergassola, M. Exploring the function of bacterial topological rather than metric distance: evidence from a field study. Proc. Natl chemotaxis. Curr. Opin. Microbiol. 45,16–21 (2018). Acad. Sci. USA 105, 1232–1237 (2008). 42. Lapidus, I. R. & Schiller, R. A model for traveling bands of chemotactic 7. Beyerinck, M. W. Ueber atmungsfiguren beweglicher bakterien. Zentr bacteria. Biophys. J. 22,1–13 (1978). Bakteriol. Parasite. 114, 827 (1893). 43. Lauffenburger, D., Kennedy, C. R. & Aris, R. Traveling bands of chemotactic 8. Adler, J. Chemotaxis in bacteria. Science 153, 708–716 (1966). bacteria in the context of population-growth. B Math. Biol. 46,19–40 (1984). 9. Adler, J. Effect of amino acids and oxygen on chemotaxis in Escherichia coli. J. 44. Wang, Z. A. Mathematics of traveling waves in chemotaxis review paper. Bacteriol. 92, 121–129 (1966). Discret. Cont. Dyn.-B 18, 601–641 (2013). 10. Keller, E. F. & Segel, L. A. Traveling bands of chemotactic bacteria: a 45. Erban, R. & Othmer, H. G. From individual to collective behavior in bacterial theoretical analysis. J. Theor. Biol. 30, 235–248 (1971). chemotaxis. SIAM J. Appl. Math. 65, 361–391 (2005). 11. Saragosti, J. et al. Directional persistence of chemotactic bacteria in a traveling 46. Tweedy, L., Knecht, D. A., Mackay, G. M. & Insall, R. H. Self-generated concentration wave. Proc. Natl Acad. Sci. USA 108, 16235–16240 (2011). chemoattractant gradients: attractant depletion extends the range and 12. Wolfe, A. J. & Berg, H. C. Migration of bacteria in semisolid agar. Proc. Natl robustness of chemotaxis. PLoS Biol. 14, e1002404 (2016). Acad. Sci. USA 86, 6973–6977 (1989). 47. Hein, A. M. et al. The evolution of distributed sensing and collective 13. Spudich, J. L. & Koshland, D. E. Jr. Non-genetic individuality: chance in the computation in animal populations. eLife 4, e10955 (2015). single cell. Nature 262, 467–471 (1976). 48. Yang, Y. et al. Relation between chemotaxis and consumption of amino acids 14. Mayor, R. & Etienne-Manneville, S. The front and rear of collective cell in bacteria. Mol. Microbiol. 96, 1272–1282 (2015). migration. Nat. Rev. Mol. Cell Biol. 17,97–109 (2016). 49. Ford,R.M.&Cummings, P. T. On therelationshipbetween cell balanceequations 15. Kussell, E. & Leibler, S. Phenotypic diversity, population growth, and for chemotactic cell populations. SIAM J. Appl. Math. 52, 1426–1441 (1992). information in fluctuating environments. Science 309, 2075–2078 (2005). 50. Berg, H. C. & Brown, D. A. Chemotaxis in Escherichia coli analysed by three- 16. Veening, J. W., Smits, W. K. & Kuipers, O. P. Bistability, epigenetics, and bet- dimensional tracking. Nature 239, 500–504 (1972). hedging in bacteria. Annu. Rev. Microbiol. 62, 193–210 (2008). 51. Demir, M. & Salman, H. Bacterial thermotaxis by speed modulation. Biophys. 17. Frankel, N. W. et al. Adaptability of non-genetic diversity in bacterial J. 103, 1683–1690 (2012). chemotaxis. eLife 3, e03526 (2014). 52. Xia, Y. W. G. Soft lithography. Annu. Rev. Mater. Sci. 28, 153–184 (1998). 18. Ackermann, M. A functional perspective on phenotypic heterogeneity in 53. Edelstein, A. D. et al. Advanced methods of microscope control using microorganisms. Nat. Rev. Microbiol. 13, 497–508 (2015). muManager software. J. Biol. Methods 1, e10 (2014). 19. Copenhagen, K., Quint, D. A. & Gopinathan, A. Self-organized sorting limits 54. Adler, M., Erickstad, M., Gutierrez, E. & Groisman, A. Studies of bacterial behavioral variability in swarms. Sci. Rep. 6, 31808 (2016). aerotaxis in a microfluidic device. Lab. Chip. 12, 4835–4847 (2012). 20. Waite, A. J., Frankel, N. W. & Emonet, T. Behavioral variability and 55. Hazel, J. R. & Sidell, B. D. A method for the determination of diffusion phenotypic diversity in bacterial chemotaxis. Annu. Rev. Biophys. 47, coefficients for small molecules in aqueous solution. Anal. Biochem. 166, 27.21–27.22 (2018). 335–341 (1987). 21. Dufour, Y. S., Gillet, S., Frankel, N. W., Weibel, D. B. & Emonet, T. Direct 56. Schellenberg, G. D. & Furlong, C. E. Resolution of the multiplicity of the correlation between motile behavior and protein abundance in single cells. glutamate and aspartate transport systems of Escherichia coli. J. Biol. Chem. PLoS Comput. Biol. 12, e1005041 (2016). 252, 9055–9064 (1977). 22. Waite, A. J. et al. Non-genetic diversity modulates population performance. 57. Wetzel R. G. in Limnology: Lake and River Ecosystems, 3rd edn (Academic Mol. Syst. Biol. 12, 895 (2016). Press, 2001). 23. Ben-Jacob, E., Cohen, I. & Levine, H. Cooperative self-organization of 58. Vollmer, A. P., Probstein, R. F., Gilbert, R. & Thorsen, T. Development of an microorganisms. Adv. Phys. 49, 395–554 (2000). integrated microfluidic platform for dynamic oxygen sensing and delivery in a 24. Tai-Chia, L. & Zhi-An, W. Development of traveling waves in an interacting flowing medium. Lab. Chip. 5, 1059–1066 (2005). two-species chemotaxis model. Discret. Contin. Dyn. Syst. A 34, 2907–2927 (2014). 25. Emako, C., Gayrard, C., Buguin, A., Neves de Almeida, L. & Vauchelet, N. Acknowledgements Traveling pulses for a two-species chemotaxis model. PLoS Comput. Biol. 12, We thank H. Salman for sharing the protocol to synthesize the Ru-micelles; Y. Dufour e1004843 (2016). for help with the microfluidics; F. Isaacs for sharing BioBrick plasmids; C. Liu, S. Par- 26. Mittal, N., Budrene, E. O., Brenner, M. P. & Van Oudenaarden, A. Motility of kinson, and H. Salman for sharing E. coli strains and plasmids; N. Clay, B. Barco, and W. Escherichia coli cells in clusters formed by chemotactic aggregation. Proc. Natl Chezem for help with HPLC experiments; D. Clark and J. Wang for help with statistics;, Acad. Sci. USA 100, 13259–13263 (2003). and A. Waite, J. Howard, C. Brannon, A. Greene, and P. Turner for comments on the 27. Celani, A. & Vergassola, M. Bacterial strategies for chemotaxis response. Proc. manuscript. This study was supported by the National Institutes of Health grant Natl Acad. Sci. USA 107, 1391–1396 (2010). 1R01GM106189 and the Allen Distinguished Investigator Program (grant 11562) 28. Si, G., Wu, T., Ouyang, Q. & Tu, Y. Pathway-based mean-field model for through The Paul G. Allen Frontiers Group. C.H. acknowledges support through the Escherichia coli chemotaxis. Phys. Rev. Lett. 109, 048101 (2012). National Natural Science Foundation of China (31570095, 11504399). 29. Dufour, Y. S., Fu, X., Hernandez-Nunez, L. & Emonet, T. Limits of feedback control in bacterial chemotaxis. PLoS Comput. Biol. 10, e1003694 (2014). Author contributions 30. Sourjik, V. & Berg, H. C. Functional interactions between receptors in X.F. and S.K. contributed equally to this work. H.M. and J.L. contributed equally to this bacterial chemotaxis. Nature 428, 437–441 (2004). work. T.E., X.F., and S.K. designed the research. X.F., S.K., H.M., and C.H. performed the 31. Long, Z., Quaife, B., Salman, H. & Oltvai, Z. N. Cell-cell communication experiments. X.F., S.K., J.L., H.M., and T.E. performed the data analysis. X.F., J.L., H.M., enhances bacterial chemotaxis toward external attractants. Sci. Rep. 7, 12855 and T.E. performed the mathematical modeling. D.C.V. and S.W.Z. contributed to (2017). mathematical modeling. T.E., X.F., S.K., H.M., and J.L. wrote the manuscript. All authors 32. Mesibov, R. & Adler, J. Chemotaxis toward amino acids in Escherichia coli. discussed the manuscript. J. Bacteriol. 112, 315–326 (1972). 33. Douarche, C., Buguin, A., Salman, H. & Libchaber, A. E. coli and oxygen: a motility transition. Phys. Rev. Lett. 102, 198101 (2009). Additional information 34. Wu, X. L. & Libchaber, A. Particle diffusion in a quasi-two-dimensional Supplementary Information accompanies this paper at https://doi.org/10.1038/s41467- bacterial bath. Phys. Rev. Lett. 84, 3017–3020 (2000). 018-04539-4. NATURE COMMUNICATIONS (2018) 9:2177 DOI: 10.1038/s41467-018-04539-4 www.nature.com/naturecommunications 11 | | | ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04539-4 Competing interests: The authors declare no competing interests. appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party Reprints and permission information is available online at http://npg.nature.com/ material in this article are included in the article’s Creative Commons license, unless reprintsandpermissions/ indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in regulation or exceeds the permitted use, you will need to obtain permission directly from published maps and institutional affiliations. the copyright holder. To view a copy of this license, visit http://creativecommons.org/ licenses/by/4.0/. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give © The Author(s) 2018 12 NATURE COMMUNICATIONS (2018) 9:2177 DOI: 10.1038/s41467-018-04539-4 www.nature.com/naturecommunications | | | http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Nature Communications Springer Journals

Spatial self-organization resolves conflicts between individuality and collective migration

Free
12 pages
Loading next page...
 
/lp/springer_journal/spatial-self-organization-resolves-conflicts-between-individuality-and-x076QqU7Vo
Publisher
Nature Publishing Group UK
Copyright
Copyright © 2018 by The Author(s)
Subject
Science, Humanities and Social Sciences, multidisciplinary; Science, Humanities and Social Sciences, multidisciplinary; Science, multidisciplinary
eISSN
2041-1723
D.O.I.
10.1038/s41467-018-04539-4
Publisher site
See Article on Publisher Site

Abstract

ARTICLE DOI: 10.1038/s41467-018-04539-4 OPEN Spatial self-organization resolves conflicts between individuality and collective migration 1,6 1,7 1,2 1 3 1,8 4,5 1,2 X. Fu , S. Kato , J. Long , H.H. Mattingly ,C.He , D.C. Vural , S.W. Zucker & T. Emonet Collective behavior can spontaneously emerge when individuals follow common rules of interaction. However, the behavior of each individual differs due to existing genetic and non- genetic variation within the population. It remains unclear how this individuality is managed to achieve collective behavior. We quantify individuality in bands of clonal Escherichia coli cells that migrate collectively along a channel by following a self-generated gradient of attractant. We discover that despite substantial differences in individual chemotactic abilities, the cells are able to migrate as a coherent group by spontaneously sorting themselves within the moving band. This sorting mechanism ensures that differences between individual chemo- tactic abilities are compensated by differences in the local steepness of the traveling gradient each individual must navigate, and determines the minimum performance required to travel with the band. By resolving conflicts between individuality and collective migration, this mechanism enables populations to maintain advantageous diversity while on the move. 1 2 Department of Molecular, Cellular, and Developmental Biology, Yale University, New Haven, CT 06520, USA. Department of Physics, Yale University, New Haven, CT 06520, USA. Institute of Synthetic Biology, Shenzhen Institutes of Advanced Technology, Chinese Academy of Sciences, Shenzhen 518055, 4 5 China. Department of Computer Science, Yale University, New Haven, CT 06520, USA. Department of Biomedical Engineering, Yale University, New Haven, CT 06520, USA. Present address: Institute of Synthetic Biology, Shenzhen Institutes of Advanced Technology, Chinese Academy of Sciences, Shenzhen 518055, China. Present address: Department of Molecular Biotechnology, Graduate School of Advanced Sciences of Matter, Hiroshima University, Higashi-Hiroshima, Hiroshima 739-8530, Japan. Present address: Department of Physics, University of Notre Dame, Notre Dame, IN 46556, USA. These authors contributed equally: X. Fu and S. Kato. Correspondence and requests for materials should be addressed to T.E. (email: thierry.emonet@yale.edu) NATURE COMMUNICATIONS (2018) 9:2177 DOI: 10.1038/s41467-018-04539-4 www.nature.com/naturecommunications 1 | | | 1234567890():,; ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04539-4 ells and larger organisms exhibit collective behaviors that quantify the distribution of phenotypes in the band (Fig. 1b and 1,2 5 are often advantageous to the participating individuals . Supplementary Fig. 1). Approximately 2 × 10 clonal E. coli cells CMany such collective behaviors dynamically emerge when grown in M9 glycerol medium (M9 salts, glycerol, and casamino a large number of individuals follow the same rules to interact acids; Methods) were introduced with fresh medium into the 3,4 with each other and the environment . Prominent examples are device and concentrated at the end of the channel by cen- 4–6 bird flocks and the collective migration of bacteria along trifugation (Methods). Following centrifugation, sequential bands 7–11 8,9,12 channels and on agar plates . At the same time, pheno- of cells collectively migrated along the channel at different but typic differences among even genetically identical individuals are nearly constant speeds (Fig. 1c), presumably consuming different a ubiquitous feature of biology . Phenotypic diversity can lead to compounds within the undefined media, as demonstrated in early useful leader–follower structures within a traveling group. For studies . example, in migrating neural crest cells and in fish shoals, many E. coli cells navigate by alternating straight “runs” with organisms may follow a few more informed individuals .In “tumbles” that randomly reorient their swimming direction microbial communities, maintaining diversity in the population (Fig. 1d). By transiently suppressing tumbles whenever attractant can enable bet-hedging strategies to survive uncertain environ- signal increases, they perform a biased random walk that allows 13,15–18 20 ments and resolve trade-offs . However, heterogeneity can them to move toward higher concentrations of attractant . In the also be disruptive, as is the case in simulated swarms where non- absence of a gradient, the fraction of time a cell spends tumbling aligners tend to be purged from the swarm . This raises a —its tumble bias (TB)—remains approximately constant and dilemma: although phenotypic diversity provides advantages, it therefore can be used as a quantitative measure of the phenotype also tends to reduce coordination. of the cell. Importantly, using the same strain and microfluidic One of the simplest cases of collective behavior is exhibited by channel depth, we previously demonstrated that the TB is a bacteria: clonal populations of motile Escherichia coli cells col- strong determinant of chemotactic performance in liquid: lower lectively migrate when placed at high density at the bottom of a TB cells drift significantly faster up a static gradient than higher 8–11 22 tube filled with nutrients . This collective behavior is mediated TB cells . To quantify the distribution of phenotypes in the by the well-characterized chemotaxis system , which enables the isogenic population that was introduced in the device, a low bacteria to follow chemical gradients, in this case generated by density of cells was loaded into the microfluidic device without their consumption of attractant present in the medium (Fig. 1a). centrifugation and individual cells were tracked to determine 21,22 However, populations of E. coli exhibit substantial cell-to-cell their TB, as previously described . TB was broadly distributed 13,21 variability in their swimming phenotypes and hence che- in the population with some cells tumbling < 10% of the time (i.e., motactic abilities, even when all cells are genetically identical . TB < 0.1) and others > 50% of the time (Fig. 1e black), consistent 13,21,22 How bacterial populations manage phenotypic heterogeneity to with previous studies . Given the functional consequences still allow coordinated collective migration remains largely of this non-genetic diversity, how can the same population of cells unknown, mainly because of the difficulties in measuring cellular migrate together as a coordinated group, as shown in Fig. 1c? behavior at both the collective and the individual levels in the To answer this question, we first considered whether all same experiment . phenotypes or only a subset of them traveled in each band. We Although the migration of traveling waves or “bands” of bac- used pressure valves to capture one band of cells at a time in the teria has also served as a classic model for the theoretical study of wide chamber of our device (Fig. 1b and Supplementary Fig. 1b). 10,23 emergent phenomena and pattern formation in biology , the After trapping cells in the wide chamber, it was perfused with effect of non-genetic diversity on this process has scarcely been fresh media to homogenize the environment and dilute the cell examined. Previous studies examined how two populations may density. We verified that perfusion of the wide chamber did not 24,25 travel together ; however, it was assumed that within each affect the distribution of TBs (Supplementary Fig. 3). Dilution population all of the individuals were identical. The mechanisms enabled us to track individual cells. Homogenization ensured that by which a continuum of phenotypes can achieve coherent cells had adapted back to a uniform environment and were not migration have not been investigated. responding to an attractant gradient when we measured their Here we used a microfluidic system that enables precise TBs. The distribution of TB was shifted toward lower TB in both quantitative measurements at the individual and collective scales traveling bands compared with the original distribution (Fig. 1e), to study the interplay of diversity and collective bacterial suggesting that it was more difficult for high TB cells to migration. Our central finding is that within the traveling band, participate in collective migration. Selection against high TB cells cells spontaneously sort themselves such that their chemotactic was stronger in the faster band (Fig. 1e, red) than in the slower abilities are matched to the local gradient steepness, enabling one (Fig. 1e, yellow). Cell density and number also varied between diverse cells to travel together with the same drift speed. the two bands (Fig. 1c), suggesting that there were interdepen- Extending the classic Keller–Segel model of traveling bands to dencies between the speed of the group, its size, and the diversity account for diversity predicts this spatial sorting and qualitatively of the individuals able to migrate with the group. We periodically recapitulates the experimental results. Our second finding is a tracked cells after they were trapped and diluted in the wide novel mechanism that reduces the rate at which cells fall off the chamber and found that the original TB distribution was recov- back of the band: when attractant consumption depends on local ered after growth (Fig. 1f). Thus, selection of low TB cells by the oxygen, oxygen limitation in the center of the band increases the collective migration was not due to genetic heterogeneity. In gradient of attractant at the back, helping cells there keep up. addition, it is unlikely that cell growth affected the TB selection Together, these two mechanisms enable populations of bacteria to while cells were traveling in a band, because the duration of the maintain diversity while migrating as a group. experiment (30 min) was shorter than the cell doubling time (~ 55 min, Supplementary Fig. 4). Results Collective migration selects against high TB cells. To determine Cells of diverse chemotactic abilities migrate as a group.To the relationship between the number of cells in the band, the quantify collective behavior and diversity in the same experiment, band speed, and diversity, we switched from casamino acids to a we designed a microfluidic device consisting of a long channel to defined M9 glycerol buffer containing aspartate (Asp) as the observe the traveling band , followed by a large chamber to only limited chemoattractant (Methods). In this condition, a 2 NATURE COMMUNICATIONS (2018) 9:2177 DOI: 10.1038/s41467-018-04539-4 www.nature.com/naturecommunications | | | NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04539-4 ARTICLE ab Group migration Attractant Cell density Channel Chamber cd Individual diversity 720 Band 2 Band 1 Tumble bias 0 3.5 7 0.1 0.2 0.3 0.4 Position (mm) ef 01234 5 Time elapsed (hrs) 6 6 4 4 2 2 0 0 0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 0.5 Tumble bias Tumble bias Fig. 1 Collective migration of a phenotypically-diverse clonal population. a When concentrated at the bottom of a nutrient channel, motile E. coli cells emerge from the high cell density region and travel in bands along the channel by following gradients of attractant produced by their consumption. b Microfluidic device used to quantify the band migration and the phenotypic diversity within the band. Control gates along the channel (black vertical lines) are initially open (top) and later closed to capture different bands of cells in the observation chamber (bottom), where single cells are tracked to quantify the distribution of phenotypes within the band (Supplementary Fig. 1). c Time-lapse imaging of E. coli cells expressing the fluorescent protein mRFP1 showing the collective migration of bands in M9 glycerol medium (M9 salts, glycerol, and casamino acids; Methods). In this undefined medium, −1 −1 8 several bands emerge that travel at different speeds (red: 0.68 mm min , yellow: 0.23 mm min ) . We verified that labeling cells did not affect band speeds nor tumble bias distributions (Supplementary Fig. 2). Scale bar, 0.6 mm. d The tumble bias (color)—average probability to tumble—of individual 21, 22 cells was quantified by tracking a cell for 2 min in a uniform environment (no gradient) and detecting tumbles (black dots) as previously described . Scale bar, 200 μm. e Collective migration selects against high TB cells. Tumble bias distributions from the first (red) and second (yellow) bands (n = 3), and from the population that was introduced in the device (black; n = 3). f The tumble bias distribution of the cells in the first wave (red in e) gradually shifted back toward the original distribution (black) during growth in the perfused chamber. Fresh M9 glycerol medium was supplied every 30 min. The TB distribution was measured every 20 min single band formed (Fig. 2a) and its speed could be tuned by density of cells in the band increased (Fig. 2c, d), and the dis- changing the concentration of supplemented Asp (Fig. 2b). To tribution of TB within the band shifted toward higher TB measure band speed and density, cells expressing mRFP1 were (Fig. 2e). In general, collective migration selected against high TB mixed with unlabeled cells at ratios of 1:20, 1:50, or 1:100 (for 50, cells, with selection being stronger in faster bands (Fig. 2e, f). It is 100, and 200 μM Asp, respectively), and their positions were noteworthy, however, that diversity was not eliminated—all 11,26 detected at various time points (Fig. 2a and Methods ). Band bands still exhibited a range of TBs. Thus, although collective speed, number of cells in the band, and density profiles were migration selected against high TB cells, it was still possible for a stable over time with a slow decay due to cells falling off at the diverse group to travel together. back of the band (Fig. 2b, c, d and Supplementary Fig. 5). Although there were variations across experimental replicates due Extending the Keller–Segel model to account for diversity.To to variations in the number of cells introduced in the device, the better understand how collective migration selects TB, we relationship between Asp concentration and band parameters was extended the classic Keller–Segel mathematical model describing consistent. Specifically, as the concentration of Asp increased, the traveling bands of bacteria to include phenotypic diversity speed of the band decreased (Fig. 2b), the number and peak (Methods Eqs. 1–3). In this model, cells consume the diffusible NATURE COMMUNICATIONS (2018) 9:2177 DOI: 10.1038/s41467-018-04539-4 www.nature.com/naturecommunications 3 | | | Time (sec) PDF PDF ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04539-4 a d 9 ×10 0 3.5 7 –1.5 –1 –0.5 0 0.5 1 1.5 Position (mm) Relative coordinate, z (mm) be *p = 0.004 0.5 *p = 0.0076 0.4 0.3 n = 5 0.2 n = 4 0.1 n = 5 Simulation 0 0 50 100 150 200 0.1 0.2 0.3 0.4 0.5 Asp conc. (µM) Tumble bias c f *p = 0.0017 *p = 0.0034 1/2 10 Simulation 1/4 50 100 150 200 0.1 0.2 0.3 0.4 0.5 Asp conc. (µM) Tumble bias Fig. 2 Relationship between band speed, density, and phenotypic diversity. a Time-lapse coordinates of cells (black dots) traveling in M9 glycerol buffer (Methods) with 200 µM aspartate. Only one band forms in each experiment. Cells (1:100) were labeled with mRFP1 and their coordinates detected (Methods). Scale bar, 0.6 mm. In remaining panels, colors are aspartate concentration in the buffer: 200 µM (blue; 1:100 cells labeled), 100 µM (yellow; 1:50 cells labeled), 50 µM (red; 1:20 cells labeled) (Methods). b Band speed decreased with aspartate concentration in the buffer. Circles: experiments in which tumble bias was measured, used in e, f. Diamonds: tumble bias was not measured. Dashed: simulations (Methods and Supplementary Fig. 5). p-value: one-sided, two-sample t-test assuming unequal variances, t-value= 4.1, df = 5.4 between 50 μM and 100 μM, t-value= 4.2, df= 3.8 between 100 μM and 200 μM. c Number of cells traveling in the band increased with aspartate concentration. Circles, diamonds, and dashed: same as in b. p-value: one- sided, two-sample t-test assuming unequal variances, t-value= − 3.9, df = 6.5 between 50 μM and 100 μM, t-value=− 4.6, df = 6.4 between 100 μM and 200 μM. d Cell density profiles averaged over time and experiments (Supplementary Fig. 5). Line: average over n = 5 (red), n = 4 (yellow), n = 5 (blue) replicate experiments. Shade: SD. For each experiment, nine profiles were measured at 1.3 min intervals. Dashed: simulations. e Tumble bias distribution of the cells that traveled with the band. Lines: average over experiments; shading: SD; dashed lines: simulations. f Ratio between the distribution of tumble bias in the traveling band (Fig. 2e, colored lines) and that of the original population (Fig. 2e, black) quantifies the enrichment of tumble bias. Lower aspartate concentrations enrich more for lower tumble bias. Solid lines: mean of the measurements from experiments; shading: standard error of the mean; dashed: simulations attractant Asp, which generates a traveling gradient that the cells indeed, for TB = 0 the cell is just diffusing and χ = 0; 21,22 follow by biasing their random walk (Fig. 3a). Faster consump- Methods) . Thus, in a gradient of attractant, cells with higher tion, more cells, or less Asp all lead to a faster traveling gra- TB do not diffuse as much and climb slower than cells with lower dient . The motion of a phenotype i depends on two parameters: TB. The dependency of f on Asp concentration has been char- its effective diffusion coefficient, μ = μ(TB ), which results from acterized as well . Moreover, we conducted high-performance i i the cells’ random walk, and its chemotactic coefficient, χ = χ liquid chromtography (HPLC) experiments to verify that the (TB ), which quantifies how effective that phenotype is at biasing chemotaxis response to Asp dominated that to amino acids 27–29 its motion to follow the perceived amount of Asp, f. Theory secreted as byproducts of Asp metabolism. Experiments con- and tracking of individual E. coli cells swimming in a static gra- ducted with mutants lacking the oxygen receptor aer or both dient of α-methylaspartate (non-metabolizable analog of Asp) aer and tsr indicated that aerotaxis was not essential and Tar have shown that μ and χ are decreasing functions of the TB (note response to Asp was sufficient for band migration (Methods that for very low values of TB < ~ 0.05, χ increases with TB; and Supplementary Fig. 6). 4 NATURE COMMUNICATIONS (2018) 9:2177 DOI: 10.1038/s41467-018-04539-4 www.nature.com/naturecommunications | | | Migration speed –1 Time (min) Cell number in the band (mm min ) –1 Fold enrichment Cell density (ml ) PDF NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04539-4 ARTICLE ab Total TB 0.1~0.2 TB 0.3~0.4 TB 0.2~0.3 TB 0.4~0.5 ×10 df df C = 2 = 1 8 dz dz 6 20 Falling off Traveling –1.5 –1 –0.5 0 0.5 1 Relative coordinate, z (mm) cd –1 TB 0.1 TB 0.2 10 20 TB 0.3 TB 0.4 0.2 0.1 Individual Group –2 –1.5 –1 –0.5 0 0.5 1 –1.5 –1 –0.5 0 0.5 1 Relative coordinate, z (mm) Relative coordinate, z (mm) Fig. 3 Mathematical modeling predicts a mechanism for consistent collective migration of diverse phenotypes. a Collective migration of diverse phenotypes at the same speed is made possible by the spontaneous spatial ordering of individual phenotypes within the band such that each individual’s chemotactic df ability is matched to the local gradient steepness . The proportion of better performers (larger χ, lower TB; red) should be enriched where the gradient is dz shallower (front), whereas the proportion of weaker performers (smaller χ, higher TB; blue) should be enriched where the gradient signal is steeper (back). The position where the perceived gradient steepness is maximum (dashed border of the gray regions) determines the highest tumble bias able to travel with the band. Cells in the gray region slowly fall out of the band. b Simulated density profiles of cells migrating in 200 μM aspartate (the same simulation is shown in Fig. 2def, blue dashed) show sorting based on tumble bias. c Chemotactic coefficient χ(z) (blue) defined by the phenotype whose density df profile peaks at position z and perceived gradient steepness (black). Red symbols correspond to the location of the peak cell density for individual dz df phenotypes. d Spatial sorting enables consistent migration velocity for traveling phenotypes. The migration velocity, χðÞ z , of the phenotype whose density dz profile peaks at position z gradually decreases toward the back of the band until the gray region is reached. In the gray region the migration velocity falls off more rapidly, preventing the high TB phenotypes located there from staying in the band To metabolize Asp, E. coli consumes oxygen . Introducing a Spatial sorting as a mechanism for inclusive migration.An fluorescent oxygen sensor in the M9 glycerol buffer revealed important feature of the experiments reproduced by the simula- that oxygen availability is reduced in the center of the traveling tions is the increasing selection against high TB cells as the band where cell density is high (Supplementary Fig. 7a-f). This amount of Asp is reduced (Fig. 2f). What mechanism enables cells results in a dependency of the average consumption rate of Asp with diverse chemotactic abilities to collectively migrate together on cell density (Supplementary Fig. 7h). We modeled this effect as one band, and what controls the upper bound on the TB such that the Asp consumption rate depended linearly on oxygen among those able to migrate together? For every phenotype that concentration and constrained the related parameters by travels with the band at constant speed c, the flux of cells must be measuring oxygen and Asp consumption rates in batch cultures approximately invariant in time and equal to the chemotactic flux (Supplementary Fig. 7g). For simplicity, we ignored possible minus that due to diffusion. Focusing on the partial differential phenotypic diversity in the Asp consumption rate, as well as the equation for the cell density of phenotype i and switching to the possible dependence of the diffusion coefficient on cell density , moving reference frame z = x − ct (here x and t are absolute which was found previously to be negligible in similar position and time), two predictions emerge (Methods Eqs. 7–8). experiments . We also omitted possible contributions of First, the traveling phenotypes must distribute themselves spa- hydrodynamics and physical interactions between cells, which tially within the band so that differences in local signal strengths, df can become important when bacteria swarm over surfaces . the slope of the black line in Fig. 3a, compensate for differences dz The resulting model (Methods and Table 1) qualitatively in the chemotactic abilities of each phenotype, χ , such that reproduced the main features of our experiments, including the df c ¼ χ ðÞ z , where z is the position of peak density of phenotype i dz i dependency on Asp concentration of the band speed (Fig. 2b), cell i (Fig. 3a). Therefore, this spatial sorting places the better per- number and density (Fig. 2c, d), TB distribution (Fig. 2e), formers (higher χ , lower TB ) in the front of the band, where the i i phenotypic selection (Fig. 2f), and average Asp consumption rate gradient is shallower and more difficult to follow (Fig. 3a, red), per cell as a function of cell number in the band (Supplementary and the weaker performers (lower χ , higher TB ) at the back, i i Fig. 7h). NATURE COMMUNICATIONS (2018) 9:2177 DOI: 10.1038/s41467-018-04539-4 www.nature.com/naturecommunications 5 | | | Chemotactic coefficient 2 –1 (mm min ) Perceived signal gradient –1 df /dz (mm ) Migration speed –1 (mm min ) –1 Cell density (ml ) Perceived signal f (a.u.) ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04539-4 where the gradient is steeper (Fig. 3a, blue). Furthermore, a result of spatial sorting. Due to experimental limitations, we could second prediction is that the gradient will reach a maximal not measure the TB distributions of all four strains in each of the steepness (Fig. 3a dashed border of gray zone), determining the experiments reported in Fig. 4a-e. To better quantify the rela- weakest phenotype that can travel (lowest χ , upper bound on tionship between peak separation and difference in ⟨TB⟩ value, TB ). Thus, we see the interplay between individuals i and the we instead mixed pairs of populations using only fluorescent effect of the community on the available resource f. strains induced with different levels of anhydrotetracycline (aTc). Analysis of simulations confirmed these analytical predictions For each pair, we measured the TB distributions right before df loading cells in the device, and then measured the distance (Fig. 3b, c, d). The steepness of the perceived signal , which dz between the fluorescence peaks in the resulting traveling band emerges dynamically from the cells’ consumption, peaked at the (Fig. 4f). This confirmed that there is a monotonic relationship back of the traveling band (low z) and decayed toward the front between peak separation and difference in ⟨TB⟩ values. Therefore, (high z) (Fig. 3c black). In contrast, the position z of the peak cells of various phenotypes appear to spontaneously sort them- density of phenotype i increased with its chemotaxis coefficient χ selves along the traveling band according to their TB, enabling (Fig. 3c blue), revealing an ordering of the phenotypes within the them to migrate collectively despite phenotypic differences. migratory band (Fig. 3b, c). Multiplying the two together gave a nearly constant velocity throughout the band, thus providing an explanation of how the various phenotypes might travel together Discussion (Fig. 3d). The rightmost points in Fig. 3c (blue line) and Fig. 3d How do organisms maintain collective behavior despite the (black line) correspond to the phenotype with the maximum potential conflicts created by phenotypic diversity among indi- chemotactic coefficient in the band. Ahead of that location, there viduals? We studied this question using traveling bands of che- are no more peaks in the cell density of any phenotype; however, motactic E. coli, which collectively migrate at the same speed there are cells due to diffusion. At the back of the band, cells with despite differences in chemotactic abilities of individuals in the TB higher than the predicted upper bound rapidly fell off of the band. Our key result is that spontaneous spatial organization of band (Fig. 3a, d gray zones). Importantly, these predictions phenotypes within a traveling band helps resolve the conflicts emerged from just the dynamics of cell density (Methods Eq. 1) between phenotypic diversity and collective migration. By and therefore hold true irrespective of whether oxygen is included matching individual abilities to the local difficulty of the navi- in the model. gation task within the band, this sorting mechanism ensures Comparing simulations with and without the oxygen- consistent migration speed across the band. This process also dependent consumption rate revealed that the oxygen depen- determines the minimum chemotactic performance required to dency reduces leakage of the highest-TB cells located at the back keep up with the band, therefore explaining how diversity can of the band (Supplementary Fig. 8b). The higher concentration of become limited by collective behavior. Thus, the mechanism oxygen at the back relative to the center of the band locally reported here enables a continuum of phenotypes to migrate increases the rate of consumption, and hence the slope of the coherently. traveling gradient of Asp, helping the cells there stay longer with In the traveling band, there is always a slow leakage of cells off the band (Methods Eq. 9). Thus, oxygen dependency has a similar the back of the band because of the finite sensitivity of the cells for cohesive effect as the secretion of a self-attractant, which also 38–41 11,25,37 the attractant they are chasing . High TB cells, in particular, helps reduce the leakage of cells . which are localized at the back of the wave, are at risk of falling off. We discovered that this leakage can be reduced (but not Cells in the band are spatially sorted by TB. To experimentally eliminated) if the consumption rate of the attractant is lower in test the prediction of spatial sorting in the band by phenotype, we the center of the band than at the back, where the consumption measured the relative position of two populations of cells with rate determines the local gradient steepness and chemotactic drift. different mean TB within the traveling band. The distribution of In our case, this arises because Asp consumption depends on TB in the population was controlled by manipulating the level of oxygen, which becomes limited in the center of the band where expression of the phosphatase CheZ, which deactivates the che- cell density is high. This mechanism provides an alternative to motaxis response regulator CheY (Supplementary Fig. 9a). We other mechanisms known to reduce cell leakage, such as the 11,25,37 generated multiple populations with different TB distributions of secretion of an attractant by the traveling cells . Note that varying mean TB (⟨TB⟩) (Fig. 4a). In each population, we labeled the spontaneous sorting mechanism discussed above helps com- 1 in 50 cells of the same genetic background and induction level pensate for differences in chemotactic abilities, irrespective of the by ectopically expressing either mRFP1 or YFP. The inducer was presence of such auxiliary mechanisms (oxygen or self- washed away when cells were resuspended in buffer before attractant). starting the migration experiment. TB distributions are stable for Traveling bands of bacteria have been studied for decades since 21 8,9 more than an hour in buffer . When cells from the lowest (red) Julius Adler’s experiments in capillary tubes . Adler reported the and highest (cyan) TB distributions were mixed in equal parts formation of multiple traveling bands in complex media; we and introduced in the device, spontaneous spatial order emerged, observe the same in casamino acids (Fig. 1) and expect multiple with cells from the low (high) TB distribution located at the front bands to be able to form when multiple consumable attractants (back) of the traveling band (Fig. 4b, c and Supplementary Fig. 9b are present. Within a migrating band the cells respond to the for replicates). The distance between the peaks of the density traveling gradient, some parts of which can be fairly steep. profiles of the two populations remained nearly constant over the Therefore, we expect the instantaneous TB of an individual cell to duration of the experiment, indicating that the two populations be dynamically changing depending on its direction of motion 11,29 traveled together at the same speed (Fig. 4e). We verified that the and position within the gradient, as previously reported . Here distance between the peak densities of the two populations was we showed that phenotypic (intrinsic) differences in adapted TB stable in longer experiments (Supplementary Fig. 9d). Mixing between cells contribute significantly to spatial structure within populations with closer TB distributions caused the peaks of the the traveling band. In future studies, we will separate the con- two traveling populations to be closer (magenta and green in tributions of phenotypic and dynamic diversity to group struc- Fig. 4a, d and Supplementary Fig. 9c for replicates), suggesting ture. It will also be interesting to examine the contribution of that distance between peaks increases with difference in ⟨TB⟩ as a initial conditions, dimensionality, and growth (necessary to 6 NATURE COMMUNICATIONS (2018) 9:2177 DOI: 10.1038/s41467-018-04539-4 www.nature.com/naturecommunications | | | NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04539-4 ARTICLE a d 9 ×10 40 ng/ml aTc 6 ng/ml aTc 4 ng/ml aTc 2 ng/ml aTc –1.5 –1 –0.5 0 0.5 1 1.5 0 0.2 0.4 0.6 0.8 1 Tumble bias Relative coordinate, z (mm) b e 560 600 2 468 2 468 Position (mm) Peak position (mm) c f ×10 0.4 0.3 0.2 0.1 0 0 –1.5 –1 –0.5 0 0.5 1 1.5 0 0.1 0.2 0.3 0.4 Relative coordinate, z (mm) Difference in mean tumble bias Fig. 4 Phenotypes spontaneously order themselves along the traveling band according to tumble bias. a The distribution of tumble bias in the population can be controlled by manipulating the level of expression of the phosphatase CheZ, which deactivates the chemotactic response regulator CheY (Supplementary Fig. 9a). b Time-lapse coordinates (colored dots) of an equal mixture of low (red in a) and high (cyan in a) TB cells traveling in 200 µM aspartate M9 glycerol buffer (see Methods). Scale bar, 0.6 mm. c Corresponding density profiles (colors) together with total cell density (black). (Line: mean over n = 34 time points measured at 40 s intervals for one experiment; shading: SD; five replicates are in Supplementary Fig. 9b). d Same as in c, but for the magenta and green populations in a. Two replicates are in Supplementary Fig. 9c. e Peak positions as a function of time for the experiment in b. f The distance between the fluorescence intensity peaks of the two populations traveling together in a single band increases with the difference between the mean TB of the two populations. For each independent experiment, two populations labeled with different fluorescent proteins (mRFP1 or YFP) were induced using different aTc levels to obtain distributions with different mean tumble biases. Dots: average over n = 4 experiments; error bars are SD 41–43 maintain traveling cell density over long times ) to this the range of sensitivity was assumed to extend to vanishing process. concentrations, as in the original Keller–Segel model , which is Previous analysis of bacterial traveling bands assumed that the not biologically realistic . 10,11,37,44,45 population consisted of identical cells or at most Following depletion of local resources, the spatial self- 24,25 two phenotypes . Here we extended these studies by taking organizing mechanism described here could enable populations into account the continuum of phenotypes that is always present of bacteria to maintain diversity while traveling toward better 13,21,22 in a population . One study that considered how two environments. This diversity increases the probability that a phenotypes might travel together made the theoretical assump- phenotype well-suited to unexpected environments will be avail- tion that cells sense only the direction of the gradient, not its able if needed during travel until a destination is reached where magnitude . This assumption causes the peak densities of the growth can replenish the population. As the range of phenotypes two phenotypes in that model to coincide in space, contrary to allowed within a traveling group depends on the spatial profile of our experimental observations. In another theoretical study , the the traveling gradient, this mechanism introduces important cells did respond to the gradient magnitude, and the two phe- feedback between the environment, cellular metabolism, and notypes in the traveling solution were spatially separated. phenotypic diversity, which together generate spatial patterns of Although not discussed in the paper, the phenotype with the phenotypes according to functional capabilities. The same higher chemotactic coefficient is in the front in that solution, in mechanism might also enable different bacterial species to travel agreement with our sorting mechanism. However, in that model, together, thus enabling migration of small ecosystems. NATURE COMMUNICATIONS (2018) 9:2177 DOI: 10.1038/s41467-018-04539-4 www.nature.com/naturecommunications 7 | | | Time (sec) –1 PDF Cell density (ml ) Distance between Time (sec) peak positions (mm) –1 Cell density (ml ) ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04539-4 order Runge–Kutta routine for temporal integration (time step 0.08 s). We used Collective migration of eukaryotes resulting from traveling no-flux boundary conditions. The initial condition was ρðÞ x; 0 ¼ gradients of attractants generated by consumption or breakdown ðÞ =x of an attractant has recently been found to be more important ρ  PðÞ TB e for 0 < x <1.6 mm and ρ (x,0) = 0for x >1.6 mm. Here, ρ is i 0 0 i than previously believed because it enables cell migration over the initial cell density scale, determined from the total initial cell number ~ 2 × 10 ; x = 0.8 mm; and P(TB ) was obtained from experimental measurements 0 i much larger distances than migration along externally imposed (Fig. 2eblack). gradients . Being able to maintain diversity within the traveling Assuming near-constant wave speed, we rewrite Eq. 2 in the moving group could be important in the context of immunity and coordinate z=x–ct and integrate from − ∞ to + ∞ to obtain cancer. cA hi α ¼ : ð4Þ The types of interactions between collective behavior and N=a phenotypic diversity reported here might also be at play beyond microbiology and cell biology, in contexts where individuals in a Here, α is the average attractant consumption rate, N is the number of cells in the band, and a is the cross-sectional area of the channel. In the absence of group respond to the cumulative effect on the environment of the oxygen-dependent consumption of Asp, the average consumption rate in the individuals ahead. Whereas a bird monitors its neighbors to band would be constant across experimental conditions . However, oxygen- benefit from the collective information acquired by the flock, here dependent consumption makes the average consumption rate decreases with individuals monitor the environmental gradient to benefit from increasing cell density in the band, which is correlated with the number of cells in the band (Fig. 2c). As shown in Supplementary Fig. 7h, the value of ⟨α ⟩ the information accumulated in the environment by the band. In A calculated from experimental data using Eq. 4 decreases as the number of cells in both cases, there is a group “memory” in the form of a spatial the band increases, which is captured in simulations. structure that individuals respond to . In general, both types of To derive the result of spatial sorting analytically from this model, we first memory are probably available and utilized by groups, with their rewrite Eq. 1 in the moving coordinate z = x − t: relative importance determined by the specific biology of the ∂ρ ∂ ρ ∂ ∂f ðAÞ i i c ¼ μ  χ ρ ; ð5Þ organisms. i i i ∂z ∂z ∂z ∂z Methods Noting that for each phenotype i to be traveling with the group, its density profile dρ Mathematical model. We extended the classic Keller–Segel model to include the i must have a peak. Around the density peak z we must have ¼ 0 and dz z¼z effect of phenotypic differences in TB. The key variables in the model are the d ρ density ρ (x, t) of cells of phenotype i as a function of position x and time t, and the i <0. dz z¼z concentration of Asp A(x, t). Cells of phenotype i are characterized by their che- i Rewriting Eq. 5 as motactic coefficient χ and diffusivity μ . As Asp consumption depends on oxy- i i gen , we also model the amount of oxygen dissolved in the media O(x, t). The 2 2 ∂ ρ ∂ρ ∂ρ ∂fAðÞ ∂ f ðAÞ i i i ð6Þ μ ¼c þ χ þ χ ; parameters of the model are in Table 1. The time-dependent evolution reads: i 2 i i 2 ∂z ∂z ∂z ∂z ∂z ∂ρ ∂ ρ ∂ ∂f ðAÞ i i ð1Þ ¼ μ  χ ρ we then have at the density peaks z : i i i i ∂t ∂x ∂x ∂x 2 2 d ρ d f μ ¼ χ ρ <0 ð7Þ i 2  i i 2 dz dz 2 z¼z z¼z X i i ∂A ∂ A A ¼ μ  α ðÞ O ρ ð2Þ A A i ∂t ∂x K þ A Eq. 7 indicates that in the front of the band the perceived gradient must be shallow and become progressively steeper as z decreases toward the back (Fig. 3c black). dρ 2 Integrating Eq. 5 and using  ¼ 0, we obtain dz ∂O ∂ O O z¼z ¼ μ  α ρ þ κðÞ O  O ð3Þ O O i ex ∂t ∂x K þ O df c ¼ χ ð8Þ dz z¼z Eqs. 7, 8 together show that the cell density peaks z of each phenotype i are Eq. 1 represents the motion of cells due to diffusion and chemotaxis. fAðÞ¼ hi monotonically ordered according to their chemotactic coefficients χ . 1þA=K Mlog is the perceived signal which depends on the local Asp concentration Examining the effect of oxygen analytically, at the back of the wave the Asp 1þA=K ∂fAðÞ M dA concentration is small so that the chemotactic drift there is χ  χ . A, where K = 3.5 μM and K = 1000 μM represent the dissociation constants of 0 1 i ∂z i K dz Asp for the inactive and active conformations of the Tar receptors, and M =6is Rewriting Eq. 2 in the moving coordinate, assuming the diffusion term is 29,48 the receptor gain . The effective diffusion coefficient and the chemotaxis negligible , and assuming A ≫ K gives an expression for dA/dz. From this, the coefficient are functions of the microscopic parameters of individual cell swimming drift becomes ðÞ 1TB 22,29 v i behavior. We draw on previous work to model them as μ ¼ i 3 ð1ΘÞλ þ2D R;i rot ∂fAðÞ χ M χ  α ðÞ O ρ : k TB ð9Þ 0 i −1 i A i and χ ¼ μ (Supplementary Fig. 8a). In these expressions, v ~36 μms is ∂z K c i TB þTB i 0 0 i i the cell swimming speed, which when projected in two dimesnions (2D) corresponds to the average speed we measure in our quasi-2D device (~ 28 μms At the back of the band there are fewer cells and therefore more oxygen than in the −1 22 ); Θ = 0.16 is the directional persistence between successive runs ; D = 0.062 middle of the band. Thus, when the Asp consumption rate depends on oxygen, rot qffiffiffiffiffiffiffiffiffiffi −1 21,50 TB α (O) becomes larger at the back than the mean over the band, ⟨α ⟩.Asa i A A s is the rotational diffusion coefficient during runs ; λ ¼ ω is the R;i 1TB consequence, the drift at the back is higher than in the case without oxygen −1 rate of switching from the run state to the tumble state; and ω = 3.8 s is the dependence, slowing down the decay of the band. We verified this by running effective switching frequency . μ and χ are monotonically decreasing functions of i i simulations with and without oxygen dependence (Supplementary Fig. 8b). In the the TB and χ /μ ≈ 22 = k , apart for the range 0 ≤ TB ≲ TB = 0.05, where χ(TB ) i i i 0 i 0 i simulations without oxygen, we set ⟨α ⟩ equal to a constant value, corresponding is increasing as recently observed in experiments where individual cells were to the average consumption rate in the band α in the simulation with oxygen tracked in a static gradient of α-methylaspartate (non-metabolizable analog of dependence. This was intended to make the wave speeds in the two simulations Asp) . similar, eliminating the effect of different wave speeds on cell leakage rates. Eq. 2 represents the change in the concentration of Asp due to diffusion, with diffusivity μ , and to consumption, with half-max constant K and maximum rate A A Strains, growth conditions, and sample preparation. E. coli RP437 was used as α ðÞ O ¼ α 1  g þ g . Here, α is the base consumption rate, O is the A0 ex A A0 A A O the wild type strain for chemotaxis in this study. Cells were grown in M9 glycerol ex −1 −1 −1 external oxygen concentration, and g is the fractional reduction of Asp medium: M9 salts (6.78 g L Na HPO , 3.0 g L KH PO , 0.5 g L NaCl, 1.0 g L A 2 4 2 4 −1 −1 consumption rate at zero oxygen. Eq. 3 describes the time-dependent evolution of NH Cl), supplemented with 4 mL L glycerol, 0.1 % casamino acids, 1.0 mM magnesium sulfate, and 0.05% w/v polyvinylpyrrolidone-40 at 30 °C. Appropriate oxygen, with diffusion coefficient μ , maximum consumption rate α , half-max O O −1 −1 constant K , and supply of oxygen through the polydimethylsiloxane (PDMS) with antibiotics were supplemented (ampicillin 100 µg mL , kanamycin 50 µg mL , −1 the mass transfer rate κ. and chloramphenicol 25 µg mL ) when necessary to maintain plasmids. Eqs. (1–3) were integrated in MATLAB using second-order centered For Fig. 1, cells were collected at mid-exponential phase (approximately an differences for the spatial derivatives (mesh size 20 µm) and an explicit fourth- OD600 of 0.3) and washed twice with fresh M9 glycerol medium, then resuspended 8 NATURE COMMUNICATIONS (2018) 9:2177 DOI: 10.1038/s41467-018-04539-4 www.nature.com/naturecommunications | | | NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04539-4 ARTICLE Table 1 Model parameters Symbol Definition and value Reference M The receptor gain for aspartate, M= 6 Ref. K The dissociation constant to aspartate for the inactive conformation of the Tar receptor, K = 3.5 μM Ref. 0 0 K The dissociation constant to aspartate for the active conformation of the Tar receptor, K = 1000 μM Ref. 1 1 2 −1 55 μ The diffusion coefficient of aspartate molecules, μ = 500 μm s Ref. A A K The aspartate concentration at half-max of its consumption, K = 0.5 μM Ref. A A −12 −1 −1 α The maximum aspartate consumption rate, α = 9.3 × 10 μmol cell Supplementary Fig. 7g A0 A0 O The external oxygen level, O = 250 μM Ref. ex ex g The basal ratio of relative consumption rate at zero oxygen, g = 0.27 Supplementary Fig. 7g A A 2 −1 55 μ The diffusion coefficient of dissolved oxygen, μ = 2500 μm s Ref. O O −11 −1 −1 α The maximum oxygen consumption rate, α = 7×10 μmol min cell Supplementary Fig. 7g O O K The dissolved oxygen concentration at half-max its consumption, K = 1 μM Supplementary Fig. 7g O O −1 58 κ The oxygen transfer rate through ~ 0.5 cm of PDMS, κ = 0.02 s Ref. ρ Initial cell density, ρ ~ 2.87 × 10 cells per mL This study 0 0 x Length scale of initial cell density profile, x = 0.8 mm This study 0 0 −1 v Cell swimming speed, v ~36 μms This study Θ Directional persistence Θ = 0.16 Ref. −1 50 D Rotational diffusion coefficient during runs, D = 0.062 s Ref. rot rot k Parameter in the expression for χ /μ , k = 22 Ref. 0 i i 0 TB Parameter in the expression for χ /μ,TB = 0.05 Ref. 0 i i 0 in fresh M9 glycerol medium to concentrate cell density at an OD600 of 0.7. These The laminated layers were then cut out and the remaining ports were punched to cells were then gently loaded into the microfluidics chamber, which was make external connections with the channels. To reduce the evaporation of the maintained at 30 °C throughout the experiment. microfluidic device, the PDMS device was soaked overnight in Millipore-filtered To generate a single traveling band, experiments were conducted in M9 glycerol water at 50 °C. buffer: motility buffer (M9 salts, 0.01 mM methionine, 0.1 mM EDTA, 0.05% w/v The assembled PDMS devices were bonded to 24 × 50 mm glass coverslips −1 polyvinylpyrrolidone-40) supplemented with 4 mL L glycerol and the indicated (#1.5). The PDMS was cleaned with transparent adhesive tape (Magic Tape, amount of Asp. This buffer was used to wash and resuspend cells instead of the Scotch) followed by rinsing with (in order) isopropanol, methanol, and Millipore- complex growth medium (M9 glycerol medium) mentioned above. The same filtered water, air-drying between each rinse. The glass was rinsed the same way, M9 salts were used in the M9 glycerol medium and M9 glycerol buffer to minimize with acetone, isopropanol, methanol, and Millipore-filtered water. The PDMS was osmolality changes. RP437 is auxotroph for leucine, histidine, methionine, and tape-cleaned an additional time, and then the two pieces were placed in a plasma threonine, and therefore does not grow in M9 glycerol buffer. bonding oven (Harrick Plasma) under vacuum, gently laminated, and then baked To control the TB distribution, we used a ΔcheZ strain derivative of RP437 on an 80 °C hotplate for 15 min to establish a covalent bond. Devices were stored at containing a chromosomally-integrated copy of the phosphatase CheZ under the room temperature and used within 24 h. control of the inducible promoter and tetR (gift from Dr. Chenli Liu). CheZ dephosphorylates the response regulator CheY, which, when phosphorylated, induces the motors to switch and the cells to tumble. Thus, low CheZ results in Band formation and imaging. Washed cells were gently loaded into the device higher CheY-P and more tumbling, and vice versa. aTc was added overnight in the which was then centrifuged for 20 min at 700 g in a 30 °C environmental room to culture when indicated to release the repressor TetR from the cheZ promoter concentrate cells at the end of the chamber (Supplementary Fig. 1b). After spin- region, inducing the expression of CheZ in this strain. To color-code strains, ning, the microfluidic device was placed on an inverted microscope (Nikon Eclipse pBca1020-r0040 carrying mRFP1 (obtained from BioBrick), pLambda driving Ti-U) equipped with a custom environmental chamber (50% humidity and 30 °C). mRFP1 (gift from Dr. Chenli Liu) and plasmids carrying YFP under constitutive A custom MATLAB script was used to control the microscope and its automated promoter were transformed into RP437 and into the inducible CheZ strain by stage (Prior) via the MicroManager interface . Time-lapse images (phase‐contrast electroporation. and fluorescence: mRFP1 or YFP) of the migrating cells were acquired using a Hamamatsu ORCA-Flash4.0 V2 camera (2,048 × 2,048 array of 6.5 × 6.5 μm pixels), a × 10 phase contrast objective (Nikon CFI Plan Fluor, N.A. 0.30, W.D. Microfluidic device design and fabrication. Microfluidic devices were con- 16.0 mm) and a LED illuminator (Lumencor SOLA light engine, Beaverton, OR) structed from the biocompatible and oxygen-permeable silicone polymer PDMS on through the mCherry block (Chroma 49008, Ex: ET560/× 40, Em: ET630/75 m) or cover glass following standard soft lithography protocols for two-layer devices . EYFP block (Chroma 49003; Ex: ET500/× 20, Em: ET535/30 m). Once the band The master molds for the device consisted of two silicon wafers with features formed, starting at the origin (closed end of the channel), the motorized stage created using ultraviolet (UV) photoresist lithography. The bottom wafer had two moved along the channel and paused every ~ 1.3 mm (the width of one frame with main parts: a large chamber created using SU-8-negative resist (thickness: 10 µm, a small overlap < 0.1 mm between consecutive positions) to take images in phase SU8 3010, Microchem) and a long channel together with two inlet/outlet channels contrast and fluorescence (exposure time 122 ms for both channels). After reaching designed to be opened and closed using pressure actuated valves. A second coat of the observation chamber, acquisition started over at the origin every 40 s (every 60 SPR-positive resist (thickness: 14 µm, SPR 220-7.0, MEGAPOSIT) on the same s for Fig. 1c). wafer was used to create a rounded channel profile that can collapse fully if In Fig. 1, all cells were expressing mRFP1. In Fig. 2, unlabeled and fluorescently depressed from above (Supplementary Fig. 1a). The second, top wafer contained labeled cells were separately grown to mid-exponential phase (OD600 of about 0.3). features for the control channels that close the collapsible features in the bottom For the experiments with 50, 100, and 200 µM Asp, the fluorescently labeled cells wafer. The top wafer was created using SU-8 negative resist (thickness: 10 µm, SU8 were diluted with unlabeled cells at the following ratios 1:20, 1:50, 1:100. The mixed 3010, Microchem). The resists were then cured using UV light exposure through cells were then washed and resuspended to the same predetermined density photomasks designed in CAD software and printed by CAD/Art Services, Inc. (OD600 of 0.7). The mixed populations were loaded into the microfluidic device (Bandon, Oregon), again following photoresist manufacturer specifications. and imaged as described above. In Fig. 4b-e, a similar procedure was followed to Subsequently, wafers were baked and the uncured photoresist was dissolved. After prepare the samples for different induction conditions. A 1:1 mixture of high and curing the SPR coat, the features were baked further to produce a rounded profile. low induction cultures was mixed and loaded in the device. After both wafers were complete, a protective coat of silane was applied by vapor In Fig.4f, a 1:1 mixture of high and low aTc induction (ranging from 1 to 10 ng −1 deposition. mL ) cultures with all cells fluorescently labeled by mRFP1 or YFP was mixed and To cast and manufacture the two-layer device, the top wafer was coated with a loaded in the device. The distances between the peaks of the density profiles of the 5 mm-thick layer of degassed 10:1 PDMS-to-curing agent ratio (Sylgard 184, Dow two populations were calculated by measuring the distance between the peaks in Corning). For the bottom layer, a 20:1 mixture was prepared and spin coated to the two fluorescent intensity profiles. The mean TB of each population were create a 100 μm-thick layer. The two layers were partially baked for 45 min at 70 °C. measured by loading a sample of the population on a cover slip and tracking The top layer was then cut and separated from the wafer, holes were punched from individual cells as previously described . the feature side using a sharpened 20-gauge blunt-tip needle to make external Once a band arrived in the perfused chamber, the gate near the chamber was connections to the control valve lines, then aligned and laminated onto the bottom closed (using 10 psi pressure) to capture the band (Supplementary Fig. 1b). To layer. The stacked layers were baked together for 1.5 h at 70 °C and allowed to cool. capture the second band in a separate experiment, the gate remained open until the NATURE COMMUNICATIONS (2018) 9:2177 DOI: 10.1038/s41467-018-04539-4 www.nature.com/naturecommunications 9 | | | ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04539-4 first band passed, then closed. The first band was immediately flushed away by fluorescence intensity, the opposite effect of decreased oxygen. Experiments were flowing M9 glycerol medium at 3 psi for 5 ~ 20 s through the chamber until all cells performed in 100 μm-deep and 14 μm-deep devices. The 100 μm-deep devices were were flushed from the chamber. The gate was then reopened to let the second band straight channels, 39 mm long and 600 μm wide. The 14 μm-deep devices were the migrate into the chamber. Once a band was captured in the perfused chamber, same as those used in the other experiments in this paper. However, as the control several pulses of fresh medium/buffer at 3 psi were flown in to reduce cell density gates interfere with the fluorescence signal, the top layer of PDMS was fabricated and homogenize the environment in the chamber. Cells were left to adapt in the without a mold for the gates. In both cases, a 1:500 dilution of concentrated perfused chamber for a few minutes. micelles was added to bacteria prepared as described above, just before loading into the microfluidic device. This dilution was intended to avoid sequestration of oxygen from the cells. Two control experiments were performed in the 14 μm TB detection. Once cell density became relatively homogeneous, the swimming device, one with dye and no cells, and one with cells and no dye. Imaging and data trajectories of individual cells were recorded for 2 min at 8 fps in order to extract analysis of all three types of experiments were performed in the same way. TB. We verified that perfusion of the observation chamber did not affect the Imaging equipment were the same as described above with a few exceptions. distribution of TBs (Supplementary Fig. 3). Each TB distribution was generated by The excitation filter from an ECFP block (Chroma 31044v2, D436/× 20), the acquiring four movies, tracking the individual cells, and determining their TB as emission filter from an mCherry block (Chroma 49008, ET630/75 m), and the 21,22 described before . The resulting number of sample trajectories longer than 10 s dichroic mirror from the mCherry block (Chroma 49008, T585lpxr) were used to (shorter tracks were discarded, because they provide poor estimate of TB ) was image the dye due to its large Stokes shift. limited by the number of independent movies one can acquire in the PDMS For the 100-μm device, the × 10 objective described above was used for imaging. chamber, the size of which corresponds to four field of views, and the density of Images were taken in 2 min intervals. For the 14-μm device, imaging was cells, which must be kept low for tracking to be possible. This procedure resulted in performed using a × 40 oil objective (Nikon CFI Plan Fluor, NA 1.3, W.D. 0.24 a minimum of 2823 and 3209 trajectories per distribution giving an error of at mm). Images were taken in 3.5 min intervals. most 1.15% and 1.00% in the determination of the cumulative distribution function To analyze the fluorescence images, fluorescence intensity was first averaged (CDF) estimated by bootstrapping for Figs. 1 and 2, respectively. over the width of the channel that was visible in the image, I(x, t). At each location x, the passing wave appeared as a brief peak in intensity during the time course. To Determination of the number of cells in the band. Image analysis was conducted separate this from slow variations in signal due to photobleaching and possible in MATLAB. We detected the position of the centroid of each fluorescent cell using global changes in oxygen concentration, we smoothed the time course of the MATLAB function bwconncomp. Figures 2a, 4b, and Supplementary Figs. 6d, e fluorescence intensity at each position, I (t), using MATLAB’s smooth function report these coordinates. The number of labeled cells was multiplied with the (smoothing method: lowess; window size: 4, corresponding to a time window of 14 dilution ratio to obtain the total number of cells in the band. min) to produce the slowly varying background signal, I ðtÞ. To extract the fast- Cell density profiles in Fig. 2d, 4c d, and Supplementary Figs. 9b, c were passing wave, we divided I (t) element wise by I ðtÞ for each position x. As a result, x x measured as follows: cell density profiles were extracted for each time point within the slowly changing background was normalized to 1, while faster changes in signal were different from 1. This also eliminated differences in illumination over space. one experiment and aligned before averaging. The cell density profile at a given time point was calculated by dividing the number of cells by the volume in one Finally, the noisy normalized profiles produced by this analysis were median filtered over space using MATLAB’s medfilt1 (window size: 10). spatial bin (~ 120 µm) along the observation channel. To reduce alignment error before averaging, the position of the peak density in each profile was identified by first smoothing the profile with a moving filter with a 5-bin span (MATLAB Oxygen consumption rate in batch cultures. RP437 cells were grown in 200 mL function smooth) and then identifying the position of the peak. To avoid boundary M9 glycerol medium to mid-exponential phase and washed twice with M9 glycerol effects, only the profiles with peak position located between 3 and 8 mm from the buffer supplemented with 200 µM Asp. Cells were then resuspended in 50 mL of origin were used to calculate the average density profile. The mean and SD shown the same defined buffer to an OD600 of 0.5. The sample was placed in a beaker at in the figures were calculated with raw cell density profiles (not smoothed). 30 °C. The surface of the buffer was sealed by overlaying mineral oil. The level of dissolved oxygen was measured with a portable dissolved oxygen meter (Mel- Measurement of amino acids by HPLC. When consuming Asp, E. coli cells waukee MW600) every 30 s. The cell sample was continuously stirred at 300 r.p.m. secrete other amino acids, which could affect group migration . We used HPLC to with a magnetic stirrer. The consumption rate per cell per minute was obtained by analyze the amino acids secreted by the cells when they are suspended in M9 dividing the reduction in oxygen by the number of cells and the sampling interval glycerol buffer supplemented with Asp. RP437 cells were grown in 200 mL M9 time (Supplementary Fig. 7g). glycerol medium up to mid-exponential phase and washed twice with M9 glycerol buffer supplemented with 500 µM Asp. Cells were then resuspended in 5 mL of the same defined buffer at an OD600 of 1 and placed in a 200 mL flask. The flask was Asp consumption rate in batch cultures. RP437 cells were grown in 200 mL M9 shaken at 200 r.p.m. to maximize aeration at 30 °C. Every 15 min, a 500 µl of glycerol medium up to mid-exponential phase and washed twice with M9 glycerol culture was sampled, filtered using 0.2 µm filter (Acrodisc 13 mm Syringe Filter buffer supplemented with 500 µM Asp. Cells were resuspended in 3 mL of the same with 0.2 µm HT Tuffryn Membrane, Pall Corporation), and analyzed by HPLC via defined buffer to make cell density at an OD600 of 1. Two samples were prepared pre-column derivatization method. The resulting derivatives were separated by in two test tubes. One tube was shaken at 200 r.p.m. to maximize aeration, whereas phase chromatography using a Dionex Ultimate 3000 HPLC, with a coupled DAD- the other tube was left on the bench with mineral oil overlaid on the liquid surface 3000RS diode array detector (Dionex) and FLD detector (Dionex) using an ACE to avoid supply of oxygen from air. These two tubes were incubated at 30 °C and C18 column (3 µm, 3 × 150 mm). Amino acid standard (AAS18 Sigma) was used as sampled every 10 min. The amount of Asp in the collected sample was then reference. measured by using Aspartate Assay Kit (Abcam). The consumption rate per cell Upon uptake of Asp, cells secreted small amounts of glutamate (Glu), per minute was obtained by dividing the reduction of Asp in the sample by the asparagine (Asn), and homoserine (HS), which are attractants (Supplementary number of cells and the sampling interval time (Supplementary Fig. 7g). Fig. 6c). To quantify the relative contribution of each amino acid to the chemotactic response in our experiment the measured concentration of each amino Statistical analysis and experimental reproducibility. No statistical methods acid was plotted in units of the corresponding EC50 of the dose response of the were used to predetermine sample size. Standard error in the CDF of the TB in chemotaxis system. The EC50 of the dose response of the chemotaxis system for each replicate experiment was determined by bootstrapping (1000 bootstrap each amino acid has been quantified by the Sourjik lab using in vivo Förster samples). One-sided unpaired two-sample Student’s t-test assuming unequal var- resonance energy transfer measurements in RP437, the strain used in this study. iances was used for comparison between two groups in Fig. 2b, c. P-values < 0.05 The EC50 values of RP437 E. coli are 0.3 μM for Asp, 50 μM for Glu, 30 μM for were considered statistically significant and marked with asterisks. For t-test, t- Asn, and 3 mM for HS . Supplementary Fig. 6c reveals that the response to Asp dominates by almost two orders of magnitude over the responses to the other values and degrees of freedom are provided in the figure legends. The error bars are defined in each figure caption and are standard deviation except in Fig. 2f. Data amino acids. It also shows that of the three secreted amino acids, Glu is the one that has the second largest effect, albeit still much smaller than the response to Asp. presented in the main figures were drawn from at least three independent repli- cates, with the exception of Fig. 4d(n = 2). The number of replicates is mentioned Finally, we also checked that chemotaxis towards oxygen does not play a significant role either. Mutant strains lacking the oxygen receptor aer or both aer and tsr form in the caption of each figure. bands under the same condition as in Fig. 2a, indicating that aerotaxis is not essential and Tar response to Asp is sufficient for the band to travel Code availability. To analyze the swimming behavior of E. coli cells we used (Supplementary Fig. 6d, e). 21,22 custom MATLAB code as reported in refs. . The code to simulate the math- ematical model is described above and available from the corresponding author Measurement of oxygen in the center of the traveling band. Ruthenium upon request. complexes are toxic to E. coli; hence, the need to encapsulate them in phospholipid micelles. This was achieved using the same protocol as in ref. . Fluorescence of the ruthenium complex is quenched by oxygen binding, so higher fluorescence Data availability. Data for each figure is provided as a MATLAB .fig file from corresponds to lower oxygen concentration. It is noteworthy that the high density which the data points can be extracted https://doi.org/10.6084/m9. of cells in the band could exclude the dye; however, this should decrease figshare.6207371. 10 NATURE COMMUNICATIONS (2018) 9:2177 DOI: 10.1038/s41467-018-04539-4 www.nature.com/naturecommunications | | | NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04539-4 ARTICLE Received: 25 July 2017 Accepted: 3 May 2018 35. Drescher, K., Dunkel, J., Cisneros, L. H., Ganguly, S. & Goldstein, R. E. Fluid dynamics and noise in bacterial cell-cell and cell-surface scattering. Proc. Natl Acad. Sci. USA 108, 10940–10945 (2011). 36. Zhang, H. P., Be’er, A., Florin, E. L. & Swinney, H. L. Collective motion and density fluctuations in bacterial colonies. Proc. Natl Acad. Sci. USA 107, 13626–13630 (2010). References 37. Brenner, M. P., Levitov, L. S. & Budrene, E. O. Physical mechanisms for 1. Sumpter D. J. T. Collective Animal Behavior. Princeton Univ. Press (2010). chemotactic pattern formation by bacteria. Biophys. J. 74, 1677–1693 (1998). 2. Liu, J. et al. Metabolic co-dependence gives rise to collective oscillations within 38. Holz, M. & Chen, S. H. Quasi-elastic light scattering from migrating biofilms. Nature 523, 550–554 (2015). chemotactic bands of Escherichia coli. Biophys. J. 23,15–31 (1978). 3. Couzin, I. D. & Krause, J. Self-organization and collective behavior in 39. Scribner, T. L., Segel, L. A. & Rogers, E. H. A numerical study of the formation vertebrates. Adv. Stud. Behav. 32,1–75 (2003). and propagation of traveling bands of chemotactic bacteria. J. Theor. Biol. 46, 4. Vicsek, T. & Zafeiris, A. Collective motion. Phys. Rep. 517,71–140 (2012). 189–219 (1974). 5. Potts, W. K. The chorus-line hypothesis of maneuver coordination in avian 40. Novick-Cohen, A. & Segel, L. A. A gradually slowing travelling band of flocks. Nature 309, 344–345 (1984). chemotactic bacteria. J. Math. Biol. 19, 125–132 (1984). 6. Ballerini, M. et al. Interaction ruling animal collective behavior depends on 41. Wong-Ng, J., Celani, A. & Vergassola, M. Exploring the function of bacterial topological rather than metric distance: evidence from a field study. Proc. Natl chemotaxis. Curr. Opin. Microbiol. 45,16–21 (2018). Acad. Sci. USA 105, 1232–1237 (2008). 42. Lapidus, I. R. & Schiller, R. A model for traveling bands of chemotactic 7. Beyerinck, M. W. Ueber atmungsfiguren beweglicher bakterien. Zentr bacteria. Biophys. J. 22,1–13 (1978). Bakteriol. Parasite. 114, 827 (1893). 43. Lauffenburger, D., Kennedy, C. R. & Aris, R. Traveling bands of chemotactic 8. Adler, J. Chemotaxis in bacteria. Science 153, 708–716 (1966). bacteria in the context of population-growth. B Math. Biol. 46,19–40 (1984). 9. Adler, J. Effect of amino acids and oxygen on chemotaxis in Escherichia coli. J. 44. Wang, Z. A. Mathematics of traveling waves in chemotaxis review paper. Bacteriol. 92, 121–129 (1966). Discret. Cont. Dyn.-B 18, 601–641 (2013). 10. Keller, E. F. & Segel, L. A. Traveling bands of chemotactic bacteria: a 45. Erban, R. & Othmer, H. G. From individual to collective behavior in bacterial theoretical analysis. J. Theor. Biol. 30, 235–248 (1971). chemotaxis. SIAM J. Appl. Math. 65, 361–391 (2005). 11. Saragosti, J. et al. Directional persistence of chemotactic bacteria in a traveling 46. Tweedy, L., Knecht, D. A., Mackay, G. M. & Insall, R. H. Self-generated concentration wave. Proc. Natl Acad. Sci. USA 108, 16235–16240 (2011). chemoattractant gradients: attractant depletion extends the range and 12. Wolfe, A. J. & Berg, H. C. Migration of bacteria in semisolid agar. Proc. Natl robustness of chemotaxis. PLoS Biol. 14, e1002404 (2016). Acad. Sci. USA 86, 6973–6977 (1989). 47. Hein, A. M. et al. The evolution of distributed sensing and collective 13. Spudich, J. L. & Koshland, D. E. Jr. Non-genetic individuality: chance in the computation in animal populations. eLife 4, e10955 (2015). single cell. Nature 262, 467–471 (1976). 48. Yang, Y. et al. Relation between chemotaxis and consumption of amino acids 14. Mayor, R. & Etienne-Manneville, S. The front and rear of collective cell in bacteria. Mol. Microbiol. 96, 1272–1282 (2015). migration. Nat. Rev. Mol. Cell Biol. 17,97–109 (2016). 49. Ford,R.M.&Cummings, P. T. On therelationshipbetween cell balanceequations 15. Kussell, E. & Leibler, S. Phenotypic diversity, population growth, and for chemotactic cell populations. SIAM J. Appl. Math. 52, 1426–1441 (1992). information in fluctuating environments. Science 309, 2075–2078 (2005). 50. Berg, H. C. & Brown, D. A. Chemotaxis in Escherichia coli analysed by three- 16. Veening, J. W., Smits, W. K. & Kuipers, O. P. Bistability, epigenetics, and bet- dimensional tracking. Nature 239, 500–504 (1972). hedging in bacteria. Annu. Rev. Microbiol. 62, 193–210 (2008). 51. Demir, M. & Salman, H. Bacterial thermotaxis by speed modulation. Biophys. 17. Frankel, N. W. et al. Adaptability of non-genetic diversity in bacterial J. 103, 1683–1690 (2012). chemotaxis. eLife 3, e03526 (2014). 52. Xia, Y. W. G. Soft lithography. Annu. Rev. Mater. Sci. 28, 153–184 (1998). 18. Ackermann, M. A functional perspective on phenotypic heterogeneity in 53. Edelstein, A. D. et al. Advanced methods of microscope control using microorganisms. Nat. Rev. Microbiol. 13, 497–508 (2015). muManager software. J. Biol. Methods 1, e10 (2014). 19. Copenhagen, K., Quint, D. A. & Gopinathan, A. Self-organized sorting limits 54. Adler, M., Erickstad, M., Gutierrez, E. & Groisman, A. Studies of bacterial behavioral variability in swarms. Sci. Rep. 6, 31808 (2016). aerotaxis in a microfluidic device. Lab. Chip. 12, 4835–4847 (2012). 20. Waite, A. J., Frankel, N. W. & Emonet, T. Behavioral variability and 55. Hazel, J. R. & Sidell, B. D. A method for the determination of diffusion phenotypic diversity in bacterial chemotaxis. Annu. Rev. Biophys. 47, coefficients for small molecules in aqueous solution. Anal. Biochem. 166, 27.21–27.22 (2018). 335–341 (1987). 21. Dufour, Y. S., Gillet, S., Frankel, N. W., Weibel, D. B. & Emonet, T. Direct 56. Schellenberg, G. D. & Furlong, C. E. Resolution of the multiplicity of the correlation between motile behavior and protein abundance in single cells. glutamate and aspartate transport systems of Escherichia coli. J. Biol. Chem. PLoS Comput. Biol. 12, e1005041 (2016). 252, 9055–9064 (1977). 22. Waite, A. J. et al. Non-genetic diversity modulates population performance. 57. Wetzel R. G. in Limnology: Lake and River Ecosystems, 3rd edn (Academic Mol. Syst. Biol. 12, 895 (2016). Press, 2001). 23. Ben-Jacob, E., Cohen, I. & Levine, H. Cooperative self-organization of 58. Vollmer, A. P., Probstein, R. F., Gilbert, R. & Thorsen, T. Development of an microorganisms. Adv. Phys. 49, 395–554 (2000). integrated microfluidic platform for dynamic oxygen sensing and delivery in a 24. Tai-Chia, L. & Zhi-An, W. Development of traveling waves in an interacting flowing medium. Lab. Chip. 5, 1059–1066 (2005). two-species chemotaxis model. Discret. Contin. Dyn. Syst. A 34, 2907–2927 (2014). 25. Emako, C., Gayrard, C., Buguin, A., Neves de Almeida, L. & Vauchelet, N. Acknowledgements Traveling pulses for a two-species chemotaxis model. PLoS Comput. Biol. 12, We thank H. Salman for sharing the protocol to synthesize the Ru-micelles; Y. Dufour e1004843 (2016). for help with the microfluidics; F. Isaacs for sharing BioBrick plasmids; C. Liu, S. Par- 26. Mittal, N., Budrene, E. O., Brenner, M. P. & Van Oudenaarden, A. Motility of kinson, and H. Salman for sharing E. coli strains and plasmids; N. Clay, B. Barco, and W. Escherichia coli cells in clusters formed by chemotactic aggregation. Proc. Natl Chezem for help with HPLC experiments; D. Clark and J. Wang for help with statistics;, Acad. Sci. USA 100, 13259–13263 (2003). and A. Waite, J. Howard, C. Brannon, A. Greene, and P. Turner for comments on the 27. Celani, A. & Vergassola, M. Bacterial strategies for chemotaxis response. Proc. manuscript. This study was supported by the National Institutes of Health grant Natl Acad. Sci. USA 107, 1391–1396 (2010). 1R01GM106189 and the Allen Distinguished Investigator Program (grant 11562) 28. Si, G., Wu, T., Ouyang, Q. & Tu, Y. Pathway-based mean-field model for through The Paul G. Allen Frontiers Group. C.H. acknowledges support through the Escherichia coli chemotaxis. Phys. Rev. Lett. 109, 048101 (2012). National Natural Science Foundation of China (31570095, 11504399). 29. Dufour, Y. S., Fu, X., Hernandez-Nunez, L. & Emonet, T. Limits of feedback control in bacterial chemotaxis. PLoS Comput. Biol. 10, e1003694 (2014). Author contributions 30. Sourjik, V. & Berg, H. C. Functional interactions between receptors in X.F. and S.K. contributed equally to this work. H.M. and J.L. contributed equally to this bacterial chemotaxis. Nature 428, 437–441 (2004). work. T.E., X.F., and S.K. designed the research. X.F., S.K., H.M., and C.H. performed the 31. Long, Z., Quaife, B., Salman, H. & Oltvai, Z. N. Cell-cell communication experiments. X.F., S.K., J.L., H.M., and T.E. performed the data analysis. X.F., J.L., H.M., enhances bacterial chemotaxis toward external attractants. Sci. Rep. 7, 12855 and T.E. performed the mathematical modeling. D.C.V. and S.W.Z. contributed to (2017). mathematical modeling. T.E., X.F., S.K., H.M., and J.L. wrote the manuscript. All authors 32. Mesibov, R. & Adler, J. Chemotaxis toward amino acids in Escherichia coli. discussed the manuscript. J. Bacteriol. 112, 315–326 (1972). 33. Douarche, C., Buguin, A., Salman, H. & Libchaber, A. E. coli and oxygen: a motility transition. Phys. Rev. Lett. 102, 198101 (2009). Additional information 34. Wu, X. L. & Libchaber, A. Particle diffusion in a quasi-two-dimensional Supplementary Information accompanies this paper at https://doi.org/10.1038/s41467- bacterial bath. Phys. Rev. Lett. 84, 3017–3020 (2000). 018-04539-4. NATURE COMMUNICATIONS (2018) 9:2177 DOI: 10.1038/s41467-018-04539-4 www.nature.com/naturecommunications 11 | | | ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04539-4 Competing interests: The authors declare no competing interests. appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party Reprints and permission information is available online at http://npg.nature.com/ material in this article are included in the article’s Creative Commons license, unless reprintsandpermissions/ indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in regulation or exceeds the permitted use, you will need to obtain permission directly from published maps and institutional affiliations. the copyright holder. To view a copy of this license, visit http://creativecommons.org/ licenses/by/4.0/. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give © The Author(s) 2018 12 NATURE COMMUNICATIONS (2018) 9:2177 DOI: 10.1038/s41467-018-04539-4 www.nature.com/naturecommunications | | |

Journal

Nature CommunicationsSpringer Journals

Published: Jun 5, 2018

References

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off