Hypoxia triggers collective aerotactic migration in Dictyostelium discoideum
Abstract
Using a selfgenerated hypoxic assay, we show that the amoeba Dictyostelium discoideum displays a remarkable collective aerotactic behavior. When a cell colony is covered, cells quickly consume the available oxygen (O_{2}) and form a dense ring moving outwards at constant speed and density. To decipher this collective process, we combined two technological developments: porphyrinbased O_{2} sensing films and microfluidic O_{2} gradient generators. We showed that Dictyostelium cells exhibit aerotactic and aerokinetic response in a low range of O_{2} concentration indicative of a very efficient detection mechanism. Cell behaviors under selfgenerated or imposed O_{2} gradients were modeled using an in silico cellular Potts model built on experimental observations. This computational model was complemented with a parsimonious ‘Go or Grow’ partial differential equation (PDE) model. In both models, we found that the collective migration of a dense ring can be explained by the interplay between cell division and the modulation of aerotaxis.
Introduction
Oxygen is the main electron acceptor for aerobic organism to allow efficient ATP synthesis. This highenergy metabolic pathway has contributed to the emergence and diversification of multicellular organism (Chen et al., 2015). While high O_{2} availability in the environment seems a given, its rapid local consumption can generate spatial and temporal gradients in many places, including within multicellular organism. Oxygen level and gradients are increasingly recognized as a central parameter in various physiopathological processes (Tonon et al., 2019), cancer and development. The wellknown HIF (hypoxiainducible factor) pathway allows cells to regulate their behavior when exposed to hypoxia. At low O_{2} levels, cells accumulate HIFα leading to the expression of genes that support cell functions appropriate to hypoxia (Pugh and Ratcliffe, 2017).
Another strategy used by organisms facing severe oxygen conditions is to move away from hypoxic regions, a mechanism called aerotaxis and first described in bacteria (Engelmann, 1881; Winn et al., 2013). Aerotaxis will occur at the interface between environments with different oxygen content, such as soil/air, water/air or even within eukaryotic multicellular organisms between different tissues (Lyons et al., 2014). In such organisms, oxygen was proposed to be a morphogen as in placentation (Genbacev et al., 1997) or a chemoattractant during sarcoma cell invasion (Lewis et al., 2016). Aerotaxis may also play a role in morphogenesis. The notion that gradients of O_{2} and energy metabolism govern spatial patterning in various embryos dates back to the classic work of Child, 1941. Such notions have mostly been abandoned due to the inability to visualize such a gradient or clarify whether they are the result or the cause of developmental patterning (Coffman and Denegre, 2007). Even at the singlecell level, in vitro experimental studies on aerotaxis are rare. One reason might be technical: gradient control and live monitoring of oxygen concentrations at the cellular level are difficult. More recently, Chang et al. found an asymmetric distribution of hypoxiainducible factor regulating dorsoventral axis establishment in the early sea urchin embryo (Chang et al., 2017). Interestingly, they also found evidence for an intrinsic hypoxia gradient in embryos, which may be a forerunner to dorsoventral patterning.
Selfgenerated chemoattractant gradients have been reported to trigger the dispersion of melanoma cells (MuinonenMartin et al., 2014; Stuelten, 2017), Dictyostelium cells (Tweedy et al., 2016) or the migration of the zebrafish lateral line primordium (Donà et al., 2013; Venkiteswaran et al., 2013). The mechanism is simple and very robust: the cell colony acts as a sink for the chemoattractant, removes it by degradation or uptake creating a gradient that, in turn, attracts the cells as long as the chemoattractant is present in the environment. Physiologically speaking, selfgenerated gradients have been demonstrated to increase the range of expansion of cell colonies (Cremer et al., 2019; Tweedy and Insall, 2020) and to serve as directional cues to help various cell types navigate complex environments, including mazes (Tweedy et al., 2020). Recently, it was demonstrated that after covering an epithelial cell colony by a coverglass non permeable to O_{2}, peripheral cells exhibit a strong outward directional migration to escape hypoxia from the center of the colony (Deygas et al., 2018). This is a striking example of a collective response to a selfgenerated oxygen gradient by eukaryotic cells. Oxygen selfgenerated gradients could therefore play important roles in a variety of contexts, such as development, cancer progression, or even environmental navigation in the soil.
Dictyostelium discoideum (Dd) is an excellent model system to study the fairly virgin field of aerotaxis and of selfgenerated gradients. Dd is an obligatory aerobic organism that requires at least 5% O_{2} to grow at optimal exponential rate (Cotter and Raper, 1968; Sandonà et al., 1995) while slower growth can occur at 2% O_{2}. However, its ecological niche in the soil and around large amount of bacteria can result in reduced O_{2} availability. During its multicellular motile stage, high oxygen level is one of the signal used to trigger culmination of the migrating slug (Xu et al., 2012). For many years, Dd has been a classical organism to study chemotaxis and has emulated the development of many models of emergent and collective behavior since the seminal work of Keller and Segel (Hillen and Painter, 2009; Keller and Segel, 1970). An integrated approach combining biological methods (mutants), technological progress, and mathematical modeling is very valuable to tackle the issue of aerotaxis.
In this article, we study the influence of O_{2} selfgenerated gradients on Dd cells. Using a simple confinement assay, microfluidic tools, original oxygen sensors and theoretical approaches, we show that oxygen selfgenerated gradients can direct a seemingly collective migration of a cell colony. Our results confirm the remarkable robustness and longlasting effect of selfgenerated gradients in collective migration. This case where oxygen is the key driver also suggests that selfgenerated gradients are widespread and a possible important feature in multicellular morphogenesis.
Results
Confinement triggers formation and propagation of a selfsustained cell ring
To trigger hypoxia on a colony of Dd cells, we used a vertical confinement strategy (Deygas et al., 2018). A spot of cells with a radius of about 1 mm was deposited and covered by a larger glass coverslip with a radius of 9 mm. We measured the vertical confinement through confocal microscopy and found the height between the bottom of the plate and the coverslip to be 50 μm (Figure 1—figure supplement 1).
Using spots containing around 2000 cells (initial density around 10^{3} cells/mm^{2}), the formation of a dense ring of cells moving outwards was observed as quickly as 30 min after initiation of the confinement (Figure 1A and Video 1). This formation time however depended nonlinearly on initial cell density (the denser, the faster, Figure 1—figure supplement 2). Once triggered, this collective migration was selfmaintained for tens of hours, even days and the ring could, at these points, span centimeters (Figure 1B).
Notably, as the ring expanded outwards, it left a trail of cells behind. This led to the formation of a central zone populated by cells which did not contribute directly to the migration of the ring (Figure 1B) but were still alive and moving albeit a clear elongated phenotype resembling preaggregative Dd cells (Figure 1C and Video 2). In comparison, cells in the ring or outside the colony were rounder, as usual vegetative cells (DelanoëAyari et al., 2008).
To study the properties of the ring, we computed density profiles using radial coordinates from the center of the colony to study cell density as a function of time and distance to the center (Figure 1D–E). We found that after a transitory period corresponding to the ring passing through the initial spot, the density in the ring, its width and its speed all remained constant over long time scales (Figure 1—figure supplement 3). The speed and density of the ring were found to be 1.2 ± 0.3 μm/min (mean ± std, N=9 independent experiments) and 1.9 10^{3} ± 0.3 10^{3} cells/mm^{2} (mean ± std, N=4 independent experiments, that is fourfold higher than behind it, Figure 1E) respectively. The density of cells left behind the ring was also found to remain constant after a transient regime (Figure 1D). As the diameter of the ring increased over time, the absence of changes in morphology implies an increase of the number of cells and thus an important role of cell division.
Overall, this selfsustained ring propagation is very robust and a long lasting collective phenotype that can easily be triggered experimentally. This shows that the spot assay is an excellent experimental system to study the response of a variety of cell types to vertical confinement and its physiological consequences (Deygas et al., 2018).
Cell dynamics in different regions
Following the reported shape differences, we questioned how cells behaved dynamically in different regions. To do so, we performed higher resolution, higher frame rate experiments to allow cell tracking over times on the order of tens of minutes. Both the cell diffusion constant and instantaneous cell speeds were fairly constant throughout the entire colony (Figure 1—figure supplement 4). Cell diffusion was 28.2 ± 1.4 μm^{2}/min (N=3 independent experiments, each containing at least 2000 cells), comparable to our measurement of activity at very low oxygen level in the microfluidic device (see below). To test the influence of motion bias, we projected cell displacements on the radial direction and computed mean speeds in this direction as a function of distance to the center. Random motion, either persistent or not, would lead to a null mean radial displacement whereas biased migration would be revealed by positive (outward motion) or negative (inward motion) values. Here, we found that significantly nonzero biases were observed only in a region spanning the entire ring and a few tens of microns behind and in front of it with the strongest positive biases found in the ring (Figure 1—figure supplement 5).
Overall, our results show that the different regions defined by the ring and its dynamics can be characterized in terms of cell behavior: (i) behind the ring in the hypoxic region: elongated shape, normal speeds, and low bias; (ii) in the ring: round shape, normal speeds and high bias.
Response of Dd cells to controlled oxygen gradients
The spot assay is experimentally simple but is not ideally suited to decipher the response of Dd cells to hypoxia since local concentrations and gradients of oxygen are coupled to cell dynamics and thus very difficult to manipulate. We hence designed a new doublelayer PDMS microfluidic device allowing to quickly generate oxygen gradients in a continuous, controlled manner (Figure 2A). Briefly, cells were seeded homogenously within a media channel positioned 500 µm below two gas channels continuously perfused with pure nitrogen on one side and air on the other. As PDMS is permeable to these gases, the gas flows imposed hypoxic conditions on one side of the media channel while the other was kept at atmospheric oxygen concentration. Of note, the distance between the two gas channels, thereafter called the gap, varied from 0.5 mm to 2 mm in order to modify the steepness of the gradients in the median region of the media channels (Figure 2A and Materials and methods).
To make sure that the gas flows were sufficient to maintain a constant O_{2} distribution against leakages and against small variation in the fabrication process, we also developed O_{2}sensing films to be able to experimentally measure O_{2} profiles both in the microfluidic devices and in the spot assay. These films consisted of porphyrin based O_{2} sensors embedded in a layer of PDMS. As O_{2} gets depleted, the luminescence quenching of the porphyrin complex is reduced leading to an increase in fluorescence intensity (Ungerböck et al., 2013). Quantitative oxygen measurements were then extracted from this fluorescent signal using a SternVolmer equation (see Materials and methods and Figure 2—figure supplements 1–4 for details).
Within 15 min, we observed the formation of a stable O_{2} gradient in the devices closely resembling numerical predictions with or without cells (Figure 2B and Figure 2—figure supplements 5–7).
We then turned our attention to the reaction of the cells to this external gradient. We first noticed that depending on local O_{2} concentrations, cell motility was remarkably different. Using cell tracking, we found that cell trajectories seemed much longer and more biased in hypoxic regions (Figure 2C). These aerokinetic (large increase in cell activity) and aerotactic responses were confirmed by quantifying the mean absolute distance travelled by cells (Figure 2D top), or the mean distance projected along the gradient direction (Figure 2D bottom) in a given time as a function of position in the device (Figure 2D). Since cells in the microfluidic devices were also experiencing oxygen gradients, we further tested if the observed was true aerokinesis. To do so, we compared cell motility in homogenous environments of either 20.95% or 0.4% O_{2}. We found cell diffusion constant to be D=40.2±9.6 µm^{2}/min (mean ± std) at 0.4% (Figure 2—figure supplement 8), comparable to our measurements in the center of the spot (Figure 1—figure supplement 4). At atmospheric oxygen concentrations though, this effective diffusion was clearly reduced as we measured it to be D=19.2±8.8 µm^{2}/min (Figure 2—figure supplement 8). The very significant difference (p<0.0001) demonstrates that Dd cells show an aerokinetic positive response to low oxygen, even in the absence of gradients. The second important observation stemming from the microfluidic experiments is an accumulation of cells at some midpoint within the cell channel (Figure 2E). Naively, one could have expected cells to follow the O_{2} gradient over its entire span leading to an accumulation of all cells on the O_{2} rich side of the channel. This did not happen and, instead, cells seemed to stop responding to the gradient at a certain point. Similarly, we observed a strong positive bias in hypoxic regions but the bias quickly fell to 0 as cells moved to oxygen levels higher than about 2% (Figure 2D), confirming that the observed cell accumulation was a result of differential migration and not, for example, differential cell division. In addition, if we inverted the gas channels halfway through the experiment, we observed that the cells responded in around 15 min (which is also the time needed to reestablish the gradient, see Figure 2—figure supplement 6) and showed the same behavior, albeit in reverse positions. We measured the bias for the different gaps and for the situation of reversed gradient and obtained a value of 1.1 ± 0.4 µm/min (N=6, three independent experiments and for each, both directions of the gradient, each value stemming from a few hundred cells).
Of note, the position at which cells accumulated and stopped responding to the gradient was still in the region were the gradient was constantly increasing. This led to the hypothesis that, in addition to gradient strength, O_{2} levels also play an important role in setting the strength of aerotaxis displayed by Dd cells.
Furthermore, when we compared experiments performed with different gaps, we found that the position of cell accumulation varied from one channel to another (Figure 2E). However, our O_{2} sensors indicated that the accumulation occurred at a similar O_{2} concentration of about 1% in all three channels (inset of Figure 2B) thus strongly hinting that the parameter controlling the aerotactic response was O_{2} levels.
Overall, these experiments in controlled environments demonstrated two main features of the response of Dd cells to hypoxia: a strong aerokinetic response and a positive aerotactic response, both modulated by local O_{2} levels regardless of the local gradient. These results reveal a subtle cross talk between O_{2} concentrations and gradients in defining cell properties and it would be very informative, in the future, to study in details the reaction of Dd cells to various, well defined hypoxic environments where O_{2} concentrations and gradients can be independently varied.
Coupled dynamics between oxygen profiles and collective motion
Thanks to these results, we turned our attention back to the collective migration of a ring of cells and asked whether similar aerotactic behaviors were observed under selfgenerated gradients. To do so, we performed spot experiments on the O_{2}sensing films described above which allowed us to image, in parallel, cell behavior and O_{2} distribution (Figure 3A, Figure 2—figure supplement 3 and Video 3).
In a first phase, preceding the formation of the ring, cell motion was limited and the structure of the colony remained mostly unchanged. As O_{2} was consumed by cells, depletion started in the center and sharp gradients appeared at the edges of the colony (Figure 3B–C).
Then, the ring formed and started moving outwards, O_{2} depletion continued and the region of high O_{2} gradients naturally started moving outwards (Figure 3B). At this point, coupled dynamics between the cells and the O_{2} distribution appeared and we observed that the position of the ring closely followed the dynamics of the O_{2} field (Figure 3D), that is it followed a constant concentration of oxygen of 0.25% (Figure 3C).
In the process, three distinct regions were created. Behind the ring, O_{2} was completely depleted and thus no gradient was visible. In front of the ring, the O_{2} concentration remained high with high gradients. Finally, in the ring region, O_{2} was low (<1%) and the gradients were strong. Based on our results in externally imposed gradients, we would thus expect cells to present a positive aerotactic bias mostly in the ring region which is indeed what we observed (Figure 1—figure supplement 5).
Minimal cellular Potts model
Based on these experimental results, we then asked whether this subtle response of Dd cells to complex oxygen environments was sufficient to explain the emergence of a highly stable, selfmaintained collective phenomenon. To do so, we developed cellular Potts models based on experimental observations and tested whether they could reproduce the observed cell dynamics. Briefly, the ingredient underlying the model are as follows (details can be found in the Materials and methods section). First, all cells consume the oxygen that is locally available at a known rate (Torija et al., 2006). Cell activity increases at low O_{2}. Cells respond positively to O_{2} gradients with a modulation of the strength of this aerotaxis based on local O_{2} concentrations, as observed in our microfluidic experiments. Finally, all cells can divide as long as they sit in a high enough O_{2} concentration (chosen at 0.7%) since it was demonstrated that cell division slows down in hypoxic conditions (Schiavo and Bisson, 1989; West et al., 2007). Of note, all parameters were scaled so that both time and length scales in the Potts models are linked to experimental times and lengths (see Materials and methods).
Although this model is based on experimental evidence, some of its parameters are not directly related to easily measurable biological properties. Therefore, we decided to fit our parameters to reproduce as faithfully as possible the results of our microfluidic experiments. Through a trial and error procedure, we managed to reproduce these results qualitatively and quantitatively (Video 4) in terms of collective behavior, cell accumulation, and individual cell behavior (Figure 4—figure supplement 1).
We then applied this model and added O_{2} consumption by cells, with initial conditions mimicking our spot assay and other ingredients mimicking the vertical confinement. We observed the rapid formation and migration of a ring (Figure 4A–B, Video 5). This ring was remarkably similar to that observed in experiments. In particular, we found its speed to be constant after an initial transitory period (Figure 4C, Figure 4—figure supplement 2). This speed was also comparable to experimental ones. Similarly, the morphology of these simulated rings was constant over time with a fixed cell density and width (Figure 4—figure supplement 2). Finally, cell behavior was qualitatively well reproduced by this model (Figure 4—figure supplement 3).
In terms of coupled dynamics between cell density and O_{2} profiles, we found here too that the driving force behind this collective phenomenon was the fact that the ring followed a constant O_{2} concentration (Figure 4D–E).
We then asked what were the key ingredients in the model to trigger this phenomenon, a question we explored by tuning our original Potts model. We started by dividing cell consumption of oxygen by a factor of 3 (Figure 5A) and found that it did not significantly change the ring speed but could change the aspect of cell density in the central region. We then turned our attention to other key elements in the model.
If we turned off cell division in our models, the formation of the ring was mostly unchanged but after a short time, the ring started slowing down and even stopped as cell density was no longer sufficient to reach highly hypoxic conditions (comparing Figure 5A and B). Second, we asked whether the observed and modeled aerokinesis was necessary to reproduce the collective migration. We found that it wasn’t as models ran at different effective temperatures applied to all cells regardless of local O_{2} concentrations all showed qualitatively similar behavior (see for example Figure 5C). Of note though, lower effective temperatures led to less dense rings as fewer cells were able to start in the ring (Figure 5—figure supplement 1). Finally, we found that modulation of aerotaxis by local O_{2} concentrations was essential. Indeed, as we increased the range of O_{2} concentration at which aerotaxis is at play (Figure 5G–H), we found that forming rings became wider and less dense (Figure 5D–E) to the point where no actual ring could be distinguished if aerotaxis was kept constant for all cells (Figure 5I).
These numerical simulations based on cellular Potts models provide a good intuition of the phenomenon and reveal that cell division and aerotactic modulation are the two key ingredients to reproduce the ring of cells. Because of their versatility, they can also be used to make some predictions on the observed phenomenon. Experimentally, we tested two such predictions to demonstrate the relevance of the underlying assumptions.
First, we show in Figure 5B the effect of turning cell division off in the simulated spot. A similar result can be achieved by placing cells in a phosphate buffer medium, lacking nutrients and thus blocking cell division (Kelly et al., 2021). In this situation, at short time scales, a ring of cells started forming and expanding outwards in a similar fashion as in nutritive medium (Figure 5—figure supplement 2). After a few hours, however, the ring started slowing down until it completely stopped and cells started dispersing again. This is in complete agreement with the predictions of the Cellular Potts Model, as one can see by comparing the density kymographs (Figure 5A and Figure 5—figure supplement 2) and firmly demonstrates the importance of cell division in this behavior. At longer time scales though (t>10h), Dd cells started forming aggregates and entering a developmental phase (Video 6). This aggregation is presumably due to the concomitant expression of cell adhesion molecules and, apparition of selforganizing secreted cAMP pulses whose timing agrees with the one reported in classical free cell spot aggregation assays (Gregor et al., 2010). Cellcell adhesion and cAMP signaling are not included in our models or numerical simulations that hence cannot predict the long times in Video 6. However, the timing is well separated from the end of the ring expansion period (t<3.5h). This still demonstrates that the phenomenon observed here is relevant for both the single cell and collective stages of Dd life cycle.
Second, we used these numerical simulations to predict the behavior of cells in more complex environments. One can see the expansion of the ring as a way for each cell to optimize its own resources. This begs the question of what happens when more than one colony is present in the environment, a problem more directly relevant for real life situations. Can the different colonies sense their respective presence and adapt accordingly by migrating preferably away from one another or, on the other hand, will the depletion of oxygen induced by a neighboring colony increase hypoxia on this side and therefore accelerate migration? In this case, what would happen when two rings come in contact? We started exploring this question by simulating two colonies put in close proximity. These simulations predict that the formed rings do not repel each other, instead they tend to rush toward one another and, when they meet, they fuse together to make an elliptical front which then relaxes towards a more circular shape (Video 7). We then performed the corresponding experiment and found very similar behavior (Video 7).
Overall, these results show that the cellular Potts model indeed recapitulates all the major experimental observations with only two key ingredients (cell division and aerotactic modulation). However, they fall short of giving an indepth quantitative description because they rely on many parameters and are not amenable to theoretical analysis per se.
'Go or Grow' hypothesis: a Meanfield approach
In order to complement the methodology of the cellular Potts model, we developed a meanfield approximation of the latter: the cell density $\rho $ is subject to a reactionadvectiondiffusion partial differential equation (PDE):
$C$ is the oxygen concentration, $a\left(C,\nabla C\right)$ corresponds to the aerotactic advection speed and $r\left(C\right)$ to the cell division rate. By assuming radial symmetry in agreement with the experiments, we propose $a\left(C,\nabla C\right)=a\left(C,{\partial}_{r}C\right)={\lambda}_{aero}\left(C\right){\partial}_{r}C$, where ${\lambda}_{aero}\left(C\right)$ is the already mentioned aerotactic strength fitting the microfluidic experiments with an upper O_{2} concentration threshold ${C}_{0}$=0.7% (Figure 5G and Material and Methods) and $r\left(C\right)=\{\begin{array}{l}{r}_{0},if\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}C>{C}_{0}\\ 0,if\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}C<{C}_{0}\end{array}$ is the division rate. When not specified, we use the same threshold ${C}_{0}$ for cell division and aerotaxis as for the cellular Potts model. Below, this assumption is coined as the ‘Go or Grow’ hypothesis. We thereby revisited the ‘Go or Grow’ model for glioma cells (Hatzikirou et al., 2012) with the transition between division and directional motion being mediated by oxygen levels rather than cell density in the mentionned study. Congestion effects such as they may arise in the cellular Potts model or in experiments have been ignored.
Oxygen is subject to a simple diffusionconsumption equation, with $b\left(C\right)$ the consumption rate of oxygen per cell (see Materials and methods):
The results obtained by numerical simulation of this meanfield model are comparable to the ones already obtained by the cellular Potts model: emergence of a high cell density area traveling at constant speed $\sigma \approx 1.0$ µm/min, leaving behind a trail of cells (Figure 6AB).
A general framework for traveling waves in cells undergoing aerotaxis and division
From there onwards, we propose a mathematical framework that investigates general conditions under which collective behavior of cells undergoing cell division and aerotaxis is triggered. The aim was to confirm conclusions already obtained experimentally or through the cellular Potts model and to decipher the contribution of cell division to the collective behavior, while also keeping the framework relatively general such that it may be applied to other types of collective cellular behavior.
We first considered models of the form given by (Equations 1 and 2), independently of the exact shape of the advection term $a\left(C,\nabla C\right)$ or division term $r\left(C\right)$. Because of its relevance for the study of planar front propagation, we studied these models in a planar symmetry ($\rho =\rho \left(t,x\right)$,$C=C\left(t,x\right)$) instead of a radial symmetry ($\rho =\rho \left(t,r\right)$,$C=C\left(t,r\right)$), neglecting thereby any curvature effects. We were interested in the study of a single front moving from left to right. Introducing the front speed $\sigma $, the front corresponds to a stationary solution in the moving frame $z=x\sigma t$, that is, a traveling wave profile, satisfying:
In the general case, the theoretical analysis of such profiles and the determination of the front speeds $\sigma $ seem out of reach due to the coupling with the reactiondiffusion equation on the O_{2} concentration. Nonetheless, it is possible to derive simple relations between the shape of the wave and the speed of propagation. By integrating (Equation 3) over the line, we obtain:
This equation balances the net flux of cells to the far lefthand with the amount of mass created by heterogeneous (oxygendependent) cell division. We illustrated this relationship with the experimental data from Figure 1E. In order to approximate the term $\int r\left(C\left(z\right)\right)\rho \left(z\right)dz$, we used an observation that we made through numerical simulations: cell division stops roughly at half of the peak, meaning that cells left to the peak do not divide, while cells right to the peak continue dividing (see Figure 4E and Figure 6— figure supplement 1). Therefore, we approximated by a rectangle method $\int r\left(C\left(z\right)\right)\rho \left(z\right)dz={\widehat{\rho}r}_{0}L/2$, where $L$ is the length spanned by the ring and $\widehat{\rho}$ is the average cell density in the ring. As there is supposedly no advection $a\left(C,{\partial}_{z}C\right)\rho =0$ at $z=\infty $ this yields the approximation $\sigma \approx \frac{{r}_{0}L\widehat{\rho}}{2\rho \left(\infty \right)}$. Quantitatively, we assume L to be on the order of 300µm (Figure 1E) and $\frac{\widehat{\rho}}{\rho \left(\infty \right)}$, the ratio between cell densities in the ring and in the bulk of cells, to be on the order of 4 (Figure 1E). This yields an estimate of the wave speed, based solely on the shape of the cell density profile, of $\sigma \approx 0.9$ μm/min.
Mathematical analysis of the ‘Go or Grow’ hypothesis
The difficulty to study (Equation 3) analytically led us to propose a simpler version of the meanfield model that recapitulates the two key ingredients, cell division and aerotaxis, in an original way. Although it deviates from the reference Potts model in the details, it has the advantage of being analytically solvable. Cells have two distinct behaviors, depending on the O_{2} concentration. Below a certain threshold ${C}_{0}$ cells move preferentially upward the oxygen gradient (go), with constant advection speed ${a}_{0}$, but they cannot divide. Above the same threshold they divide (grow) and move randomly without directional bias. This model may be considered as a strong simplification of (Equation 1), here restricted to the onedimensional space, where:
The coupling between (Equations 1 and 2) then goes merely through the location of the oxygen threshold ${C}_{0}$. This elementary ‘Go or Grow’ model was meant to 1 demonstrate that its simple ingredients suffice to trigger a collective motion and 2 determine the relative contributions of cell division and aerotaxis on the speed of the ring in a general framework.
Interestingly enough, in this case (Equation 3) admits explicit traveling wave solutions (see more details in the Materials and methods section). Moreover, an explicit formula for the wave speed was obtained (Figure 6C and Materials and methods for a detailed derivation):
To the best of our knowledge, this analytical formula is new and captures basic features of a wave under a single selfgenerated gradient. It is remarkable that Formula (Equation 6) does not depend on the dynamics of oxygen consumption and diffusion. Furthermore, Formula (Equation 6) presents a dichotomy according to the relative size of aerotaxis strength ${a}_{0}$ and the quantity $\sqrt{{r}_{0}D}$: in the case of smallbias (i.e. ${a}_{0}\le \sqrt{{r}_{0}D}$), the wave speed $\sigma $ is independent of aerotaxis and coincides with the socalled Fisher’s wave speed $2\sqrt{{r}_{0}D}$. This speed is related to the FisherKPP equation (Aronson and Weinberger, 1975; Fisher, 1937; Kolmogorov et al., 1937), which describes front propagation under the combined effects of diffusion and growth (without advection). However, in the case of largebias (i.e., $a}_{0}>\sqrt{{r}_{0}D$), aerotaxis is strong enough to contribute to the speed and the wave speed increases $\sigma >2\sqrt{{r}_{0}D}$.
Based on these observations, we propose the fraction $\phi =2\sqrt{{r}_{0}D}/\sigma $ as a measure of the relative contribution of cell division and diffusion to the overall wave speed. Indeed, when aerotaxis is absent (or as in the smallbias case not contributing to the wave speed), the value of $\phi $ is 1 and the wave is driven by cell division and unbiased random motion, that is, a reactiondiffusion wave. In the largebias case, $1/\phi $ describes how much faster the wave travels, compared to if it were only driven by diffusion and division. We illustrated the behavior of $\phi $ with a heatmap (Figure 6D) as a function of the parameters ${a}_{0}$ and $\mathrm{l}\mathrm{n}\left(2\right)/{r}_{0}$ (the doubling time of the cell population), the diffusion coefficient being fixed to its experimental value D=30 µm^{2}/min.
We confront this reasoning with the experimental data: as a rough approximation with ${a}_{0}=1$ µm/min in experiments, assuming a doubling time of 8 hr for Dd cells, ${r}_{0}=ln2/480$ min^{−1}, we are clearly in the case of large bias ($\sqrt{{r}_{0}D}=0.2$ µm/min) and (Equation 4) yields $\sigma =1.04$ µm/min while the fraction $\phi =40\text{\%}$. The wave travels 2.5 times faster than a wave merely driven by cell division, showing that in this case the dominant ingredient to set the wave speed is aerotaxis. Still, our results can similarly be applied to other systems in which this balance could be different. Finally, note that the density profile of the model (Figure 6E) does not present a sharp front peak as in experiments (Figure 1D,E), Potts simulations (Figure 3B,C) or in the mean field model (Figure 6A,B). We will show below that it can be slightly modified to change the profile of the fronts while keeping the analytical results relevant thus describing a whole class of systems (Figure 6E and Figure 6—figure supplement 1).
Inside dynamics of the wave front
The wave speed of the elementary ‘Go or Grow’ model coincides with Fisher’s speed, that is $\sigma =2\sqrt{{r}_{0}D}$, in the regime of small bias (${a}_{0}\u2a7d\sqrt{{r}_{0}D}$). This is the signature of a pulled wave, meaning that the propagation is driven by the division and motion of cells at the edge of the front, with negligible contribution from the bulk, and little diversity in the expanding population. In contrast, when the bias is large ($a}_{0}>\sqrt{{r}_{0}D$) then the wave speed in (Equation 6) is greater than Fisher's speed. This is the signature of a pushed wave, meaning that there is a significant contribution from the bulk to the net propagation, with an expanding population maintaining diversity across expansion, see Birzu et al., 2018; Stokes, 1976 for insights about the dichotomy between pulled and pushed waves. In particular, it was conjectured that the ratio $\phi =2\sqrt{{r}_{0}D}/\sigma $ proposed above controls the transitions between different regimes of diversity loss in a wide class of reactiondiffusion models (Birzu et al., 2018; Birzu et al., 2019).
In order to explore this dichotomy between pulled and pushed waves, we used the framework of neutral labeling (Roques et al., 2012) in the context of PDE models. We colored fractions of the density profile during wave propagation to mimic labeling of cells with two colors. Then, we followed numerically the dynamics of these fractions, and quantified the mixing of the two colors. Our results were in perfect agreement with (Roques et al., 2012), extending their results beyond classical reactiondiffusion equations to equations which also include advection (see Materials and methods). In the case of large bias (Figure 7A–C), the wave is pushed and the profile is a perfect mixture of blue and yellow cells at long times. Contrarily, the wave is pulled in the regime of small bias: only cells that were already initially in the front, here colored in blue (Figure 7B–D), are conserved in the front, whilst yellow cells at the back cannot catch up with the front.
In the absence of associated experimental data, we explored the cellular Potts model with such neutral labeling. The results were in agreement with the PDE simulations (Figure 7—figure supplement 1) showing a clear, rapid mixing of the two cell populations under the propagation of a pushed wave in the regime of experimental parameter values.
Robustness of the conclusions to structural variations of the model
We voluntarily defined our elementary ‘Go or Grow model’ as a rough simplification of our original meanfield model in order to keep it solvable and extract a general formula for the front speed and an analysis of the relative contribution of diffusion/division and aerotaxis in that respect. However, many experimental systems will not conform to the hypothesis underlying this model (in particular the shapes of the aerotactic response and cell division modulation). In order to investigate whether the conclusions drawn from the elementary ‘Go or Grow’ model extend to more general situations, we decided to submit it to structural variations and check if the results obtained above still held. First, we made the hypothesis of a second oxygen threshold $C}_{0}^{\prime}<{C}_{0$, below which cells are not sensitive to gradients any longer (Figure 6F). In the general case, we were not able to do a thorough analysis of this model, but through numerical exploration we found that the propagation speed remained close to the value given by formula (Equation 6) (at most 15% of relative difference in a relevant range of parameters, Figure 6C). Intuitively, the main contribution to the collective speed is the strong bias inside the highdensity area at intermediate levels of O_{2}, whereas cells at levels below the second threshold ${C}_{0}\text{'}$, where the dynamics of both models diverge, do not contribute much to the collective speed. We also noticed that cell density profiles (Figure 6F) were much closer to experimental observations and results obtained through the cellular Potts model or the original meanfield approach. Moreover, the wave speed is no longer independent of the oxygen dynamics. In the Materials and methods section, we pushed further the analysis with a specific form of oxygen consumption and developed a specific case of such a ‘Go or Grow’ model with a second threshold, where we were able to conduct its complete analysis. Figure 7—figure supplement 2 shows that the conclusions concerning the contribution of growth to the wave speed are robust. Finally, we show on this modified ‘Go or Grow’ model that our conclusions regarding how the behavior can switch from a pulled to a pushed wave remain true as well (Figure 7—figure supplement 3) demonstrating that our results can be generalized to a variety of different systems showing the propagation of a front in response to a single selfgenerated gradient.
To go beyond this first variation with two oxygen thresholds, we also investigated the influence of the shape of the aerotactic response such as linear or logarithmic gradient sensitivity. Figure 6—figure supplement 1 shows the qualitative outcomes of these different models. This numerical exploration indicates that a wide combination of the two key ingredients, aerotaxis and cell division, can drive the propagation of a stable wave with various density profiles. Cell division at the edge yields a net flux of cells, backward in the moving frame, that sustains the wave propagation in the long term, but may have a relatively small contribution to the wave speed.
Discussion
Overall, our results demonstrate the ability of Dd cells to respond to hypoxia through both aerotactic and aerokinetic responses. Both of these behaviors could be very important to help Dd cells to navigate complex hypoxic environments they encounter in the soil. In addition, our results are a confirmation of the ability of selfgenerated gradients to serve as very robust, longlasting directional cues in environmental navigation, a property which has recently emerged in a variety of systems (Cremer et al., 2019; Tweedy and Insall, 2020). Finally, our work goes beyond theses results as it demonstrates that oxygen can play the role of the attractant in selfgenerated gradients therefore potentially extending the physiological relevance of the use of such cues in collective migration.
In addition, although our experimental results were obtained on simple, 2d experiments, our findings can generalize to more complex cases. The fact that the dense front of cells follows a constant oxygen concentration (Figure 3E, Figure 4E) provides a hint that any situation in which cell density is locally high enough to trigger hypoxic conditions will also lead to a similar behavior. Then, depending on the dimensionality of the system, its architecture and the position of possible oxygen sources, we hypothesize that a similar front will develop and follow isoconcentration lines. Indeed, the original experiments of Adler (Adler, 1966) and more recent developments (Cremer et al., 2019; Fu et al., 2018; Saragosti et al., 2011) on bacteria demonstrated that similar ingredients as the ones presented here can lead to front propagation in both 1d and 2d situations. Similarly, using an under agarose assay, it was demonstrated that selfgenerated gradients of degraded folate induce a group migration of cells in bands (in 1D) or rings (in 2D spots) up to 4 mm (Tweedy and Insall, 2020). Beyond the dimensionality of the system, it was also shown that selfgenerated gradients allow cells to solve mazes by locally degrading an attractant that has a source at the exit of the maze (Tweedy et al., 2020). Our results are in total agreement with these past examples. To further show the generality of the underlying principles, we ran some 3D Potts simulations using a qualitative version of our model. Briefly, we show that in three dimensions, if oxygen is provided on all sides, a spherical front of cells starts moving outwards (Video 8). However, if we assumed that the bottom of the space was completely deprived of oxygen (i.e. a symmetry breaking situation that can be encountered in various physiological situations), this front was migrating upwards only in a halfspherical shape (Video 8). Our 2d results can therefore be extended to any other situations and they show that the key to proper steering are high enough cell densities and the creation of robust selfgenerated gradients.
While aerotaxis is well established for bacteria, its role is often invoked in multicellular organisms to explain various processes in development or cancer progression but very few in vitro studies were conducted to prove it is an efficient and operating mechanism or to understand the molecular mechanisms at play during aerotaxis. Deygas et al. showed that confined epithelial colonies may trigger a selfgenerated O_{2} gradient and an aerotactic indirect response through a secondary ROS selfgenerated gradient (Deygas et al., 2018). Gilkes et al. showed that hypoxia enhances MDAMB231 breast cancer cell motility through an increased activity of HIFs (Gilkes et al., 2014). HIFs activate transcription of the Rho family member RHOA and Rho kinase 1 (ROCK1) genes, leading to cytoskeletal changes, focal adhesion formation and actomyosin contractions that underlie the invasive cancer cell phenotype. This study suggests a role for aerotaxis in tumor escape, but it only demonstrates aerokinesis as O_{2} gradients were not imposed to probe a directed migration toward O_{2}. Using a microfluidic device, the same cancer cell line was submitted to various oxygen levels as well as oxygen gradients (Koens et al., 2020) but the observed aerotactic response was not clear.
By contrast, the experimental results presented here with Dd show a strong response to hypoxia. Within 15 min, cells exhibit an aerokinetic and aerotactic response when exposed to externally imposed O_{2} gradients (Figure 2). Selfgenerated O_{2} gradients are produced within 20 min (Figure 3 and Figure 1—figure supplement 2). But this cellular response is within the equilibration time of the oxygen distribution (Figure 2—figure supplement 6). Hence we can consider the cellular response as almost instantaneous with Dd. The difference with previously studied cells is probably due to the extreme plasticity of the rapidly moving amoeboid cells (Dd) and their almost adhesion independent migration mechanism (Friedl et al., 2001) while mesenchymal cancer cells move slower by coordinating cytoskeleton forces and focal adhesion (Palecek et al., 1997).
The quick response of Dd in directed migration assays has been largely exploited to decipher the molecular mechanisms at play during chemotaxis (Nakajima et al., 2014). The molecular mechanisms used for O_{2} sensing and its transduction into cellular response are for the moment unknown but we can expect that the O_{2} molecular sensors modulate cytoskeleton organization, particularly localized actin polymerization/depolymerization through some of the molecular components involved in classical chemotaxis toward folate or cAMP (Pan et al., 2016; van Haastert et al., 2017). However, new and unexpected mechanisms cannot be excluded.
The finding that migrating cells can influence the direction of their own migration by building chemoattractant gradients is not new. Several species of bacteria can move preferentially toward oxygen or nutrient as reported by Adler, 1966. However, this mechanism was only recently reported in eukaryotes (Stuelten, 2017): melanoma cells that break down lysophosphatidic acid (LPA) and generate a LPA gradient in the vicinity of the tumor (MuinonenMartin et al., 2014), Dd colonies that generate folate gradients (Tweedy et al., 2016) or for the migration of the zebrafish lateral line primordium through a selfgenerated chemokine gradient (Donà et al., 2013; Venkiteswaran et al., 2013). The dispersal of melanoma cells is particularly instructive. The stroma surrounding the tumor acts as a source of LPA. The tumor cells act as a sink for LPA. As long as LPA is present in the environment a steady wave of migrating melanoma cells propagates away from the initial tumor over long distances and long time periods.
The selfgenerated LPA (melanoma) and folate (Dd) gradients were modeled with a simple numerical model that was able to predict the steady wave. In particular, it predicted an invasive front where cells are exposed to a steep chemoattractant gradient, followed by a ‘trailing end’ where the gradient is shallow and fewer cells migrate with poor directionality (Tweedy and Insall, 2020). It also predicted that the wave may have a less marked front, and/or a smaller speed, or even vanishes if the cell density was too low due to insufficient chemoattractant removal. All these features are surprisingly similar to our experimental measurements of cell density and O_{2} profiles (Figure 1E, Figure 3C). The atmospheric O_{2} that diffuses through the culture medium and eventually the plastic surfaces is the chemoattractant. The O_{2} consumption triggers hypoxia that in turn generate an aerotactic response toward O_{2} in a very narrow range of O_{2} concentrations (0–1.5%) (Figure 3C). The exact value of the lower O_{2} threshold value will deserve future investigations. The exact nature of the cellular response at these extremely low O_{2} levels, and in a very shallow gradient, also has yet to be clarified.
Our different models unveil a set of basic assumptions which are sufficient for collective motion of cells without cellcell interactions (attractive or otherwise), in contrast with (Sandonà et al., 1995). Cell growth is necessary to produce a longstanding wave without any damping effect. However, it may not be the main contribution in the wave speed, depending on the relative ratio between directional motion (the bias ${a}_{0}$), and reactiondiffusion (the Fisher halfspeed $\sqrt{{r}_{0}D}$). In the case where the former is greater than the latter, the wave is due to the combination of growth and directional motion and it is pushed. This result differs particularly from the FisherKPP equation with constant advection (meaning with uniform migration and division) where the wave speed is ${a}_{0}+2\sqrt{{r}_{0}D}$ and the wave is pulled. In the experiments under study, we estimate directional motion to contribute the most to the cell speed, ruling out the possibility of seeing a pulled wave driven by cell division and diffusion at the edge of the front.
In conclusion, we demonstrate the remarkable stability of collective motion driven by selfgenerated gradients through depletion of oxygen. Through coupled dynamics, these gradients give rise to long lasting, communicationfree migrations of entire colonies of cells which are important both from ecological and developmental points of view. In the case presented here where oxygen plays the role of the depleted attractant, this could prove to be a very general mechanism as oxygen is ubiquitous and always consumed by cells.
Materials and methods
Cell preparation
Request a detailed protocolThe AX2 cell line was used and cultured in HL5 media (Formedium, Norfolk, UK) at 22°C with shaking at 180 rpm for oxygenation (Sussman, 1987). Exponentially growing cells were harvested, counted to adjust cell density to the desired one, typically 2000 cells/μL.
Observations and analysis of selfgenerated aerotaxis by cell confinement (spot assay)
Request a detailed protocolFor the spot assay, 1 µL of cell suspension containing 1000–8000 cells (typically 2000) was carefully deposited on a dry surface using a 1 µL syringe (Hamilton, Reno, NV, USA). The dry surface was either the nunclon treated surface of Nunc six wells polystyrene plates for usual experiments (ThermoFisher Scientific, Waltham, MA, USA) or polydimethylsiloxane surface (PDMS, Sylgard 184, Dow Corning, Midland, MI, USA) for experiments on oxygensensing films. The drop was incubated for 5 to 7 min in humid atmosphere at 22°C before gently adding 2 ml of HL5 medium without detaching the spotted cells forming a microcolony. A 14 mm or 18 mm diameter round glass coverslip cleaned in ethanol, thoroughly rinsed in HL5 was kept wet and deposited on top of it. In some experiments, fluorescein FITC at 16 µM was added to the HL5 medium and confocal slices were taken, showing that confined Dictyostelium cells were not compressed by the coverglass but separated from it by a layer of medium of about 50 µm (Figure 1—figure supplement 1).
The outward spreading of the Dictyostelium microcolony was observed at 22°C in transmission with three types of microscope: (i) a TE2000E inverted microscope (Nikon, Tokyo, Japan) equipped with motorized stage, a 4x Plan Fluor objective lens (Nikon) and a Zyla camera (Andor, Belfast, Northern Ireland) using brightfield for most of the experiments lasting 16 hr (Figure 1A), (ii) a binocular MZ16 (Leica, Wetzlar, Germany) equipped with a TL3000 Ergo transmitted light base (Leica) operated in the onesided darkfield illumination mode and a LC/DMC camera (Leica) for experiments over days (Figure 1B) and finally (iii) a confocal microscope (Leica SP5) with a 10x objective lens for a few larger magnification experiments (Figure 1C).
For computing densities, cell positions were determined using the builtin Find Maxima plugin in ImageJ (National Institutes of Health, Bethesda, MD, USA) through a custom made routine. Data analysis and plotting was performed in Matlab (Mathworks, Natick, MA, USA). For density profiles (Figure 1E) and kymographs (Figure 1D), the center of the colony was defined as the center of mass of all cells detected at all times. Cell positions were then turned into radial coordinates and cells were counted within concentric crown regions. Densities were calculated by dividing this count by the area of each crown.
Density profiles such as the ones showed in Figure 1E were treated to automatically extract the position, width and density of a ring in various experiments and at various time points. Density profiles were first stripped of values lower than 500 cell/mm^{2} in order to avoid asymmetric baselines behind and in front of the ring. Resulting profiles were then fitted in Matlab by a Gaussian function with a nonzero baseline. The nonzero baseline corresponds to density in the bulk, the maximum of the Gaussian gives ring position, its height added to the nonzero baseline gives the cell density in the ring and its width the width of the ring.
Cell tracking, diffusion coefficients, and aerotactic biases
Request a detailed protocolAfter retrieving cells’ positions with optimized ImageJ macros based on Find Maxima, the individual trajectories were reconstructed with a squareddisplacement minimization algorithm (http://site.physics.georgetown.edu/matlab/). Data were analysed using inhouse Matlab programs. Timelaspse microscopy experiments devoted to cell tracking in the spot assay experiments was acquired at a high frame rate (1 frame every 15 s) (Figure 1—figure supplements 4–5) due to the very high cell density in the ring region (up to 2000 cells/mm^{2}, Figures 1E and 3C). For the microfluidic experiments, as cells were plated at a lower density (less than 200 cells/mm^{2}), 1 min time intervals was used to track cell trajectories (Figure 2C). In order to highlight aerotactic biases, cells displacements over various time lags dt (dt up to 60 min) were projected in the radial direction for spot assays and in the gradient direction X for microfluidic experiments and eventually divided by dt to obtain velocity biases. Individual biases were then averaged within bins of equal distance (Figure 2D, Figure 1—figure supplement 5, Figure 4—figure supplement 1). Individual effective cell diffusion constants were measured as the square of their displacement over their entire trajectory divided by the trajectory time length and divided by 4. These measurements were then similarly averaged over bins (Figure 1—figure supplement 4).
Microfluidicbased oxygen gradient generator: design, fabrication, and cell injection
Request a detailed protocolA schematic of the present doublelayer microfluidic device is shown in Figure 2A. It is made of several layers of PDMS mounted on a bottom glass coverslip. The overall diameter D of the microfluidic device is 27 mm and the overall thickness H is 4 mm. Three parallel media channels are positioned for cell culture, and two gas channels are positioned at a height H_{g}=0.5 mm above the media channels to allow gas exchange between the channels during cell culture. The horizontal distance between the two gas channels narrows stepbystep (2 mm, 1 mm, and 0.5 mm gaps, respectively), thus yielding to generate different gradients of oxygen concentration along the media channels simultaneously. All channels are 125 µm high and 2 mm wide, and therefore, the media and gas channels are separated by a PDMS wall of 375 µm thickness. A polycarbonate (PC) film (26 mm in diameter and 0.5 mm thickness) is embedded inside the device at a height Hf=1 mm from the bottom coverslip to prevent oxygen diffusion from the atmosphere. The cartesian coordinate origin was set at the center of the media channel (median axis), and the x and ydirections were defined as parallel to media and gas channels respectively (Figure 2A). The zdirection was set to the vertical direction from the top of the bottom coverslip.
The manufacturing steps are as follow. The media channel and gas channels were drawn with AutoCAD (Autodesk, Mill Valley, CA, USA) and replicated in SU8 photoresist using classical photolithography techniques. These SU8 molds were silanized to make it nonadherent and reusable. PDMS was mixed at a 10:1 ratio of base:curing agent, poured over each mold to a thickness of H_{g}, and cured in an oven at 60°C for more than four hours. On top of the cured PDMS layer of the gas channels, the abovementioned PC film with 3 mm port holes punched at the locations of the media and gas channel ports was positioned. Additional PDMS was then poured over the PC film until the total PDMS layer became 3.5 mm thick, then the PDMS layer was cured in an oven at 60°C overnight. The PDMS layers of the media and gas channel patterns were peeled off the silicon wafers and cut into 27 mm diameter circles. The PDMS layer with the gas channel pattern was punched to form inlets and outlets 2 mm in diameter to allow the infusion of gas mixtures. The channelpatterned surface of the PDMS layer with the gas channels and the top surface of the other PDMS layer with the media channels were plasma treated (PDC001HP; Harrick Plasma Inc, Ithaca, NY, USA) to bond with each other. After incubating the bonded PDMS mold overnight in an oven at 60°C, 2 mm diameter inlets were punched to allow access to the media channels, respectively. Finally, the channelpatterned side of the PDMS mold and a 35 mmdiameter glass bottom dish with or without covered by an oxygen sensing film were plasma treated and bonded each other.
Measurements of aerokinesis in homogenous environments were performed using a homemade glassduralumin environmental chamber to perform random motility assays (d’Alessandro et al., 2018; Golé et al., 2011).
Dictyostelium cells were seeded in the media channels at density of 2x10^{6} cells/ml, and the cell culture medium was filled in the glass bottom dish up to the height covering the PDMS mold. Cells were allowed to adhere to the bottom surface (bare glass or coverglass covered with a sensing film) for 15 min.
Gas control and injection
Request a detailed protocolWe used a controlled oxygen concentration for three types of experiments: (i) to calibrate oxygen sensing films (see below), (ii) to create the oxygen gradients within microfluidic devices (see below) or (iii) to insure a pure hypoxic condition (pure N_{2}) at the end of the spot assay experiment. The gas mixture (0% to 21% O_{2} in N_{2}) was prepared in a gas mixer (Okolab 2GFMIXER to mix compressed AIR with 100% N2 or HORIBA STEC MU3405, Kyoto, Japan to mix pure O2 and N2) by mixing pure O_{2} (or air) and pure N_{2}. Free sensing films for calibration (i) or for the spot assay (iii) were placed inside 6wells plates and the multiwells were placed in an environmental chamber fitting our microscope stage (H301Kframe, Okolab, Pozzuoli, Italia). Gas was injected at about 500 mL/min in this chamber. Eventually, multiwells were drilled to a diameter of 25 mm and the sensing films were glued with a silicone adhesive on the plate bottom to reduce the background noise from fluorescence. For microfluidic experiments, the tubes from the mixer were connected to the gas channels and gas was injected at a controlled flowrate (between 60 and 180 mL/min) into the device.
Oxygensensing film preparation
Request a detailed protocolOxygensensing films were prepared by inserting the luminescent O_{2}sensitive dye 5,10,15,20Tetrakis(2,3,4,5,6pentafluorophenyl)porphyrinPt(II) (PtTFPP, PorLab, PorphyrinLaboratories, Scharbeutz, Germany) in a 4:1 PDMS:curing agent thin layer spincoated on 30 mm to 35 mm rounded coverglasses. Briefly, 17 mg of PtTFPP was dissolved in 5 mL chloroform and thoroughly mixed with 2.8 mg of PDMS and 0.7 mg curing agent. The mixture was degassed in a vacuum chamber for 5 hr. About 0.5 mL to 1 mL of this solution was spread on the coverglass and spincoated for 2 min at 500 rpm with a final speed of 2000 rpm during 10 s to flatten the edge bead. Chloroform was allowed to evaporate overnight while the polymer cured at 60°C. The final PtTFPP sensor film had a dye concentration of 4 mmol/L and was 25 μm thick. This thickness was measured using a ContourGTK 3D Optical microscope (Bruker, Billerica, MA, USA) after removing a piece of film with a surgical blade. Sensing films were stored in dark. They were used to measure the oxygen concentration in selfgenerated O_{2} gradients (spot assay) and for microfluidic experiments with controlled O_{2} gradients.
Fluorescence microscopy for oxygen measurements
Request a detailed protocolFluorescence images of O_{2}sensing films (either for film calibration or for in situ oxygen measurements in the spot assay or in microfluidic devices) were taken with two inverted epifluorescence microscopes: (i) a TE2000E inverted microscope (Nikon) equipped with motorized stage, a 4x Plan Fluor objective lens (Nikon), a XCite Series 120PC illumination lamp, a TRITC bandpass filter cube and a Zyla camera (Andor) (Figure 3Aii, Figure 2—figure supplement 3), (ii) a IX83 inverted microscope (Olympus, Tokyo, Japan) equipped with a motorized stage, a UPlanSApo 4x objective lens (Olympus), a UHGLGPS lamp (Olympus), a RFP bandpass filter and a Zyla camera (Andor). This second microscope was used for mosaic imaging, in order to scan the whole dimension of the three media channels (about 1 cm in length) thanks to the dedicated imaging software cellSens (Olympus) (Figure 2—figure supplement 1).
Oxygensensing film calibration
Request a detailed protocolCalibration was carried out with the sensing films in air, in water or in HL5 culture medium. We applied gas concentration ramps with steps of 5 min for calibration in air (time to exchange fully the gas composition of the chamber and tubes, as O_{2} almost instantaneously diffuses within the 25 µm thick sensing film) and with much longer steps (i.e. 2–4 hr) for calibrations in liquid. There is indeed an additional diffusion time in the PDMS intermediate layer of our microfluidic devices or in the medium height of a Petri dish: typically, a few minutes for a 0.5mmthick PDMS layer and 1h30 for a 2.7mmthick liquid layer in a dish.
Timelapse fluorescence images we recorded and signal intensity I was measured in ROIs of typically 64x64 pixels in various positions of the image, especially along a line scanning the middle of the image (Figure 2—figure supplement 1B, Figure 2—figure supplement 3A). The response of the sensing film in the presence of oxygen can be modeled by a linear SternVolmer relationship:
where C is the oxygen concentration expressed as a percentage of oxygen in the injected gas phase (nearly 21% for atmospheric conditions), I_{0} is the reference intensity in the absence of oxygen, B_{g} is the background intensity independent from the oxygen sensitive signal of the PtTFPP molecules and K is the SternVolmer constant used as an indicator of the sensing film sensitivity.
Notice that the background is usually not included in the SternVolmer relationship but a representative background image (O_{2} independent) is subtracted prior to intensity measurements (Nock et al., 2008; Thomas et al., 2009). This O_{2}independent background value can be deduced from the fluorescence of a plain PDMS film prepared in the same conditions than the sensing film but devoid of PtTFPP molecules (i.e. a standard) (Thomas et al., 2009). We tested that procedure that is basically working but we choose to include the background as a fitting parameter because illumination conditions may change between the sample and the standard (especially the focus plane that affects the focused height of autofluorescent medium above the surface). The slight changes in thickness and PtTFPP composition of sensing films at the large spatial scales we are interested here (3–6 mm wide images, Figure 2—figure supplement 1A–D, Figure 2—figure supplement 3A–C) are another source of heterogeneity especially for the K value. For those reasons, we apply the SternVolmer relation with K and B_{g} as a fitting parameters in many different small regions of interest (ROI) of the surface. For each ROI, we found that the measured intensities follow perfectly the Stern–Volmer relation (i.e. C linearly increases with the SternVolmer parameter (I_{0}B_{g})/(I(C)B_{g})−1, Figure 2—figure supplement 2B) and that K and B_{g} are clearly uncorrelated, in particular B_{g} depends on the illumination pattern but not K.
The illumination pattern is clearly visible on fluorescence images at 21% O_{2}. For instance, the large field of view of Figure 2—figure supplement 1B reconstituted by the multiarea module of the microscope displays an up and down landscape in the 21% intensity (Figure 2—figure supplement 2C) due illumination changes in the periphery of each overlapped area but also due to the slightly different fluorescence in the gap region of the microfluidic device and especially at its interfaces. The single large field of view of Figure 2—figure supplement 3A taken with our second microscope setup displays a dome shaped pattern with 20% intensity difference between the center and borders (Figure 2—figure supplement 3D). The background value B_{g} is very correlated with the 21% signal variations in both imaging configurations (Figure 2—figure supplement 2C, Figure 2—figure supplement 3D) and we can for each experiment calibrate a linear relationship between between B_{g} and I(21%) (Figure 2—figure supplement 2D, Figure 2—figure supplement 3F).
As already stated, the background is due to the autofluorescence of the various media (PDMS layer and coverglass for the sensing film, bottom plastic dish if any, and surrounding fluid). For Figure 2—figure supplement 1, the microfluidic device was filled with pure water and for Figure 2—figure supplement 3, the calibration was performed in the autofluorescent HL5 medium before spotting the cells. Sometimes the sensing film was placed in a nondrilled well and in that case the strong autofluorescence of the plastic bottom of the plate becomes a major source of background (not shown). Finally, B_{g} also includes the read noise RN of the camera which is a constant independent of the light output or exposure time. For Figure 2—figure supplement 2, we measured RN=108 A.U. and hence a light background B_{g}* = B_{g}RN=15 A.U. which is half the ‘true oxygen dependent signal’ at 21% O_{2}, I(21%)RN=30 A.U. The maximum deviation of B_{g} from the linear fit in Figure 2—figure supplement 5D is about 1.5 A.U. Hence a relative error 1.5/15=10% for B_{g}* will be taken in the following. For Figure 2—figure supplement 3, we measured RN=100 A.U. and at the top of the bell curve, B_{g}*=1400–100=1300 A.U. while I(21%)*=1600–100=1500 A.U. (Figure 2—figure supplement 3D). Hence, the background is nearly 87% of the signal due to the HL5 autofluorescence. Nevertheless, the maximum relative deviation from the linear fit is smaller at about 25/1300≈2% of B_{g}*. All these values will be used for the error analysis of the oxygen profiles below.
For uncovered (free) sensing films the sensibility K ranges between 3 and 5 %^{−1} and is very constant, weakly dependent on the illumination pattern (Figure 2—figure supplement 2D and blue points in Figure 2—figure supplement 3E). When films are covered with a coverglass, the fluorescence under hypoxic conditions (0%) increases significantly on the covered region (Figure 2—figure supplement 3C) but not at 21% (Figure 2—figure supplement 3B). As a results, K, which is proportional to this ratio, increases significantly (Figure 2—figure supplement 3E) but not B_{g} and I(21%) (Figure 2—figure supplement 3D) confirming that K and B_{g} are independent. This increase in K is probably due to a local temperature increase: the coverglass adsorbs more heat from the microscope illumination light and this heat is difficult to evacuate due to the confinement. In principle, for the spot assay experiment, it would be necessary to perform an independent calibration with the SternVolmer relation in the covered situation. However, this is difficult due to the very long time required to equilibrate the oxygen level under the confinement far from the coverglass boundary (this is why we started the SternVolmer fit at ROI_{7} in Figure 2—figure supplement 3B,D,E). A too long procedure causes other problems such as medium evaporation, stage or focus drift… To avoid that, we decided to apply the protocol described in the image analysis pipeline of Figure 2—figure supplement 4. First, we perform a gas calibration ramp and do a SternVolmer analysis in various points of an uncovered sensing film (the subscript U is for uncovered) where I_{0U} is reliable in order to get the linear background relation B_{gU}=α I(21%)_{U} +β and to measure the measure the ratio R=I_{0U}/I(21%)_{U}. Reliable means here any point if gas mixture is applied uniformly when calibrating in a dish or just underneath the gas channel in microfluidic devices. Second, we choose the reference fluorescence image I(21%) immediately before starting any experiment (i.e. just after covering the spot, or just before applying the gradient in the gas channels). From that image, we build a B_{g} image as α I(21%)+β (as B_{g} is the same for uncovered and covered case) and eventually we build a reconstituted I_{0} image as R I(21%). Finally, we subtract and divide images with the ‘Image calculator’ of ImageJ (i.e. pixel by pixel) following the SternVolmer model and hence get a Kvalue image map and subsequently an oxygen map (Figure 2—figure supplement 4C–D). This enables to correct nonhomogeneous illumination conditions or nonhomogeneous sensing film properties as well to quickly estimate error bars on the oxygen map from the estimated errors on B_{g}, I(21%) and I_{0} detailed below.
We already discussed the error on the background. In principle, I(21%) is a reference image (hence error free), however as an experiment (especially the spot assay) may run overnight we need somewhere to evaluate drift in the absolute intensity for instance by running an overnight timelapse experiment with a sensing film under ambient gas conditions. This error is added on I(21%) and was estimated as 2% of the true fluorescence signal corrected from read noise I(21%)RN. Error on I_{0} could be much larger. As the intensity I(C) is strongly nonlinearly increasing with the oxygen concentration C, if we measure an intensity I_{0}* corresponding to a residual small oxygen level C_{0}*, we need to correct the true I_{0} value using the relation ${I}_{0}{B}_{g}\approx {I}_{0}\approx {{I}_{0}}^{}\left(1+K{{C}_{0}}^{}\right)$. Due to oxygen leakage along tubes and within our environmental chamber, when applying 100% N_{2}, we measured a residual ${{C}_{0}}^{}\approx 0.15$ O_{2} in a culture medium dish using a bare fiber oxygen sensor coupled to its commercial oxymeter (Firesting, Pyroscience, Aachen, Germany). Hence with a typical K=5 value, we obtain a very large discrepancy between the measured fluorescence I_{0}* and the ideal one: ${I}_{0}\approx {{1.75I}_{0}}^{}$, but finally this discrepancy is not really dramatic on the measured error again due to the nonlinearity.
The effect of these different error sources on the measured oxygen map is presented in Figure 2—figure supplement 1E and Figure 2—figure supplement 3H for a typical microfluidic and spot experiments. Even if we make a 1.75 error on I_{0}, this has little effect on the profiles except in the very hypoxic region when C<0.25% where the error exceed 50%. But even in the region around C=1%, the error is less than 10%. The error on I(21%) on the other hand has a significant effect on the high oxygen regions but less on the hypoxic regions. Finally, the background error is relatively visible in the intermediate oxygen concentration region (very visible on the side of the spot in Figure 2—figure supplement 3H, but also to some extend around the median axis at C~10% in the microfluidic experiment, Figure 2—figure supplement 1E). Finally, we defined error bars with (maxmin)/2 values of the calculated C when exploring the estimated errors discussed above. In the 0.5–1.5% region were we observe most of the interesting aerotactic behaviors with Dictyostelium cells, the precision on the oxygen concentration ΔC/C is less than 0.3. For the purpose of this paper, we can conclude that aerotaxis and aerokinesis occurs undoubtedly between C=0% and C=2% (Figure 2).
Numerical simulation of oxygen tension. Oxygen tension inside the device was computed using commercial finite element software (COMSOL Multiphysics 5.5; COMSOL, Inc, Burlington, MA, USA). The gas flow in the individual channels were simulated by solving the NavierStokes equations coupled with mass continuity for an incompressible fluid:
where u is the velocity vector, p is the pressure, and ${\mathit{\rho}}_{\mathit{G}}$ and ${\mathit{\mu}}_{\mathit{G}}$ are the gas density and viscosity (taken as 1 kg/m^{3} and 10^{−5} Pa.s, respectively). The spatial and temporal distribution of oxygen inside the device was then calculated by solving the convectiondiffusion equation:
where c is the oxygen concentration, D is the diffusion coefficient of oxygen, and T is the time.
The device was assumed to be in an atmosphere containing 21% O_{2}. Medium at 21% O_{2} concentration was supplied to media channels. Gases containing 0% and 21% O_{2} were respectively supplied to the lefthand and righthand side gas channels at 30 ml/min to generate an oxygen gradient. Zero pressure and convection flux conditions were set at the outlets of the gas channels, and a noslip condition was applied on the channel walls for fluid flow analysis. Boundary conditions for oxygen concentration were set according to Henry’s law. Oxygen concentration at the interfaces between the PDMS and gas phase (atmosphere and gas mixture in the gas channel) was set correspondingly to the product of the solubility coefficient of oxygen in PDMS and the partial pressure of oxygen. At the interfaces between PDMS and media or gel, a partition condition was applied, which balanced the mass flux of oxygen to satisfy continuity of partial pressure of oxygen:
where c_{PDMS} and S_{PDMS} are the oxygen concentration and the solubility of oxygen in the PDMS, respectively, and c_{channel} and S_{channel} are those in the media and gel channels. Moreover, oxygen consumption by cells was considered by setting an outward flux of oxygen of 6x10^{−8} [mole/(m^{2}.s)] on the bottom of the media channels (calculated as bρ* where b=1.2.10^{−16} mole/(cell.s) is the oxygen molar consumption per Dd cell per unit of time (Torija et al., 2006) and ρ*=500 cell/mm^{2} is the highest density used in the device).
The diffusion constants of oxygen in the various media were taken to be 2.10^{−9}, 4.1.10^{−9} and 2.10^{−12} m^{2}/s for culture medium, PDMS and PC, respectively. Oxygen solubility at 1atm were taken to be 219 (close to the measured value, see below), 1666 and 1666 μM for culture medium, PDMS and PC, respectively (PDMS and the PC values were assumed to be the same since they are reportedly within the same range Merkel et al., 2000; Moon et al., 2009). The computational models consisted of approximately 1,135,000 computational elements. The initial condition of oxygen concentration in each material was set to 21% O_{2} everywhere (219 μM).
Potts models
Request a detailed protocolPotts model simulations were run using CompuCell3D (Swat et al., 2012) with a mix of prebuilt modules and homemade Python steppables in particular to implement the modulation of aerotactic strength by local oxygen levels. Most parameters were fitted to experimental measurements and both time and length scales were also adapted to achieve quantitative simulations. In all simulations, we used Compucell’s Volume module which applies to all cells a Hamiltonian of the form:
where $V$ is the volume of a cell and ${V}_{cell}$ a target volume set to 2 pixels. This already set the length scale of our simulations to 1pixel = 10µm. ${\lambda}_{v}$ was set to 800. These values were adapted to reproduce the cell speeds observed in the microfluidic experiments. To achieve this relationship, we also decided to fix that one step of the simulation (Monte Carlo Step) was meant to represent 0.1s.
Aerotaxis was modeled using CompuCell’s builtin chemotaxis plugin. This leads to a new term in the Hamiltonian of the form
where $\Delta C$ is the difference in oxygen concentration $C$ between the source and target pixels of a flip and ${\lambda}_{aero}$ is the aerotactic strength. Key to our model is thus the fact that we made ${\lambda}_{aero}$ different for each cell and dependent on the local oxygen concentration. This modulation was fitted to the microfluidic experiments and set, in the general model as:
where C is the oxygen concertation at the center of mass of a cell. Figure 5 shows variations on that relationship which are:
Based on experimental observations, we also made the effective temperature of the model different for each cell and dependent on local oxygen concentrations. This allowed us to reproduce the aerokinetic effect and the modulation of the temperature was fitted to reproduce the cell diffusion constants measured in the microfluidics experiment. The main model thus uses the following relationship for temperature T
Figure 5 and Figure 5—figure supplement 1 show variations on this relationship, which is simply replaced by a constant value of T (115, 135 or 155).
Another key aspect of the models is the oxygen field, which is implemented using CompuCell’s DiffusionSolverFE module. In the case of the microfluidics experiments, the oxygen field was made to be constant (no diffusion, no consumption by cells) and fitted on experimental measurements of this gradient shown in Figure 2B. The oxygen concentrations were expressed in % giving 0 and 21 as natural boundaries. The actual oxygen profile varied in the x direction only as:
where x is the position in pixels in the simulation which was run on a 400 by 800 xy grid.
For spot assays, oxygen was allowed to diffuse freely. In our time and length units, the diffusion constant of oxygen in liquids (2.10^{3} µm^{2}.s^{−1}) is two pixel^{2} step^{−1}.
In terms of consumption, we took the oxygen consumption by Dicty cells to be 1.2.10^{−16} mole/(cell.s) (Torija et al., 2006) and we measured the oxygen solubility in HL5 medium as 250 µM (measurement with a bare fiber sensor plugged to a a Firesting oximeter, Pyroscience, Germany). Taking the measured vertical confinement of 50 µm (Figure 1—figure supplement 1), the amount of oxygen available, at maximum, above a single pixel of the simulation is 1.25.10^{−15} moles, which we define, in arbitrary units, to be 21. We can then turn the consumption of a single cell into a consumption per pixel given that the typical size of a cell is two pixels and per time step, each representing 0.1 s. We end up with a consumption, in our arbitrary units of 0.1 pixel^{−1} step^{−1} which is only applied to pixels occupied by a cell. Of note, in case an occupied pixel had a remaining oxygen level of less than this values, then consumption was set at this oxygen level so that all oxygen was consumed. The last ingredient in oxygen dynamics is the leak of oxygen coming from the bottom of the multiwall plate. Assuming complete hypoxia on the cells’ side, this would lead to a net flux of oxygen of DC/e where D is the oxygen diffusion constant in polystyrene, C is the oxygen concentration on the outside and e is the thickness of the polystyrene bottom. This leads to a flux by unit surface, in our Potts units of 0.001 pixel^{−1}.step^{−1}. We therefore implement a source of oxygen for all pixels in the simulation, whether they are occupied by a cell or not, of the form:
where C is the local oxygen concentration at the considered pixel.
This was sufficient to faithfully reproduce the formation time of the rings. Finally, the spot simulations were run on a 500 by 500 pixels grid and we imposed boundary conditions to the oxygen field as a constant concentration of 21, the borders acting as a source of oxygen just like the edges of the coverslip in the experiments.
Cell division was also set to experimental observations. Given a doubling time of 8 hr, we implemented random divisions at each time point, each cell having a 1/ (8h * 3600 s/h * 10step/s) = 3.10^{−6} chance of dividing. However, cell division was turned off at low oxygen concentrations (<0.7%). In Figure 5, a simulation is shown were this probability was set to 0 for all cells in all conditions.
In terms of initial conditions, the microfluidic simulations were started from a homogenous cell density, each cell being initialized on a grid: two pixels per cells and a six pixel gap to the next neighbor in all directions. For the spot simulations, cells were seeded in three circular, concentric regions of decreasing density. The first region was set to be 30 pixels (300 µm) in radius with a gap of 1 pixel between each cell, the second one spanned the radii between 30 and 60 pixels with a gap of two pixels between each cell and the last one spanned between 60 and 90 pixels with a gap of 3 pixels. This lead to an initial colony with a radius of 900 µm and between 1900 and 2000 cells, both very similar to experiments.
Meanfield model, Go or Grow model and simulations
Request a detailed protocolBoth diffusion equations (Equations 1 and 2) were discretized through a timebackward spacecentered difference scheme with an upwind discretization for the advection operator. In the case of the meanfield model, we were considering (Equations 1 and 2) in a radial symmetry, which lead to the following discretization for $\rho $:
, where ${a}_{i}^{n}=\lambda \left({C}_{i}^{n}\right)\frac{{C}_{i}^{n}{C}_{i1}^{n}}{\Delta x},{r}_{i}^{n}=r\left({C}_{i}^{n}\right)$ and ${R}_{i}$ the distance from the center.
In the case of the Go or grow model with its planar symmetry, the discretization for $\rho $ was:
$\frac{{\rho}_{i}^{n+\frac{1}{2}}{\rho}_{i}^{n\frac{1}{2}}}{\mathrm{\Delta}t}D\frac{{\rho}_{i1}^{n+\frac{1}{2}}2{\rho}_{i}^{n+\frac{1}{2}}+{\rho}_{i+1}^{n+\frac{1}{2}}}{\mathrm{\Delta}{x}^{2}}+\frac{1}{\mathrm{\Delta}x}\cdot \left\{\begin{array}{c}{a}_{i}^{n}{\rho}_{i}^{n+\frac{1}{2}}{a}_{i1}^{n}{\rho}_{i1}^{n+\frac{1}{2}}\text{}if\text{}{C}_{i}^{n}\ge {C}_{i1}^{n}\\ {a}_{i}^{n}{\rho}_{i}^{n+\frac{1}{2}}{a}_{i+1}^{n}{\rho}_{i+1}^{n+\frac{1}{2}}\text{}if\text{}{C}_{i}^{n}{C}_{i1}^{n}\end{array}\right\}={r}_{i}^{n}{\rho}_{i}^{n+\frac{1}{2}}$, with ${a}_{i}^{n}=a\left({C}_{i}^{n}\right)$.
Concerning the equation on Oxygen concentration Equation 2, the consumption term $b\left(C\right)\rho $ was expressed in two different manners: either $b\left(C\right)={b}_{0}$ and a nonnegativity constraint was added on the Oxygen concentration, just as it was the case in the cellular Potts model, or $b\left(C\right)=min\left({b}_{0},{b}_{0}\frac{C}{{C}_{0}\text{'}}\right)$ which leads to an oxygen consumption that goes to zero in the region of very low concentration $C<{C}_{0}^{\prime}$ and therefore ensures nonnegativity for $C$ under a sufficiently small time step $\Delta t$. Both expressions led to qualitatively similar results, but we opted for the latter in all the simulations presented here. Finally, the discretization scheme for $C$ in the planar symmetry was:
$\frac{{C}_{i}^{n+1}{C}_{i}^{n}}{\Delta t}{D}_{oxy}\frac{{C}_{i1}^{n+1}2{C}_{i}^{n+1}+{C}_{i+1}^{n+1}}{{\Delta x}^{2}}={b}_{i}^{n}{\rho}_{i}^{n+\frac{1}{2}}$, with ${b}_{i}^{n}=b\left({C}_{i}^{n}\right)$. For the radial symmetry:
The schemes were coded in Python language. All simulations of (Equations 1 and 2) shown in this article were carried out with a mesh size $\Delta t=0.02min$ and $\Delta x=1\mu m$. The values used for the constants are: $D=30{\mu m}^{2}\cdot {min}^{1}$ (effective cellular diffusion constant), ${C}_{0}=0.7\text{\%}{O}_{2}$ (threshold for cell division), ${C}_{0}\text{'}=0.1\text{\%}{O}_{2}$ (lower threshold in the two threshold ‘Go or Grow’ model, below which cells stop aerotaxis), ${D}_{oxy}=1.2\cdot {10}^{5}{\mu m}^{2}\cdot {min}^{1}$ (oxygen diffusion constant in medium), ${r}_{0}=ln2/480{min}^{1}$ (rate of cell division) and ${b}_{0}=0.01\text{\%}{O}_{2}{min}^{1}{cell}^{1}$ (using the equivalence 1.25.10^{−15} moles = 21 %O_2 discussed above in Potts model section).
The scheme on $C$ was supplemented with the boundary condition $C\left(L\right)=21\text{\%}{O}_{2}$. We have chosen $L=9mm$ for the Go or grow model to match experimental conditions. For the meanfield model, we have chosen $L=2.5mm$ in order to match the cellular Potts model for which size was a concern for computation time.
In the meanfield model, the initial condition for $\rho $ was taken the same as in the cellular Potts model. For the other simulations, initial conditions for $\rho $ and $C$ were chosen such that they were already close to the expected stationary profile.
We measured the speed of the wave $\sigma $, once the wave profile was qualitatively stable, by considering the evolution of the point $\overline{x}\left(t\right)$ such that $C\left(t,\overline{x}\left(t\right)\right)={C}_{0}$.
Mathematical analysis of the ‘Go or Grow’ model
Request a detailed protocolWe present below a preliminary analysis of the ‘Go or Grow’ model. A more detailed mathematical investigation of this model will be carried out in a separate article.
1.The ‘Go or Grow’ model admits explicit traveling wave solutions.
We recall that $z=x\sigma t$ is the spatial variable in the moving frame at (unknown) speed $\sigma >0$. We seek a pair of stationary profiles, resp. the density $\rho \left(z\right)$ and the oxygen concentration $C\left(z\right)$. We assume that $C\left(z\right)$ is an increasing function. By translation invariance, we set without loss of generality that $C\left(0\right)={C}_{0}$, so that Equation 3 becomes:
Furthermore, the function $\rho \left(z\right)$ must satisfy at $z=0$ the following relation (i.e. the continuity of the flux) :
Thus the equation becomes a second order differential equation with piecewise constant coefficients on each halfline, that can be solved explicitly.
For $z<0$, the solution is of the form $A+B{e}^{\frac{{a}_{0}\sigma}{D}z}$. From Equation 4 we observe that $\sigma \ge {a}_{0}$ equivalently $\frac{{a}_{0}\sigma}{D}\le 0$ and as $\rho $ is bounded, it implies that $B=0$.
For $z>0$, we look at the roots of the characteristic polynomial $P\left(\mu \right)=D{\mu}^{2}+\sigma \mu +{r}_{0}$. We note that to yield a nonnegative solution, we need ${\sigma}^{2}\ge 4{r}_{0}D$.
If $\sigma =2\sqrt{{r}_{0}D}$, then the solution is of the form $\left(Cz+D\right){e}^{\sqrt{\frac{{r}_{0}}{D}}z}$ and with relation Equation 8, we obtain $\rho \left(z\right)=A\left(\frac{\sqrt{{Dr}_{0}}{a}_{0}}{D}z+1\right){e}^{\sqrt{\frac{{r}_{0}}{D}}z}$ and observe that in this case, we necessarily have ${a}_{0}\le \sqrt{{r}_{0}D}$.
If $\sigma >2\sqrt{{r}_{0}D}$, the solution is of the form $A\text{'}{e}^{{\mu}_{}z}+B\text{'}{e}^{{\mu}_{+}z}$, with ${\mu}_{\pm}=\frac{\sigma \pm \sqrt{{\sigma}^{2}4{r}_{0}D}}{2D}$. By arguments exposed in Van Saarloos, 2003, solutions with initial datum localized cannot decrease exponentially at a rate $\mu >\sqrt{\frac{{r}_{0}}{D}}$, where $\sqrt{\frac{{r}_{0}}{D}}$ corresponds to the exponential decay parameter when $\sigma =2\sqrt{{r}_{0}D}$. This leads to $B\text{'}=0$, as $\mu}_{+}>\sqrt{\frac{{r}_{0}}{D}$. Then $A=A\text{'}$, but in order to satisfy the C^{1}discontinuity jump relation Equation 7, it must be that :
Equation 9 can be solved algebraically for $\sigma $, which yields $\sigma ={a}_{0}+\frac{{Dr}_{0}}{{a}_{0}}$. Furthermore, we can rewrite Equation 9 as follows $\left(2{a}_{0}\sigma \right)=\sqrt{{\sigma}^{2}4{r}_{0}D}$, multiplying by $2{a}_{0}+\sigma $, we find that $\left(4{a}_{0}^{2}{\sigma}^{2}\right)=\left(2{a}_{0}+\sigma \right)\sqrt{{\sigma}^{2}4{r}_{0}D}>0$, which leads to ${a}_{0}^{2}>\frac{{\sigma}^{2}}{4}>{r}_{0}D$.
Thus, we have disclosed all possible profiles. In the case ${a}_{0}\le \sqrt{{r}_{0}D}$ the profile travels at speed $\sigma =2\sqrt{{r}_{0}D}$, whilst for $a}_{0}>\sqrt{{r}_{0}D$ the profile travels at speed $\sigma ={a}_{0}+\frac{{r}_{0}D}{{a}_{0}}$.
One needs to verify that each of these profiles admits an associated oxygen profile that satisfies the condition $C\left(0\right)={C}_{0}$, but the preceding profiles were defined up to the multiplicative constant $A$, by linearity of Equation 7. The differential equation on $C$ becomes with $\stackrel{~}{\rho}$ the solution given above for $A=1$:
One concludes by checking that by monotonicity there exists a unique constant $A$ such that the solution to the differential Equation 10 equation satisfies $C\left(0\right)={C}_{0}$.
2. The wave is pushed in the case $a}_{0}>\sqrt{{r}_{0}D$.
A neutral fraction ${v}^{k}$ is defined as satisfying the following linear equation in the moving frame $z=x\sigma t$:
with ${v}^{k}\left(0,z\right)={v}_{0}^{k}\left(z\right)$ where we identify $a\left(z\right)=a\left(C\left(z\right)\right)$ and $r\left(z\right)=r\left(C\left(z\right)\right)$ for the sake of clarity. This corresponds biologically to staining the cells given by the initial distribution ${v}_{0}^{k}$ at time $t=0$ with a neutral label (Roques et al., 2012).
Defining $U\left(z\right):=\frac{\left(\sigma a\left(z\right)\right)}{D}z$, then we note that $Lf=D\frac{\partial}{\partial z}\left({e}^{U}\frac{\partial}{\partial z}\left({e}^{U}f\right)\right)r\left(z\right)f$. This leads to setting $w:={e}^{\frac{U}{2}}{v}^{k}$ that satisfies the parabolic equation $\frac{\partial w}{\partial t}+\stackrel{~}{L}w=0$, with $\stackrel{~}{L}g:={De}^{\frac{U}{2}}\frac{\partial}{\partial z}\left({e}^{U}\frac{\partial}{\partial z}\left({e}^{\frac{U}{2}}g\right)\right)r\left(z\right)g=D\frac{{\partial}^{2}g}{\partial {z}^{2}}+\left(\frac{U{\text{'}}^{2}}{4}r\left(z\right)\frac{{a}_{0}}{2}\delta \right)g$. The operator $\stackrel{~}{L}$ is selfadjoint in ${L}^{2}\left(\mathbb{R},dz\right)$ on the appropriate domain. Then by setting $\gamma :=min\left(inf\left\{\frac{{U}^{\prime 2}}{4}r\left(z\right)\right\},{\left(\frac{{r}_{0}D}{{a}_{0}}\right)}^{2}\right)>0$, one can first show that every element of the spectrum of $\lambda \in \sigma \left(\stackrel{~}{L}\right)$ such that $\lambda <\gamma$ is an eigenvalue of $\stackrel{~}{L}$. Second, one shows that the only eigenvalue $\lambda $ of $\stackrel{~}{L}$ such that $\lambda <\gamma$ is $\lambda =0$. Finally by standard theory of selfadjoint operators and semigroup theory, one obtains that $w\left(t\right)={Pw}_{0}+{e}^{t\stackrel{~}{L}}\left(IP\right){w}_{0}$, where ${\Vert {e}^{t\stackrel{~}{L}}\left(IP\right){w}_{0}\Vert}_{{L}^{2}\left(\mathbb{R},dz\right)}\u2a7d{e}^{\gamma t}{\Vert {w}_{0}\Vert}_{{L}^{2}\left(\mathbb{R},dz\right)}$. Translating these properties onto the neutral fraction ${v}^{k}$, we have that ${v}^{k}\left(t\right)\to \frac{{\u27e8{v}_{0}^{k},\rho \u27e9}_{{L}^{2}\left(\mathbb{R},{e}^{U}dz\right)}}{{\u27e8\rho ,\rho \u27e9}_{{L}^{2}\left(\mathbb{R},{e}^{U}dz\right)}}\rho $ at an exponential rate, where $\rho $ is the traveling wave profile calculated in the previous section. Therefore, each fraction converges to a fixed proportion of the whole population. We conclude that after some time the wave becomes a perfect mix of each neutral fraction. This corresponds to the definition of a pushed wave according to Roques et al., 2012.
3. The wave is pulled in the case ${a}_{0}\le \sqrt{{r}_{0}D}$.
The preceding reasoning does not apply to this case and the intuition is clear, as the wave speed coincides with Fisher’s $\sigma =2\sqrt{{r}_{0}D}$, which is typically the signature of a pulled reaction diffusion front. In order to prove the pulled nature of the front, we consider ${w}^{k}=\frac{{v}^{k}}{\rho}$, where $\rho $ is the corresponding wave profile. By computation, ${w}^{k}$ then satisfies the following PDE:
$\frac{\partial {w}^{k}}{\partial t}\beta \left(z\right)\frac{\partial {w}^{k}}{\partial z}D\frac{{\partial}^{2}{w}^{k}}{\partial {z}^{2}}=0$ with $\beta \left(z\right)=\{\begin{array}{l}\frac{2\sqrt{{r}_{0}D}{a}_{0}}{D},ifz<0\\ 2\frac{\sqrt{{r}_{0}D}{a}_{0}}{\left(\sqrt{{r}_{0}D}{a}_{0}\right)z+1},ifz\ge 0\end{array}$ and set $\eta $ a positive solution to the differential equation $\eta \text{'}=\beta \eta $. As $\beta \left(z\right)\ge 0$ and $\beta \text{'}$ bounded above, it can be shown by arguments similar to Roques et al., 2012, that under the integrability condition $\int {\left({w}_{0}^{k}\left(z\right)\right)}^{2}\eta \left(z\right)dz<\mathrm{\infty}$, the neutral fraction goes extinct, that is $\underset{t=+\infty}{\mathit{lim}}\Vert {\left({w}^{k}\right)}^{2}\eta {\Vert}_{\infty}=0$, which characterizes a pulled wave in the framework of neutral fractions.
Mathematical analysis of a specific ‘Go or Grow’ model with a second threshold
Request a detailed protocolWe present quickly a specific case for a ‘Go or Grow’ model with a second threshold, that is completely analytically solvable. We consider the advection term of the form $a\left(C\right)sign\left({\partial}_{x}C\right)$ with $a\left(C\right)=\{\begin{array}{l}{a}_{0},\text{}\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}if\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\text{}{C}_{0}^{\prime}C{C}_{0}\\ 0,\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}otherwise\end{array}$, the division rate $r\left(C\right)=\{\begin{array}{l}{r}_{0},\text{}\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}if\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\text{}C{C}_{0}\\ 0,\text{}\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}if\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\text{}C{C}_{0}\end{array}$ and the O_{2} consumption rate per cell $b\left(C\right)={b}_{0}$, without including the constraint that the O_{2} concentration $C$ be nonnegative. Although this hypothesis seems physically non relevant, it is consistent with the fact that cells are not sensitive to O_{2} concentration gradients below the threshold ${C}_{0}\text{'}$.
Given a traveling wave profile $\rho ,C$ and the corresponding front speed $\sigma $, we suppose $C\left(0\right)={C}_{0}$ and we introduce the spatial gap $h>0$ between the two thresholds, i.e. $C\left(h\right)={C}_{0}\text{'}$, so that (Equation 1) becomes:
Introducing a multiplicative constant $A$, $\rho $ is then of the shape:
With $B=\frac{\sigma \mu D}{\sigma {a}_{0}}$, $E=\frac{\mu D{a}_{0}}{\sigma {a}_{0}}$ and $\mu =\frac{\sigma +\sqrt{{\sigma}^{2}4D{r}_{0}}}{2D}$. We obtain the following condition, that establishes a onetoone correspondence between $\sigma $ and $h$:
The equation on $C$ becomes:
With the assumption that that $C$ be continuously differentiable, we can solve Equation 15 for $C$:
With $F=\frac{{b}_{0}}{\sigma}\left(B+{Ee}^{\frac{\sigma {a}_{0}}{D}h}\right)$, $G={C}_{0}\text{'}+FAh$, $H=\frac{{b}_{0}B}{\sigma}$, $I=\frac{{b}_{0}{D}^{2}E}{\left(\sigma {a}_{0}\right)\left(D\sigma {D}_{oxy}\left(\sigma {a}_{0}\right)\right)}$, $L=\frac{{b}_{0}}{\mu \left({\mu D}_{oxy}\sigma \right)}$, $M={C}_{0}LA{C}_{init}$, $J={C}_{0}IAK$ and, by setting $\Delta =\frac{\sigma}{D}\left(HFI\left(\frac{\sigma {a}_{0}}{D}\right){e}^{\frac{\sigma {a}_{0}}{D}h}\right)\left(H+L\left(\mu \frac{\sigma}{{D}_{oxy}}\right)I\left(\frac{\sigma {a}_{0}}{D}\right)\right)\left(\frac{\sigma}{{D}_{oxy}}{e}^{\frac{\sigma}{{D}_{oxy}}h}\right)$, we have that $K=\frac{1}{\Delta}\left(HFI\left(\frac{\sigma {a}_{0}}{D}\right){e}^{\frac{\sigma {a}_{0}}{D}h}\right)\frac{\sigma}{{D}_{oxy}}\left({C}_{init}{C}_{0}\text{'}\right)$ and $A=\frac{1}{\Delta}\left(\frac{\sigma}{{D}_{oxy}}{e}^{\frac{\sigma}{{D}_{oxy}}h}\right)\frac{\sigma}{{D}_{oxy}}\left({C}_{init}{C}_{0}\text{'}\right)$. This closes the system, but one more constraint remains, which is:
The front speed $\sigma $ of a traveling wave must therefore satisfy the implicit Equation 16. Finding a closed form for the solutions of Equation 16 seems out of reach. Nevertheless, we can approximate the roots numerically, especially by noticing through numerical observation that Equation 16 is monotone on the interval $\left(2\sqrt{{r}_{0}D},{a}_{0}+\frac{{r}_{0}D}{{a}_{0}}\right)$, where the root $\sigma $ is located. Hence through a dichotomy search algorithm we can find the speed $\sigma $ of the traveling wave with arbitrary accuracy.
Data availability
All data generated or analysed during this study are included in the manuscript and supporting files.
References

BookNonlinear Diffusion in Population Genetics, Combustion, and Nerve Pulse PropagationBerlin, Heidelberg: Springer.https://doi.org/10.1007/BFb0070595

Mitochondria, redox signaling and axis specification in metazoan embryosDevelopmental Biology 308:266–280.https://doi.org/10.1016/j.ydbio.2007.05.042

Properties of germinating spores of Dictyostelium discoideumJournal of bacteriology 96:1680–1689.https://doi.org/10.1128/jb.96.5.16801689.1968

Changes in the magnitude and distribution of forces at different Dictyostelium developmental stagesCell Motility and the Cytoskeleton 65:314–331.https://doi.org/10.1002/cm.20262

Collective regulation of cell motility using an accurate densitysensing systemJournal of the Royal Society Interface 15:45.https://doi.org/10.1098/rsif.2018.0006

Neue methode zur untersuchung der sauerstoffausscheidung pflanzlicher und tierischer organismen (New method for investigation of oxygensearching plant and animal organisms)Pflugers Arch. Gesammte Physiol 25:285–292.

The wave of advance of advantageous genesAnnals of Eugenics 7:355–369.https://doi.org/10.1111/j.14691809.1937.tb02153.x

Amoeboid leukocyte crawling through extracellular matrix: lessons from the Dictyostelium paradigm of cell movementJournal of Leukocyte Biology 70:491–509.

'Go or grow': the key to the emergence of invasion in tumour progression?Mathematical Medicine and Biology 29:49–65.https://doi.org/10.1093/imammb/dqq011

A user’s guide to PDE models for chemotaxisJournal of Mathematical Biology 58:183–217.https://doi.org/10.1007/s0028500802013

Initiation of slime mold aggregation viewed as an instabilityJournal of Theoretical Biology 26:399–415.https://doi.org/10.1016/00225193(70)900925

Etude de L’équation De La Diffusion Avec Croissance De La Quantité De Matière Et Son Application À Un Problème BiologiqueBull. Univ. Moskow,Ser. Internat 1:1–25.

Gas sorption, diffusion, and permeation in poly(dimethylsiloxane)Journal of Polymer Science Part B: Polymer Physics 38:415–434.https://doi.org/10.1002/(SICI)10990488(20000201)38:3<415::AIDPOLB8>3.0.CO;2Z

Outgassing of oxygen from polycarbonateACS Applied Materials & Interfaces 1:1539–1543.https://doi.org/10.1021/am900206e

Rectified directional sensing in longrange cell migrationNature Communications 5:1–14.https://doi.org/10.1038/ncomms6367

New horizons in hypoxia signaling pathwaysExperimental Cell Research 356:116–121.https://doi.org/10.1016/j.yexcr.2017.03.008

Expression of cytochrome c oxidase during growth and development of DictyosteliumJournal of Biological Chemistry 270:5587–5593.https://doi.org/10.1074/jbc.270.10.5587

Oxygen influences the subunit structure of cytochrome c oxidase in the slime mold Dictyostelium discoideumJournal of Biological Chemistry 264:7129–7134.https://doi.org/10.1016/S00219258(18)832112

On two types of moving front in quasilinear diffusionMathematical Biosciences 31:307–315.https://doi.org/10.1016/00255564(76)900870

Moving in and out: dispersion of cells in SelfGenerated gradientsJournal of Clinical & Cellular Immunology 8:507.https://doi.org/10.4172/21559899.1000507

Multiscale modeling of tissues using CompuCell3DMethods in Cell Biology 110:325–366.https://doi.org/10.1016/B9780123884039.000138

A noninvasive thin film sensor for monitoring oxygen tension during in vitro cell cultureAnalytical Chemistry 81:9239–9246.https://doi.org/10.1021/ac9013379

In vitro metabolic zonation through oxygen gradient on a chipScientific Reports 9:13557.https://doi.org/10.1038/s41598019494126

Selfgenerated chemotactic gradientscells steering themselvesCurrent Opinion in Cell Biology 42:46–51.https://doi.org/10.1016/j.ceb.2016.04.003

SelfGenerated gradients yield exceptionally robust steering cuesFrontiers in Cell and Developmental Biology 8:133.https://doi.org/10.3389/fcell.2020.00133

Coupled excitable Ras and Factin activation mediates spontaneous pseudopod formation and directed cell movementMolecular Biology of the Cell 28:922–934.https://doi.org/10.1091/mbc.e16100733

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

Tatjana PiotrowskiReviewing Editor; Stowers Institute for Medical Research, United States

Naama BarkaiSenior Editor; Weizmann Institute of Science, Israel

Roeland MH MerksReviewer; Leiden University, Netherlands

Robert InsallReviewer
In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.
Acceptance summary:
The manuscript describes a very interesting novel observation of the collective response of cells of cellular slime mold Dictyostelium discoideum: it forms dense rings in response to oxygen gradients. Using carefully executed experimental work combined with thorough numerical and mathematical modeling the authors plausibly show that the phenomenon can be explained through a combination of concentrationdependent taxis to oxygen gradients, and increased motility at low oxygen concentrations. The work is of interest to developmental biologists and mathematical biologists, as well as to cancer researchers due the similarity to cancer cell migration under hypoxia.
Decision letter after peer review:
Thank you for submitting your article "Hypoxia triggers collective aerotactic migration in Dictyostelium discoideum" for consideration by eLife. Your article has been reviewed by 2 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Naama Barkai as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Roeland MH Merks (Reviewer #1); Robert Insall (Reviewer #2).
The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.
Summary:
The manuscript describes a highly interesting novel observation of the collective response of cells of cellular slime mold Dictyostelium discoideum: it forms dense rings in response to oxygen gradients. Using experimental work combined with numerical and mathematical modeling the authors plausibly show that the phenomenon can be explained through a combination of concentrationdependent taxis to oxygen gradients, and increased motility at low oxygen concentrations. The findings are relevant to a surprisingly wide range of biomedicine – especially tumour growth, but also carnivore foraging (particularly amoebas). However, the biological relevance of the phenomenon is not clear and the model predictions are not evaluated experimentally. The major conclusions are largely justified, but there is little insight in the correctness or explanatory power of the model predictions. With one or a few simple validation experiments, the model predictions in Figure 5 could probably be tested. This, or even an indepth discussion of the implications of the model predicts would strengthen the support for the major conclusions. In particular this would complete the paper, which currently seems to divert a little into every more detailed and thorough modeling.
Essential revisions:
Strong points of the work are the interesting novel phenomenon of aerotaxis and aerokinesis in Dd, the thorough quantitative data, including quantitative characterization of the oxygen gradients and cell velocities, as well as the detailed mathematical modeling. However, there are some weaknesses that need to be addressed:
1) The biological relevance of the aerotactic behavior is not discussed. Is ring formation observed under more 'natural' laboratory conditions, or is ring formation a consequence of the aerotactic and aerokinetic behavior that is relevant under other biological conditions? At what life stage of Dd is the behavior important: the single cell or collective stage? Is ring formation true 'collective behavior' like fruiting body formation or e.g. thermo or phototactic behavior of slugs, or should one see it as a side effect of many individual cells performing the same singlecell behavior together under unusual laboratory conditions (as the only interaction between cells for this phenomenon is through the oxygen levels in the microenvironment). Or can it be considered a model system, e.g., for cancer cell migration under hypoxic conditions.
2) The authors first present a cellbased model, and then they study a "mean field" model that is based on the same underlying biological assumptions as the cellbased model. Interestingly, however, the cellbased model explains the data well, whereas for the meanfield model additional assumptions (a second threshold below which no aerotaxis occurs or alternative response functions) are required. What explains this discrepancy between the two models? We would expect that the meanfield model has the same overall behavior as the cellbased model upon which it is based, unless the assumptions differ or the discrete interactions between cells introduce new dynamics. For example, one could imagine that cell adhesion represented in the CPM strengthens the formation of the ring (i.e. reduces cell diffusion at higher cell densities). Could this or other "hidden" features of the CPM be introduced into the PDE and suffice to explain the discrepancy? What is the effect of the CPM parameters on ring formation?
3) The cellbased models and mathematical model produce interesting predictions as outlined in Figure 5 and in the expressions of wave speed for example. Yet the modeling is an end point. The models are compared with the goorgrow models and with the FisherKPP reactiondiffusionmodel, but the predictions are not used to further analyse the system, which makes the story not 'round' – the modeling seems to get ever more detailed until it reaches a sort of 'dead end'. Yet many of the predictions can be tested. E.g. divisions (Figure 5B) can be inhibited. The oxygen gradients can be manipulated. Motility (D) can be changed etc.
4) The authors are making conclusions about "aerokinesis" because cells experiencing a gradient but behind the crest are highly motile. Another explanation is that the motile cells' behaviour is due to different (greater?) resolvable directional information (Tweedy and Endres, Scientific Reports, around 2014). Could the authors perform experiments to test whether cells at the "aerokinetic" O2 levels behave the same way if they are not experiencing a gradient? Can the authors keep cells in a homogeneous O2 at that level and see if they still polarise and move fast? If so it is indeed aerokinesis. Alternatively, the authors could discuss the issue.
5) Around line 497 the authors make some strong conclusions about the generality of the system. But it seems they are all based on a strictly 2D system, and some of the most interesting conditions are nearly 1D (cells moving down a pipe) or 3D (real life colonies of amoebas or bacteria in soil; there's lots of 2D action in soil but O2 is likely freely available in most). Could the authors (a) try a 1D experiment (b) consider a 3D thought experiment or model? In any case they should discuss how the 1D and 3D cases differ and how important they might be. Maybe even semi3D (thickish layers?).
6) An open request, which would help future work to pick the molecular bones of this work – is for an overt description of the concentrations of O2 that provoke the strongest responses.
Reviewer #1:
The manuscript by CochetEscartin et al. describes a novel, aerotactic behavior of Dictyostelium discoideum (Dd) and analyses in it in great detail using experimental and mathematical tools, leading to compelling quantitative data of the phenomenon. When Dd colonies are covered, they migrate outwards towards highoxygen conditions in a dense ring. Inside the central zone surrounded by the ring, Dd cells assume a quiescent, elongated phenotype, while within the ring and outside cells are more rounded and migratory. The behavior of the cells was hypothesized to be regulated by oxygen, and therefore the response of the Dd cells was studied inside a microfluidic device allowing tight control of the oxygen gradient. In particular it was found that a low oxygen levels (<2 %) the cells became sensitive to oxygen gradients, where they moved to higher levels of oxygen. After that using an oxygen sensitive substrate, it was found that oxygen gradients are selfgenerated by Dd cells due to oxygen consumption, leading to the formation of the ring.
Based on these observations, the authors proposed two types of mathematical models. The first model is a cellbased simulation model based on a hybrid Cellular Potts and PDE model. It assumes a field of oxygen that is formed through leakage of oxygen through the polystyrene bottom of and the cells' consumption of oxygen. The cells switch from a quiescent to a motile state if the oxygen level drops below a threshold, at the motile state they become aerotaxic and they proliferate. These assumptions suffice to reproduce the observed ring formation, oxygen gradient and the low density of motile cells in the center of the colony.
The second model is a continuum model based on the same assumptions, leading to a 'go or grow' type reactionadvectiondiffusion equation. Aerotaxis again occurs only at low oxygen concentration (i.e. the advection term) and there is an oxygendependent proliferation term that lets cells proliferate only at abovethreshold oxygen levels. Interestingly the model can be solved in part analytically, leading to an expression for the wave speed that depends on the growth rate as well as the diffusion rate and for sufficiently strong aerotaxis also on the aerotaxis. The numerical simulations do not match the observations well: in particular the concentration of cells in the core of the ring remains much higher than what is observed in the data (e.g. Figure 6C). In order to "fix" this discrepancy between model and data, a further assumption is introduced, namely a second, lower threshold C_{0}' below which the cells become insensitive to the gradients. I.e. only in the range [C_{0}' – C_{0}] do the cells become aerotactic. Also other assumptions lead to a better match with the data.
Strong points of the work are the interesting novel phenomenon of aerotaxis and aerokinesis in Dd, the thorough quantitative data, including quantitative characterization of the oxygen gradients and cell velocities, as well as the detailed mathematical modeling.
However, the biological relevance of the aerotactic behavior is not discussed at all. It is not clear whether ring formation is observed under more 'natural' laboratory conditions, or whether ring formation is a consequence of the aerotactic and aerokinetic behavior that is relevant under other biological conditions. It is also not clear at what life stage of Dd this behavior is important, the single cell or collective stage, and whether ring formation is true 'collective behavior' like fruiting body formation or e.g. thermo or phototactic behavior of slugs. Or should we see it as a side effect of many individual cells performing the same singlecell behavior together under unusual laboratory conditions? I would tend towards the latter (lab conditions) as the only interaction between cells for this phenomenon is through the oxygen levels in the microenvironment. Or should we really see this as a model system, e.g., for cancer cell migration under hypoxic conditions?
Furthermore, the authors first present a cellbased model, and then they study a "mean field" model that is based on the same underlying biological assumptions as the cellbased model. Interestingly, however, the cellbased model explains the data well, whereas for the meanfield model additional assumptions (a second threshold below which no aerotaxis occurs or alternative response functions) are required. What explains the discrepancy between the two models is not clarified. We would expect that the meanfield model has the same overall behavior as the cellbased model upon which it is based, unless the assumptions differ or the discrete interactions between cells introduce new dynamics.
The cellbased models and mathematical model produce interesting predictions as outlined in Figure 5 and in the expressions of wave speed for example. Yet the modeling is an end point. The models are compared with the goorgrow models and with the FisherKPP reactiondiffusionmodel, but the predictions are not used to further analyse the system, which makes the story not 'round' – the modeling seems to get ever more detailed until it reaches a sort of 'dead end'. Yet many of the predictions can be tested. For example, divisions (Figure 5B) can be inhibited, the oxygen gradients can be manipulated and motility (D) can be changed.Reviewer #2:
The authors combine wet experiments, cellular automaton modelling and analytical maths in a dazzling way to show that cells collectively run away from hypoxia (that they have collectively generated). This is a really cool paper. It's relevant to a suprisingly wide range of biomedicine – especially tumour growth, but also carnivore foraging (particularly, of course, amoebas). It's in an exciting and visible area of science. In a sense it's obvious, but – after reading this paper. The data and models are mutually supportive and combine to make a thorough and scholarly exploration of the area.The combination of wet experiments, models and math is this paper's strength; bringing different informative figures together to make a clear conclusion. Its weakness to some will be a lack of description of the underlying pathways.
Two areas where there appears to be lack of clarity:
I'm not sure it's wise to conclude aerokinesis (that is, particular O2 concentrations without a gradient promote nondirectional movement) without an explicit test. If I've understood right the cells have not been seen moving faster at low O2 with no gradient; they've just been seen going particularly rapidly at one point in the gradient.
Another explanation (Tweedy and Endres, Scientific Reports, around 2014; the same author as the later selfgenerated gradient papers) is that the motile cells' behaviour is due to different (greater?) resolvable directional information.
Around l. 497 the authors make some strong conclusions about the generality of the system. But it seems that they are all based on a strictly 2D system, and some of the most interesting conditions are nearly 1D (cells moving down a pipe) or 3D (real life colonies of amoebas or bacteria in soil; there's lots of 2D action in soil but I think O2 is freely available in most).
https://doi.org/10.7554/eLife.64731.sa1Author response
Summary:
The manuscript describes a highly interesting novel observation of the collective response of cells of cellular slime mold Dictyostelium discoideum: it forms dense rings in response to oxygen gradients. Using experimental work combined with numerical and mathematical modeling the authors plausibly show that the phenomenon can be explained through a combination of concentrationdependent taxis to oxygen gradients, and increased motility at low oxygen concentrations. The findings are relevant to a surprisingly wide range of biomedicine – especially tumour growth, but also carnivore foraging (particularly amoebas). However, the biological relevance of the phenomenon is not clear and the model predictions are not evaluated experimentally. The major conclusions are largely justified, but there is little insight in the correctness or explanatory power of the model predictions. With one or a few simple validation experiments, the model predictions in Figure 5 could probably be tested. This, or even an indepth discussion of the implications of the model predicts would strengthen the support for the major conclusions. In particular this would complete the paper, which currently seems to divert a little into every more detailed and thorough modeling.
We thank the Reviewers and Editors for their comments. In response, many aspects of the papers have been modified to address the concerns raised. We agree that our findings are “relevant to a wide range of biomedicine” questions but that our modelling work was somewhat too focused on our own experiments and, in places, failed to demonstrate their generality. As a result, particular attention was given to this part of the paper which has been almost entirely reworked. We now show that our analytical approach can not only reproduce our own experimental data but give insight into front propagation under a single selfgenerated gradient regardless of the underlying details of the process. Along the same line, we have also added new experiments and simulations to further test the predictions of our models. The introduction and Discussion section have been significantly reworked to include discussions on the generality and biological relevance of our results. Finally, we now demonstrate the aerokinesis of Dd cells in homogenous environments through new experiments. Below, we offer a point by point response to the essential revisions which recapitulated all concerns raised by the reviewers.
Essential revisions:
Strong points of the work are the interesting novel phenomenon of aerotaxis and aerokinesis in Dd, the thorough quantitative data, including quantitative characterization of the oxygen gradients and cell velocities, as well as the detailed mathematical modeling. However, there are some weaknesses that need to be addressed:
1) The biological relevance of the aerotactic behavior is not discussed. Is ring formation observed under more 'natural' laboratory conditions, or is ring formation a consequence of the aerotactic and aerokinetic behavior that is relevant under other biological conditions? At what life stage of Dd is the behavior important: the single cell or collective stage? Is ring formation true 'collective behavior' like fruiting body formation or e.g. thermo or phototactic behavior of slugs, or should one see it as a side effect of many individual cells performing the same singlecell behavior together under unusual laboratory conditions (as the only interaction between cells for this phenomenon is through the oxygen levels in the microenvironment). Or can it be considered a model system, e.g., for cancer cell migration under hypoxic conditions.
We thank the Reviewer for these comments and apologize for not being clearer. We agree that our results are based on an artificial situation under controlled laboratory conditions and do not directly describe a real life situation. However, we demonstrate throughout the paper that this behavior is governed by aerotaxis and selfgenerated gradients. Our paper demonstrates this ability in a new system, Dd cells. By demonstrating the ability of Dd cells to perform both aerokinesis and aerotaxis, we agree that “ring formation is a consequence of aerotactic and aerokinetic behavior that is relevant under other biological conditions.”
Indeed, aerotaxis is a strategy to avoid unfavorable O_{2} conditions probably conserved among all eukaryotes. It may play a role in placentation, in tumor progression by orienting the migration and invasion of the cells of the primary hypoxic tumor to the blood capillaries, promoting metastatic spread. We believed these examples discussed in the second paragraph of the manuscript introduction demonstrate that aerotaxis is relevant in various physiological situations. We added another interesting notion that gradients of O_{2} and energy metabolism govern spatial patterning in various embryos. This notion dates back to the classic work of Child (Child 1941). Such notions have mostly been abandoned due to inability to visualize such a gradient or clarify whether they are the result or the cause of developmental patterning (Coffman and Denegre 2007). More recently, Chang et al. found an asymmetric distribution of hypoxiainducible factor regulating dorsoventral axis establishment in the early sea urchin embryo (Chang et al. 2017). Interestingly, they also found evidence for an intrinsic hypoxia gradient in embryos, which may be a forerunner to dorsoventral patterning. They did not measure directly oxygen nor explained the possible origin of this gradient. We believe new tools, experimental and theoretical, such as those developed in this work, are now available to tackle this question. We hence added the following text to the introduction of the paper:
Addition to the text: “Aerotaxis may also play a role in morphogenesis. […] Interestingly, they also found evidence for an intrinsic hypoxia gradient in embryos, which may be a forerunner to dorsoventral patterning.”
Regarding selfgenerated gradients, besides the case of O_{2}, such gradients of different molecules consumed or degraded by cells are increasingly encountered in the literature. We already cited the important contribution of Reviewer #2 in the introduction (paragraph 3) and in the conclusion (paragraph 4). We again insist on the very diverse patterns including rings that have been found for a long time with bacteria when food or the diffusion of some metabolite including O_{2} is limited (J Adler 1966)(Budrene and Berg 1991). Another interesting system is in plants. Negative plant soil feedback also explains ring formation in clonal plants (Carteni et al., 2012). One key aspect of selfgenerated gradients that is currently emerging is their ability to drive longlasting migrations in complex environments and to enhance the robustness of collective migrations. These physiologically relevant aspects of selfgenerated gradients were demonstrated by the team of Reviewer #2 (Tweedy and Insall, 2020, Tweedy et al., 2020) and others (Cremer et al., 2019). We have expanded on this aspect in the introduction of the paper to underline the physiological relevance of collective migrations driven by selfgenerated gradients.
Addition to the text: “Physiologically speaking, selfgenerated gradients have been demonstrated to increase the range of expansion of cell colonies (Cremer et al. 2019; Tweedy and Insall 2020) and to serve as directional cues to help various cell types navigate complex environments, including mazes (Tweedy et al. 2020). Oxygen selfgenerated gradients could therefore play important roles in a variety of contexts, such as development, cancer progression or even environmental navigation in the soil.”
We have also added a few sentences at the beginning of the Discussion section to more directly highlight the possible physiological relevance of our observations.
Addition to the text: “Overall, our results demonstrate the ability of Dd cells to respond to hypoxia through both aerotactic and aerokinetic responses. […] Finally, our work goes beyond these results as it demonstrates that oxygen can play the role of the attractant in selfgenerated gradients therefore potentially extending the physiological relevance of the use of such cues in collective migration.”
Still, “Natural laboratory conditions” might also mean “close to physiological conditions” and the reviewer might wonder whether Dd cells actually encounter hypoxia in their normal conditions. In the natural habitat of Dd, the soil, depletion of nutrients or of oxygen is generally gradual and growth phase cells have mechanisms to sense when hard times are approaching by secreting factors (Maeda 2005). But, abrupt oxygen changes such as the one induced by our vertical confinement with the coverglass might also happens after a flood (Pedersen, Perata, and Voesenek 2017).
Regarding the question about the life stage of Dicty, we believe that this may indicate that the reviewer has in mind classical Dd development conditions after abruptly washing nutrients and putting cell in a starvation buffer. To approach that question and to reply as well to item #3 below (i.e., about the influence of cell division), we performed our confinement spot assay in a phosphate starvation buffer. A ring indeed develops quickly, in the same timing as in the nutrient buffer HL5 but it quickly dissipates as Dd cells stop dividing in these conditions. Still, this new experiment shows that ring formation happens during both the early development stage (after abrupt starvation) and during the growth phase. Description of this experiment was added at the end of the Potts simulation section of the paper and we added a new figure supplement and a new video to report its results.
Addition to the text: “First, we show in Figure 5B the effect of turning cell division off in the simulated spot. […] This still demonstrates that the phenomenon observed here is relevant for both the single cell and collective stages of Dd life cycle.”
The reviewer also questions whether the observed phenomenon is truly a collective one. We agree that the expansion of the ring is due to “many individual cells performing the same singlecell behavior”. However, we argue that ring formation is a response to a collective stress. The response is clearly collective as well and we can talk about an emergent behavior, because a single confined cell or even a too loose colony will not produce a ring. The ring is the consequence of the O_{2} gradient generated by the high cell density and the limited available O_{2}. Said differently, we agree that the “only interaction between cells (…) is through the oxygen levels” but these oxygen levels are themselves the result of the presence of a large population of cells. For these reasons, we have made the choice of describing this phenotype as a collective migration but specifying that it occurs without direct cellcell communication.
Finally, we also agree with the reviewer that our confined spot assay is a good model system to study the areotaxis of many cell types including cancer cells. Some of the authors developed it on epithelial cells (Deygas et al. 2018) and this manuscript demonstrates the robustness of the readouts (ring speed, shape…) for amoeboid cells. We believe that invasive cells might also be studied under this assay. We have added a short sentence to the end of the first Results subsection to underline this aspect:
Addition to the text: “This shows that the spot assay is an excellent experimental system to study the response of a variety of cell types to vertical confinement and its physiological consequences (Deygas et al. 2018).”
2) The authors first present a cellbased model, and then they study a "mean field" model that is based on the same underlying biological assumptions as the cellbased model. Interestingly, however, the cellbased model explains the data well, whereas for the meanfield model additional assumptions (a second threshold below which no aerotaxis occurs or alternative response functions) are required. What explains this discrepancy between the two models? We would expect that the meanfield model has the same overall behavior as the cellbased model upon which it is based, unless the assumptions differ or the discrete interactions between cells introduce new dynamics. For example, one could imagine that cell adhesion represented in the CPM strengthens the formation of the ring (i.e. reduces cell diffusion at higher cell densities). Could this or other "hidden" features of the CPM be introduced into the PDE and suffice to explain the discrepancy? What is the effect of the CPM parameters on ring formation?
The reviewer is correct in that Potts simulations and the elementary Go or Grow model led to different predictions, especially in terms of the shape of the advancing front. The referee is also correct in pointing out that both approaches are based on the same underlying assumptions. However, one goal of this elementary Go or Grow model was to remain analytically solvable so that we could have a more detailed understanding of key parameters than with Potts simulations. For that reason, the elementary Go or Grow model wasn’t a mean field representation of the Potts model as many details were different. For instance, the aerotactic response was set to be proportional to the local gradient in Potts simulations but based solely on its direction in the elementary Go or Grow model. The form of oxygen consumption by cells was also different in both approaches. Finally, our previous version of the Potts model didn’t have reduced cell division at low O_{2} concentrations unlike the elementary Go or Grow model. Still, we agree with the reviewer that these differences were not clearly put forward and that the discrepancies between the two models were thus problematic. Based on this comment, we have made numerous changes to both the Potts and the Go or Grow section of the paper. In summary:
– Potts simulations now incorporate the fact that cells stop dividing below an O_{2} threshold, a biologically sound hypothesis (Giampietro Schiavo and Bissons 1989; West, van der Wel, and Wang 2007), thus turning it into a Go or Grow model as well. Figures showing the results of these simulations have been updated accordingly.
Addition to the text: “Finally, all cells can divide as long as they sit in a high enough O_{2} concentration since it was demonstrated that cell division slows down in hypoxic conditions (G. Schiavo and Bisson 1989; West, van der Wel, and Wang 2007)”.
– The mathematical section of the paper now starts with a mean field representation of these Potts simulations. We show that this analytical approach can therefore recapitulate the numerical simulations when all ingredients are similar and hence recapitulate with a very satisfactory agreement the experimental data.
– Since this mean field approach is not solvable, we then turn to our original elementary Go or Grow model which gets further away from our experiments but allows us to make quantitative, general predictions in particular about the wave speed, its independence towards oxygen dynamics and the relative role of division and aerotaxis in setting the wave speed.
– We now end by showing variations on this Go or Grow model not in order to reproduce our experiments but to show the applicability of our method and to underline the fact that the analytical predictions are robust with respect to small modifications of the model, hence showing the general relevance of our results beyond the aerotaxis of Dd cells.
As a result, the mathematical section of the model has been largely rewritten and reorganized. Figures 6 and 7 have been modified to include this new structure. One issue we pointed out with the Potts model (formation of a dense core of cells) was removed as it no longer is true and this actually improves the quality of these simulations. Finally, Methods sections have been modified to reflect these changes.
We thank the reviewer for this comment as we believe it has made our approach more sound and easier to follow for the reader while insisting more on general results rather than simply reproducing the experiments. In that sense, these changes also partly address Comment #3 below.
3) The cellbased models and mathematical model produce interesting predictions as outlined in Figure 5 and in the expressions of wave speed for example. Yet the modeling is an end point. The models are compared with the goorgrow models and with the FisherKPP reactiondiffusionmodel, but the predictions are not used to further analyse the system, which makes the story not 'round' – the modeling seems to get ever more detailed until it reaches a sort of 'dead end'. Yet many of the predictions can be tested. E.g. divisions (Figure 5B) can be inhibited. The oxygen gradients can be manipulated. Motility (D) can be changed etc.
We agree with the reviewer that both our models led to interesting predictions that were not tested experimentally to further justify their generality.
We have now incorporated two such comparisons between the Potts simulations and experiments.
First, we performed spot assay experiments in phosphate buffer (see response to comment #1). This experiment not only demonstrate that ring formation also occurs under starvation conditions, it allows us as well to test the behavior of this ring in the absence of cell division (see figure 5—figure supplement 2). The agreement with the simulations (Figure 5B) is very satisfactory as the ring first propagated but quickly stopped after 3h and finally dispersed as cell density was no longer sufficient to maintain hypoxic conditions.
In addition, we also asked what would be the behavior of rings under more complex oxygen environments. As a first step, we asked what would be the behavior emerging from two colonies put in close proximity. Would the rings repel each other, attract each other, avoid each other? There too, we found good qualitative agreement between the simulations and experiments: these rings will fuse, form an oblong shape that will relax towards a circular one thus forming one large ring as a result. We have added a new figure supplement and a new video to document these two results. We have also added the paragraph below to the end of the Potts section of the paper to introduce these tests.
Addition to the text:
“Because of their versatility, they can also be used to make some predictions on the observed phenomenon. […] However, they fall short of giving an indepth quantitative description because they rely on many parameters and are not amenable to theoretical analysis per se.”
We also agree that our analytical efforts reached a ‘dead end’ and we believe this was due to our attempt at tuning this model to simply reproduce our experiments. As described in response to comment #2, we have now completely restructured this part of the paper and focused on the applicability of our results to different systems and experiments rather than focusing solely on our own.
– The mathematical section of the paper now starts with a mean field representation of these Potts simulations. We show that this analytical approach can therefore recapitulate the numerical simulations when all ingredients are similar.
– Since this mean field approach is not solvable, we then turn to our original elementary Go or Grow model which gets further away from our experiments but allows us to make quantitative, general predictions in particular about the wave speed, its independence towards oxygen dynamics and the key role of division and aerotaxis.
– We now end by showing variations on this Go or Grow model to show the applicability of our method and to underline the fact that the analytical predictions are robust with respect to small modifications of the model, hence showing the general relevance of our results beyond the aerotaxis of Dd cells.
We agree that other comparisons between predictions of the models and experiments would be beneficial but many are difficult to achieve and we feel go beyond the scope of this paper. In the near future, we intend to pursue the fine manipulation of oxygen gradients and oxygen levels in parallel to refine the aerotactic response of Dd cells.
4) The authors are making conclusions about "aerokinesis" because cells experiencing a gradient but behind the crest are highly motile. Another explanation is that the motile cells' behaviour is due to different (greater?) resolvable directional information (Tweedy and Endres, Scientific Reports, around 2014). Could the authors perform experiments to test whether cells at the "aerokinetic" O2 levels behave the same way if they are not experiencing a gradient? Can the authors keep cells in a homogeneous O2 at that level and see if they still polarise and move fast? If so it is indeed aerokinesis. Alternatively, the authors could discuss the issue.
We thank the reviewer for this remark and we agree that our claim about aerokinesis was based solely on data where cells were also experiencing an oxygen gradient. We thus followed the reviewer’s advice and checked for aerokinesis in homogenous environments.
We built a glass environmental chamber compatible with videomicroscopy studies at very low oxygen concentrations in particular with the random motility assay we have performed for many years in the laboratory (d’Alessandro et al. 2018; Golé et al. 2011). Our current chamber and gas blender enables the control of O_{2} down to 0.4%. Cell diffusion constant was very significantly higher at 0.4% (D=40.2±9.6 µm^{2}/min) than at the 20.95 % control condition (D=19.2±8.8 µm^{2}/min, Figure 2figure supplement 8 AD). We could not notice any difference between the diffusion constant of cells that stayed 3h or 20h at 0.4% (Figure 2—figure supplement 8) while control cell motility tends to decrease with time as quorum sensing factors accumulates in the medium (d’Alessandro et al. 2018; Golé et al. 2011). Noticeably, cell trajectories were more homogeneous at 0.4% than at 20.95%. A reason could be that cells almost stop dividing at 0.4% while the motility of control cells might be affected by cellcycle position. About cell shape and cell polarization, we did not image classical polarization markers but this is something we plan to do for future studies. However, the shapes of the cells in the hypoxic condition are more elongated than in normoxic conditions (not shown but similar to the shape of cells respectively at the left, i.e., hypoxic region, and right side, i.e., normoxic region, of the ring in Figure 1C). All together, we showed for the first time, to our best knowledge, an aerokinetic effect in Dd at very low O_{2} concentrations. This work offers very interesting perspective for future work to measure motility and phenotype changes as a function of oxygen level in the 0.1%2% range where we observe an aerotactic response. We have added a figure supplement to Figure2 to recapitulate these results along with a description of these experiments into the microfluidics section of the text.
Addition to the text: “Since cells in the microfluidic devices were also experiencing oxygen gradients, we further tested if the observed was true aerokinesis. […] The very significant difference (p<0.0001) demonstrates that Dd cells show an aerokinetic positive response to low oxygen, even in the absence of gradients.”
5) Around line 497 the authors make some strong conclusions about the generality of the system. But it seems they are all based on a strictly 2D system, and some of the most interesting conditions are nearly 1D (cells moving down a pipe) or 3D (real life colonies of amoebas or bacteria in soil; there's lots of 2D action in soil but O2 is likely freely available in most). Could the authors (a) try a 1D experiment (b) consider a 3D thought experiment or model? In any case they should discuss how the 1D and 3D cases differ and how important they might be. Maybe even semi3D (thickish layers?).
We agree that this comment deserves a short discussion in the manuscript and we thank the reviewer to point this lack. Our spot assay is indeed almost purely 2D as the vertical confinement thickness e~50µm is negligible as compared to radial directions. We also agree that in soil the geometry is probably hybrid with cells restricted to move on complex surfaces (grains or soil fibers) but with O_{2} diffusing from the whole surrounding space in 3D. As such, our experiments are clearly simplifications of the environment actually encountered by Dd cells in the wild.
However, our experimental measurements as well as our numerical simulations show that an important property of the migrating front is the fact that it follows a constant oxygen concentration within the dynamic gradients. This result doesn’t depend on the dimensionality of the problem and simply requires that hypoxic conditions are reached and that selfgenerated gradients are present, i.e. as long as cell density is high enough and oxygen sources are further away from the cells. For instance, a similar behavior in different geometries was already observed by Reviewer #2 using the under agarose assays (Tweedy and Insall 2020). The ~1mm thick agarose layer soaked with folate restricts convection without greatly affecting diffusion. The reviewer shows that selfgenerated gradients of degraded folate induce a group migration of cells in bands (in 1D) or rings (in 2D spots) up to 4 mm. This work, just as ours, points to the importance of a sufficient cell density in the initial colony to degrade enough of the attractant and of cell division to maintain this density throughout the migration over long distances. We therefore believe that our results can be generalized to any dimension and that knowledge of the oxygen profiles arising from cell consumption and diffusion from a source are sufficient to predict the existence and behavior of a migrating front.
In 3d: under spherical symmetry, we expect the formation and migration of a dense spherical ring of cells. The cell density required to achieve hypoxic conditions could however differ from 2d situations because diffusion will be more efficient at replenishing the oxygen consumed by the cells for geometrical reasons.
In more complex geometries, such as the soil, we believe it is difficult to speculate as to the distribution and diffusion of oxygen in such a complex environment but we believe that a dense enough colony to create hypoxic conditions will trigger the migration of a dense front whose geometry will depend on the shape of the complex oxygen profiles.
Still, to fully demonstrate the generality of our results, we have now added Potts model numerical simulations in 3d, one showing the formation of a dense sphere of cells when oxygen sources are present in all directions. In a second one, we model a slightly more complex geometry in which oxygen sources are present on all sides and on the top of the original colony but not on the bottom. We then observe the migration of a halfspheric front moving upwards towards the source but no front was observed moving downwards. Of note, for computational reasons, these simulations were made using a minimal Potts model that is no longer quantitative but possesses all the key ingredients described in the text. We have added a paragraph in the Discussion section addressing this role of dimensionality and introducing these new results which we present as an additional video.
Concerning the analytical results, under the assumption of spherical symmetry, the equation reduces to the 1D case when curvature effects can be neglected. This is actually the natural mathematical framework to construct traveling wave solutions which are stationary in the moving frame. Indeed, curvature transiently slow down the speed of propagation as compared to the 1D case. Therefore, the 1D analysis can be viewed as the long time asymptotic of planar front propagation beyond the curvature effects.
Addition to the text:
“In addition, although our experimental results were obtained on simple, 2d experiments, our findings can generalize to more complex cases. […] Our 2d results can therefore be extended to any other situations and they show that the key to proper steering are high enough cell densities and the creation of robust selfgenerated gradients.”
6) An open request, which would help future work to pick the molecular bones of this work – is for an overt description of the concentrations of O2 that provoke the strongest responses.
We thank the reviewers and editor for this suggestion to which we clearly adhere. We believe it also echoes one of the reviewers’ comment: “The combination of wet experiments, models and math is this paper's strength; bringing different informative figures together to make a clear conclusion. Its weakness to some will be a lack of description of the underlying pathways.”
In our opinion, our work calls for two main prolongations in the future:
– An indepth exploration of the response of individual Dd cells to more diverse oxygen environments. For instance, as described in the response to comment #3, we will now focus on developing tools to independently control oxygen levels and gradients to better study the crosstalk between these two aspects, crosstalk which we find to be central to the collective behavior we observed. For instance, it would be very informative to observe Dd cells in the 02% range of oxygen where they are aerotactic and tune the steepness of the gradient within this range to pinpoint the exact aerotactic modulation by local concentrations.
– A molecular characterization of the reponse of Dd cells to hypoxia and oxygen gradients. The spot assay which is easy to use experimentally could be used in this context as a screening tool for mutants.
Although we fully agree that these are natural steps stemming from this paper, we also feel like they both require long developments and/or experiments and that they would produce data which could be drowned in this paper that already touches on many different fields and aspects of the problem. Both would probably be better suited for systematic analysis.
References:
Adler, J. 1966. “Chemotaxis in Bacteria.” Science 153: 708.
Adler, Julius. 1966. “Chemotaxis in Bacteria.” Science 153(3737): 708–16.
https://pubmed.ncbi.nlm.nih.gov/4957395/ (December 9, 2020).
Budrene, Elena O, and Howard C Berg. 1991. “Complex Patterns Formed by Motile Cells of Escherichia coli.” 349(February): 630–33.
Chang, Weilun et al. 2017. “Asymmetric Distribution of HypoxiaInducible Factor α Regulates Dorsoventral Axis Establishment in the Early Sea Urchin Embryo.” Development 144: 2940–50.
Child, C. M. 1941. “Formation and Reduction of Indophenol Blue in Development of an Echinoderm.” Proceedings of the National Academy of Sciences 27(11): 523–28. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1078374/ (April 8, 2021).
Coffman, J.A., and J.M. Denegre. 2007. “Mitochondria, Redox Signaling and Axis Specification in Metazoan Embryos.” Dev Biol. 308: 266–80.
Cremer, Jonas et al. 2019. “Chemotaxis as a Navigation Strategy to Boost Range Expansion.” Nature 575(7784): 658–63. https://doi.org/10.1038/s415860191733y (April 8, 2021). d’Alessandro, Joseph et al. 2018. “Collective Regulation of Cell Motility Using an Accurate DensitySensing System.” Journal of The Royal Society Interface 15(140): 20180006. https://royalsocietypublishing.org/doi/10.1098/rsif.2018.0006 (June 5, 2020).
Deygas, Mathieu et al. 2018. “Redox Regulation of EGFR Steers Migration of Hypoxic Mammary Cells towards Oxygen.” Nature communications 9(1): 4545. http://www.nature.com/articles/s41467018069883 (December 21, 2018).
Golé, Laurent, Charlotte Rivière, Yoshinori Hayakawa, and Jean Paul Rieu. 2011. “A QuorumSensing Factor in Vegetative Dictyostelium Discoideum Cells Revealed by Quantitative Migration Analysis.” PLoS ONE 6(11).
Gregor, Thomas, Koichi Fujimoto, Noritaka Masaki, and Satoshi Sawai. 2010. “The Onset of Collective Behavior in Social Amoebae.” Science 328(5981): 1021–25.
Kelly, Beth et al. 2021. “Sulfur Sequestration Promotes Multicellularity during Nutrient Limitation.” Nature.
Maeda, Yasuo. 2005. “Regulation of Growth and Differentiation in Dictyostelium.” International Review of Cytology 244: 287.
Pedersen, O., P. Perata, and L. Voesenek. 2017. “Flooding and Low Oxygen Responses in Plants.” Functional Plant Biology 44: iii–vi.
Saragosti, J et al. 2011. “Directional Persistence of Chemotactic Bacteria in a Traveling Concentration Wave.” Proceedings of the National Academy of Sciences 108(39): 16235–40. http://www.pnas.org/cgi/doi/10.1073/pnas.1101996108.
Schiavo, G., and R. Bisson. 1989. “Oxygen Influences the Subunit Structure of Cytochrome c Oxidase in the Slime Mold Dictyostelium Discoideum.” The Journal of biological chemistry 264(13): 7129–34. http://www.jbc.org/article/S0021925818832112/fulltext (April 8, 2021).
Schiavo, Giampietro, and Roberto Bissons. 1989. “Oxygen Influences the Subunit Structure of Cytochrome c Oxidase in the Slime Mold Dictyostelium Discoideum ”.” c.
Tweedy, Luke et al. 2020. “Seeing around Corners: Cells Solve Mazes and Respond at a Distance Using Attractant Breakdown.” Science 369(6507). https://science.sciencemag.org/content/369/6507/eaay9792 (April 8, 2021).
Tweedy, Luke, and Robert H. Insall. 2020. “SelfGenerated Gradients Yield Exceptionally Robust Steering Cues.” Frontiers in Cell and Developmental Biology 8. /pmc/articles/PMC7066204/?report=abstract (July 9, 2020).
West, Christopher M., Hanke van der Wel, and Zhuo A. Wang. 2007. “Prolyl 4Hydroxylase1 Mediates O2 Signaling during Development of Dictyostelium.” Development 134(18): 3349–58. http://www.dictybase.org (April 8, 2021).
https://doi.org/10.7554/eLife.64731.sa2Article and author information
Author details
Funding
IFS (J19Ly11)
 Kenichi Funamoto
 JeanPaul Rieu
Université de Lyon (STARMAJ)
 Satomi Hirose
European Research Council (639638)
 Vincent Calvez
Agence Nationale de la Recherche (ANR19CE45000202)
 JeanPaul Rieu
GDR Imabio (AMI)
 Christophe Anjard
 JeanPaul Rieu
Centre National de la Recherche Scientifique (MITI 2019)
 Vincent Calvez
 JeanPaul Rieu
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank R Fulcrand for his expert help in microfabrication; G Torch, C Zoude, G Simon and A Piednoir for technical assistance. This study was supported by the CNRS  Mission pour les Initiatives Transverses et Interdisciplinaires – « Défi Modélisation du vivant  2019», by the GDR ImaBio (AMI fellowship) and by the IFS LyC Collaborative Research Project 2019 (Tohoku University). S Hirose was supported by the STARMAJ Program (FranceJapan: Research Internships for Master’s Students, Université de Lyon) and K Funamoto by the CNRS (invited researcher position). This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No 639638 to V Calvez) and from Agence Nationale de la Recherche (grant 2019, ADHeC project, No ANR19CE45000202 to JP Rieu).
Senior Editor
 Naama Barkai, Weizmann Institute of Science, Israel
Reviewing Editor
 Tatjana Piotrowski, Stowers Institute for Medical Research, United States
Reviewers
 Roeland MH Merks, Leiden University, Netherlands
 Robert Insall
Publication history
 Preprint posted: August 18, 2020 (view preprint)
 Received: November 9, 2020
 Accepted: July 30, 2021
 Version of Record published: August 20, 2021 (version 1)
Copyright
© 2021, CochetEscartin et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 1,273
 Page views

 101
 Downloads

 12
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Cancer Biology
 Computational and Systems Biology
Lung squamous cell carcinoma (LUSC) is a type of lung cancer with a dismal prognosis that lacks adequate therapies and actionable targets. This disease is characterized by a sequence of low and highgrade preinvasive stages with increasing probability of malignant progression. Increasing our knowledge about the biology of these premalignant lesions (PMLs) is necessary to design new methods of early detection and prevention, and to identify the molecular processes that are key for malignant progression. To facilitate this research, we have designed XTABLE (Exploring Transcriptomes of Bronchial Lesions), an opensource application that integrates the most extensive transcriptomic databases of PMLs published so far. With this tool, users can stratify samples using multiple parameters and interrogate PML biology in multiple manners, such as two and multiplegroup comparisons, interrogation of genes of interests, and transcriptional signatures. Using XTABLE, we have carried out a comparative study of the potential role of chromosomal instability scores as biomarkers of PML progression and mapped the onset of the most relevant LUSC pathways to the sequence of LUSC developmental stages. XTABLE will critically facilitate new research for the identification of early detection biomarkers and acquire a better understanding of the LUSC precancerous stages.

 Computational and Systems Biology
 Neuroscience
Humans make a number of choices when they walk, such as how fast and for how long. The preferred steady walking speed seems chosen to minimize energy expenditure per distance traveled. But the speed of actual walking bouts is not only steady, but rather a timevarying trajectory, which can also be modulated by task urgency or an individual’s movement vigor. Here we show that speed trajectories and durations of human walking bouts are explained better by an objective to minimize Energy and Time, meaning the total work or energy to reach destination, plus a cost proportional to bout duration. Applied to a computational model of walking dynamics, this objective predicts dynamic speed vs. time trajectories with inverted U shapes. Model and human experiment (N=10) show that shorter bouts are unsteady and dominated by the time and effort of accelerating, and longer ones are steadier and faster and dominated by steadystate time and effort. Individualdependent vigor may be characterized by the energy one is willing to spend to save a unit of time, which explains why some may walk faster than others, but everyone may have similarshaped trajectories due to similar walking dynamics. Tradeoffs between energy and time costs can predict transient, steady, and vigorrelated aspects of walking.