We develop a coarse-grained stochastic model for the influence of signal relay on the collective behavior of migrating Dictyostelium discoideum cells. In the experiment, cells display a range of collective migration patterns, including uncorrelated motion, formation of partially localized streams, and clumping, depending on the type of cell and the strength of the external, linear concentration gradient of the signaling molecule cyclic adenosine monophosphate (cAMP). From our model, we find that the pattern of migration can be quantitatively described by the competition of two processes, the secretion rate of cAMP by the cells and the degradation rate of cAMP in the gradient chamber. Model simulations are compared to experiments for a wide range of strengths of an external linear-gradient signal. With degradation, the model secreting cells form streams and efficiently transverse the gradient, but without degradation, we find that model secreting cells form clumps without streaming. This indicates that the observed effective collective migration in streams requires not only signal relay but also degradation of the signal. In addition, our model allows us to detect and quantify precursors of correlated motion, even when cells do not exhibit obvious streaming.
Eukaryotic cells frequently transduce external chemical gradients into directed cell migration [1], a phenomenon known as chemotaxis. Seminal work in the last few decades has identified components of the intracellular biochemical networks mediating cell response to external chemical gradients and found that responsive components such as the phosphoinositide lipids (PIPs), PI3K, and PTEN are highly conserved across cell types. In these efforts, our model organism (Dictyostelium discoideum) has been a useful source for discovery of network components and the development of quantitative models exploring plausible mechanisms for mediating directional sensing. Despite the vast similarities in gradient detection among D. discoideum and mammalian cells including neutrophils and neurons, D. discoideum chemotaxis displays a striking collective phenomenon not often found in other cell types where D. discoideum cells responding to the extracellular chemical signal cyclic-AMP (cAMP) tend to migrate in a head-to-tail fashion termed streams. In response to an external cAMP cue, D. discoideum cells synthesize and secrete cAMP relaying the initial signal to nearby cells. Many cell types, including neutrophils, macrophages, and epithelial cells, have potential signal relay loops, but they do not tend to migrate in streams in a standard chemotaxis assay. Building on previous work [2]–[5], we develop a minimalistic model for D. discoideum migration and signal relay in a linear gradient. Our model incorporates recent experimental measurements on cell migration persistence [6], independence of signal strength [5], and migration mechanism and lag in reorienting to signals [7]. We use the model to ask what aspects of the signal relay loop promote streaming. We find that a balance between fast secretion and degradation are needed to match experimental observations. To constrain the migration parameters, we measure the time autocorrelations and the fluctuations of the cell motion from our experimental systems and we suggest the possible use of these metrics to find evidence of signal relay in cells that do not display streams. Our efforts are motivated by recent experiments on D. discoideum, that show a notable visual distinction between cells that relay signals, and cells that both relay and degrade the signal. Wild-type cells, which emit cAMP and degrade cAMP, can form streams where cells are aligned head to tail, while mutant PDE1- cells that are unable to degrade cAMP form transient, aberrant streams that lead to clusters [8]. When food is plentiful, D. discoideum cells exist as single cells and chemotax towards the bacterial metabolic product folic acid. When food is removed, D. discoideum transitions from single cell to collective behavior - through the spontaneous secretion and detection of cAMP. The cooperative behavior of this spontaneous transition was found to follow Winfield synchronization and, the emergence of pulsatile, signaling centers is beautifully described in [5]. These pulses travel through a population of D. discoideum in spiral waves [9], [10]. Secretion of the extracellular phosphodiesterase (PDE1) is essential for the spontaneous transition [11]. Each pulse of external cAMP detected by cells results in an increase in gene expression promoting collective behavior [12], and after 4–6 hours of cAMP mediated development, cells begin to aggregate. In order to determine the essentials for chemotaxis and streaming separate from those needed for development, researchers often provide exogenous pulses of cAMP [12], [13]. From these studies, it has been found that cAMP secretion is essential for streaming, but not for chemotaxis. Cells lacking adenyl cyclase A, the enzyme primarily responsible for internal cAMP production during aggregation, will chemotax to cAMP without forming streams [14]. Development and chemotaxis to cAMP in cells lacking the gene for PDE1 can be rescued through periodic addition of partially purified PDE1. Cells lacking PDE1 secretion will chemotax to cAMP and form transient streams to a central source of cAMP, though in linear gradients, such as the under agar assays, the streams appear thicker than wild type [8]. Spontaneous aggregation by developed PDE1 null cells can be recovered with the addition of a uniform bolus of exogenous PDE1, though the bolus is insufficient to recover the spatial extent of the streams. Because we intend to examine a minimalistic model, we include continuous, local cAMP secretion and a constant background of cAMP degradation. The dynamics of the pre-aggregation stage of D. discoideum development was analyzed by Potel and Mackay [15], where they observed the motion of cells and calculated various dynamic quantities, such as the mean speed and the mean square displacement of cells and used Furth' s persistent motion model [16], [17] to explain their observations. Futrelle et al. [18] investigated chemotactic response to an external signal for early, middle and late developed cells for different duration and frequencies of cAMP pulses. In particular, the chemotactic index and the speed of the cells during development were analyzed, and significant timescales that define the dynamics were extracted, including the response time to a change in cAMP gradient which they estimated to be on the order of seconds. Gingle [19] measured the smallest cell density, above which collective motion occurs. Gingle and Robertson [20] showed that this limit density depends on the development time of the cells. The spontaneous emergence of traveling waves in a population of D. discoideum cells has attracted interest of the mathematics and physics communities and lead to the development of several computational models to test hypothesis for mechanisms involving signal transduction, signal relay, and gradient sensing. Pioneering work by Martiel and Goldbeter used a differential equation approach based on the receptor activation and desensitization dynamics [21] to explain the pulses of cyclic AMP. In addition to modeling the receptor dynamics, following models studied mechanisms in D. discoideum chemotaxis including wave propagation of cAMP signals in an inhomogeneous excitable medium [9], [22]–[25], directional sensing via receptor activation followed by further intracellular signaling [26]–[28], and physical forces that regulate cell-cell or cell-surface interactions [29]–[32]. Other models of chemotaxis focus on stochastic aspects of the cellular processes. These models discuss mechanisms including stochastic dynamics of directional sensing and speed control [2], [33]–[36], ″memory″ associated with membrane deformations [37]–[39], extension of new pseudopods conditional on the locations of existing ones [40]–[41]. Recent models of chemotaxis study the effects of noise due to fluctuations in receptor binding as well as the noise arising from subsequent internal responses [4], [42]–[46]. In the simplest models directional sensing is represented as stochastic dynamics of a single angular variable (which represents the density asymmetry of both the occupied receptors and further downstream processes such as regulation). Schienbein et al. [33] showed that the dynamics of the stochastic angle agrees very well with the directional sensing dynamics of granulocytes during galvanotaxis. The stochastic angle model was also implemented for D. discoideum chemotaxis by including receptor kinetics and chemical gradient steepness [4]. In this work we choose to capture the stochastic effects by associating the stochasticity of the previously described angular variable with the measured fluctuations in the direction of motion. The focus of our study is on modeling, simulating, and analyzing collective motion arising from chemotaxis and signal relay. While collective motion, chemotaxis, and signal relay have all been investigated before, this work focuses on collective behavior in the presence of a linear gradient without fluid flow. The linear, no-flow gradient geometry has been used in conjunction with Zigmond chambers and underagar assays but was cumbersome and often replaced with point sources, such as a micropipette, which leads to convergent cell trajectories even in the absence of signal relay. A linear gradient has been recently incorporated into a microfluidic system which can simultaneously monitor multiple gradient conditions and cell lines (using EZ-TAXIScan system (ECI, Japan) [47]). By monitoring many parallel conditions we are able to clearly analyze signal relay and differentiate different types of collective motion. It also allows us to validate metrics for detection of collective behavior that should be useful for the analysis of a number of other investigations of cell signaling that are starting to be carried out in this signal geometry. Linear gradients have been introduced for quantitative studies of gradient sensing, but recent work in microfluidics devices has been carried out in chambers with fluid flow which flushes out signal relay (e. g. , in Refs. [45], [46]). The controlled linear gradient allows us develop a quantitative phenotype for the onset of signal relay between cells. We are able to tune the relative strength of signal relay continuously, by varying the linear gradient strength. This allows us to measure collective behavior based on correlations between cell trajectories. We anticipate that our systematic studies will be valuable for a broad range of investigations of collective cell behavior. Indeed cell trajectories in such linear gradient chambers are starting to be collected to investigate signaling pathways that regulate chemotaxis in various types of cells (e. g. , D. discoideum [48], neutrophils [49], [50], eosinophils [51], and osteoclasts [52]).
The EZ-TAXIScan system uses an etched silicon chip to form separate channels for chemotaxis experiments in a linear geometry [47]. Each channel contains two buffer wells on opposite sides of a thin, terraced gap (microns long, mm wide and microns deep). Cells are gently pipetted into one well and allowed to settle to the glass surface. The opposite channel is filled with cAMP and diffusion sets a linear gradient in the channel within minutes. Cells, responding to the external signal enter the terraced region and travel microns towards the other side. Parallel to the edges of the terrace are small columns (microns long, microns apart) that set the vertical spacing, but provide little impedance to cell motion. If not modulated by cAMP or by PDE1 secreted by the cells, the imposed gradient stays constant at least for minutes [47], [51]. This type of setup provides a good signaling geometry for separating the effect of intercellular communication and an imposed gradient. Fig. 1A and Fig. 1B show time lapse images of wild-type cells and mutant cells under the influence of a linear (downward in the figures) cAMP gradient. At cells placed in a reservoir without cAMP begin to move into the chamber (at the top boundary in the figures). Although the cells are initially introduced uniformly in the horizontal direction (min panel of Fig. 1A and Fig. 1B), wild-type cells are attracted to each other and form streams (min panel of Fig. 1A), which in this example evolve to swirling groups (min. panel of Fig. 1A). The mutual attraction of the cells is due to the enzyme adenyl cyclase A (ACA) localized at the back of the cells [14]. ACA synthesizes intracellular cAMP, which diffuses into the extracellular medium. As shown in Fig. 1B, mutant cells (aca-), lacking ACA, do not exhibit collective motion and, throughout the time-course of the experiment, move without streaming or clumping in the direction of the external cAMP gradient. To analyze these observed migratory behaviors, we use a cell tracking algorithm to determine cell displacement vectors over a short time interval of the position of the center of the imaged intensity of each cell. We define a motion angle as the angle of a cell' s displacement vector with respect to the imposed cAMP gradient. Fig. 1C shows representative tracks of cells during chemotaxis (color coded according to real time). Fig. 1D shows the distributions of the angle for aca- cells, subject to four different external cAMP gradient strengths, increasing by a factor of from panel to panel. The spread of reflects the competition between noise and the ability of cells to sense and react to the gradient. Note that the width of the distributions first decreases with increasing gradient strength then decreases indicating an optimum. This finding agrees with observations of Fuller et al. [45], which shows that the chemotactic response is limited by external noise (noise due to receptor-ligand binding) for small local cAMP concentration and by internal noise (noise due to subsequent internal signaling) for higher local cAMP concentration. The distributions in Fig. 1D show that the cells do not always orient in the direction of the extracellular gradient. As discussed in [53] the gradient-sensing mechanism is stochastic with many sources of noise that can cause random deviation from the direction of the external gradient. Our data for the angular distributions suggest that above a threshold gradient the cell orientation is independent of the gradient strength. Below this threshold (e. g. , see the panel of Fig. 1D), the width of the distribution increases with decrease of the gradient [45]. In what follows we focus on the regime where the cell migration is less sensitive to the gradient strength. For several representative cells, Figs. 2A–C show the time autocorrelation of, where the angle brackets denote an average over time for cells that are located in the region between the cell exit plane and the mid plane of the gradient chamber (lower half of the panels in Figs. 1A and B, (the number of cells are, , and, respectively) ). The reason for restricting the averaging to the half of the chamber on the cell exit side is to eliminate any bias of the cell orientation angle distribution due to influence of the process of entry into the chamber. For small angles () the autocorrelation is. The variance of, , is plotted as a function of the distance from the starting point of the cells in Fig. 2D for the three different gradient strengths. In the next section we develop a model which estimates the level of the fluctuations in the displacement (dashed line in Fig. 2D). Previous studies on eukaryotic HaCaT cells highlight the dependence of velocity autocorrelations on two time scales [37]. Nevertheless, we see from Figs. 2A–C that can be well fitted to a dependence of the form parametrized by the single characteristic time. The fits for the average correlations for the individual gradient strengths are displayed in Figs. 2A–C. The single time scale, , is approximately constant over the two orders of magnitude in the external cAMP gradient strengths (min, min and min for Fig. 2A, Fig. 2B, and Fig. 2C). This time scale is roughly consistent with the dynamics of contractions of cells [31]. The characteristic size of eukaryotic cells is an order of magnitude larger than that of bacterial cells, and, in contrast with the sensing by bacterial cells, eukaryotic cells can sense the difference in chemoattractant concentration between the front and the back of a cell, thus detecting spatial gradients without moving. For D. discoideum, gradient sensing is accomplished via a G-protein coupled receptor and downstream signaling pathways [54]. Models of chemotaxis treating the cAMP signal transduction mechanism, including the biochemical details such as receptor desensitization [21] and adaptation [55], demonstrate the emergence of the experimentally observed cAMP waves. In the present paper our modeling approach will differ somewhat from past works (e. g. , Refs. [9], [22], [24], [56]) in that we seek a model that is simple enough that its relatively few parameters can be inferred from experiments, yet is still capable of capturing the distinctions between streams and clumps seen in our experiments on D. discoideum. We model cells as self-propelled soft disks of radius. For each cell we specify the location of its center and its orientation by the two-dimensional vectors and (by definition). We specify locations of the cells using a rectangular coordinate system, where the chamber in which the cells move is located in. In the experiment, the chamber boundaries, and, have perforations and are thus permeable to transport of cells and cAMP. The speed of each cell is assumed to be well-approximated as constant in time, independent of signal strength, in agreement with controlled chemotaxis experiments [6]. The cAMP concentration field is denoted. In the experiment the cells are deposited in a large reservoir (corresponding to in the model) where there is no externally injected cAMP. This experimental condition is modeled by a Dirichlet boundary condition on the cAMP concentration, at, and by introducing individual discrete cells at with a uniform flux cells per unit time per unit length in (each newly introduced cell' s orientation is initially in the). In addition, the experiment has an aqueous solution of cAMP in a large reservoir on the other side of the chamber (corresponding to in the model), and the cAMP concentration in this reservoir stays constant during the course of the experiment. This is modeled by a Dirichlet boundary condition at, , along with the removal of cells when they reach. We applied periodic boundary conditions in, such that and each cell that leaves the chamber at a lateral boundary, or at, reenters the chamber at the other end. Using these definitions, we propose the following minimal, agent-based model for cell motion in our experimental setup, (1) (2) (3) The first equation corresponds to the constant speed assumption. The second equation dictates that the unit vector specifying the cell' s orientation is attracted toward the direction of the vector, (4) with relaxation time. This relaxation time may be thought of as including both the chemically determined time for a cell to ‘perceive’ the gradient, as well as the time it takes the cell to mechanically turn its orientation. The first term in is a unit vector in the direction of the cAMP gradient. Note that, in accord with the observed similarity of the second, third, and fourth panels of Fig. 1D, this term is independent of the level of cAMP (i. e. , invariant to the transformation). The second term in is white noise, (5) The third term in is a repulsive ‘force’ modeling a soft two-body contact interaction between neighboring cells, (6) where is the region. In Eq. (6) we have taken the form of the repulsive force to decrease linearly with distance from the center of the cell. We have also tried other forms for the dependence of the repulsive force and found that no qualitative differences occurred. Szabo et al. [57] and Chate et al. , [58] discuss the effect of adding cohesive (i. e. , attractive) forces in modeling tissue cells. The parameter determines the strength of the repulsion force. Eq. (3) is the diffusion equation governing the evolution of the distribution of the cAMP density, with constant diffusivity [59]. The parameter is the cAMP secretion rate of a cell. The cAMP decays at a rate which can be spatially nonuniform and is approximately proportional to the concentration of the degradation enzyme phosphodiesterase PDE1 [25]. We introduce a degradation inhomogeneity suitable for our experimental setup in the following section. The cAMP degradation rate in Eq. (3) is meant to account for the presence of the cAMP degrading enzyme PDE1 with assumed proportional to the enzyme density. Since PDE1 is secreted by the cells themselves and then diffuses, we can expect that and hence are time and space dependent quantities obeying an equation similar to Eq. (3) for the cAMP density, but with the term analogous to the degradation in Eq. (3) omitted. In the interest of simplicity, for our minimalist model, we wish to circumvent a full time-dependent diffusion equation model for. Instead, we assume that a time-independent steady state that is homogeneous in is established for the (we show in Text S1 that this is justified for the conditions of our experimental setup). This corresponds to depending on but not and, . Furthermore, in steady state, the -averaged cell flux in the -direction must, by conservation of cell number, be independent of in the linear gradient chamber, and its value everywhere in the chamber must be the same as the cell injection flux at. In the simplest case, without clumps, the averaged density of cells in the external linear gradient region will thus be roughly uniform in and of the order of. Thus the averaged PDE1 density, satisfies a one-dimensional, time-independent diffusion equation of the form (7) Here we approximate as constant in and given by where is the production rate of the PDE1 per cell per unit time; is the diffusivity of the PDE1 and is approximately [60]. In addition, we will argue that the appropriate boundary conditions on the PDE1 density are at and. Solution of Eq. (7) with these boundary conditions leads to the model, (8) That is, varies parabolically in; , and has its maximum value in the center of the chamber, . In our numerical explorations we mostly use the model Eq. (8). We also note that in other experiments, depending on the experimental setup, may have different dependence on. For comparison, we repeated our numerical runs with the spatially constant form, where the numerical prefactor is chosen so that the total amount of PDE1 in is the same as for Eq. (8) (i. e. , is the same). The spatially constant form for was used in other models of D. Discoideum chemotaxis [9], [21], [22], [24]. The results (shown in Text S1) are qualitatively similar to the results presented here. We now outline how we motivate the use of the boundary conditions (more detailed quantitative justification is given in Text S1). In our experiments, cells are placed in the reservoir located in. The cells then rapidly sink to the bottom of the reservoir (). The reservoir has a vertical thickness that is more than times larger than the vertical thickness of the linear gradient chamber. The same dimensions apply for the reservoir in. The bottom glass surface () of the reservoir in extends into, where it forms the bottom plane of the linear gradient chamber and of the reservoir in. Cells that are on the bottom of the reservoir supply a source of cells for entry at into the linear gradient chamber. The cAMP-degrading-enzyme PDE1, secreted by cells in the reservoir are assumed to be transported vertically upward by small convection flows in the reservoir fluid into the vertically large region of the reservoir. In contrast, the distribution of the PDE1 emitted by the cells in is constrained to the much thinner vertical region defined by the chamber dimensions. Thus, in the linear gradient chamber the PDE1 density cannot be attenuated to low levels by spreading vertically. As quantitatively shown in Text S1, based on this consideration, the enzyme density in and is much less than in the interior of the chamber. This leads to our previously stated approximate boundary conditions, , used in obtaining Eq. (8). In order to systematically determine the essential dependence of the behavior of the model on its parameters, we introduce appropriate nondimensionalizations. We define the dimensionless spatial coordinates by and. The dimensionless time scale is defined as, and the dimensionless cAMP density is defined as. With the rescaled variables, the cAMP boundary conditions become, and. Additionally, the white noise is transformed to, where. The model equations with the rescaled variables and Eq. (8) for can now be written as (9) (10) (11) where, , , , and. The integral of the summation over the square, is the number of cells in the unnormalized square, and is roughly. In the situations we investigate is always large compared to unity. Thus the term roughly plays the role of a normalized density whose nominal value is one. With these normalizations, the parameters in our model are, , and. We wish to explore the variation of the system behavior as a function of parameters. This is clearly an impossible task to carry out for the full set of dimensionless parameters just listed. Thus we now seek to restrict our detailed considerations to just a few of these parameters whose influence is, we think, the most interesting. If we regard for the cells as fixed, then the parameter is dictated by the experimental setup. Experimentally, the typical cell speed and hence is observed to be roughly the same for wild type, and mutant cells [6], and we therefore take as fixed. The noise term will be fixed by the experimental observations (e. g. , Fig. 1D) which imply that it does not vary significantly across the different experimental conditions investigated (see Text S1). Thus we will keep, and fixed at the appropriate estimated values. Furthermore, we suspect that the qualitative behavior of the system will be insensitive to the precise value of so long as (the situation we are interested in). Thus our main numerical model explorations will focus on how the model behavior depends on and. We now further discuss our reason for interest in varying and. First, with respect to, in reference [8] a genetic perturbation to the cells results in mutants lacking the ability to produce the degradation enzyme PDE1 (but still emitting cAMP). In our model this corresponds to setting. In our numerical experiments we will explore a continuous dependence on, partly because is not well determined, but also to understand the difference between mutant cells that do not emit PDE1 (i. e. , pdsA-/PEC cells) and wild-type cells. We also suggest that it may be useful for future experiments to explore continuous dependence on PDE1 secretion rate (i. e. ,) which might be realized by introducing a mixture of wild-type and mutant PDE1- cells. Regarding variations of, we note that the secretion of cAMP from cells, is biologically inhibited for another type of mutant, the aca- cells. Also, in our experiments, we change the external concentration of cAMP, . The biological and chemical changing of the parameters, and, both yield change of. (Also, could be tuned by changing the reservoir cell density and hence, but we have kept constant in our experiments.) Aside from and the parameters we used in our simulations are summarized in Table 1. We assume that the cell parameters in this table (i. e. , , , , , ,) are the same for wild-type cells () and mutant cells (). In the absence of mutual attractions through cell' s secretion of cAMP, a Fokker-Planck version of Eqs. (1–6) can be solved analytically (see the Text S1), and in Eq. (5) can be determined by matching the analytical result to experimental observations of mutant cells. Also, we estimate as being of the order of determined from our experimentally observed time-autocorrelation of the orientation vector (Fig. 2A), where is defined at the end of the previous section. This time scale is comparable to the contraction rate of D. discoideum cells which in the work of Satulovsky et al. [31] is considered as the bulk relaxation time. We note that the real cells' secretion rates of cAMP and of PDE1 are not well quantified and can be varied by drug treatment or by the use of mutant cells. Thus we will regard and the PDE1-level-dependent parameter as variable parameters and investigate how the dependence of the collective cell dynamics depends on them. The model equations, Eqs. (1–6) are simulated numerically. Figs. 3A–3C show representative cell tracks for three different values of the normalized cAMP secretion rate. For all three of these cases is fixed at, which we estimate to be consistent with previous experimental measurements [21]. The color at a given point on a cell track in Figs. 3A–3C indicates the time that the cell making the track was at that point, where red corresponds to the beginning of the simulation and blue corresponds to the end of the simulation. Figs. 3D–3F show representative snapshots, where the position and the orientation of the cell is indicated by an ellipse (at normalized time for D, E, and F). In the top panels of Fig. 3 (Figs. 3A and 3D), the relative cAMP secretion rate is very small (i. e. ,). This regime mimics the aca- mutant cells, and our numerical results qualitatively agree with the experimental observations of aca- cells (cf. , min panel of Fig. 1B). For larger values of, and depending on, our numerical results can be classified under two main categories, streams (Fig. 3E) and clumps (Fig. 3F). At moderate (Fig. 3E) streams are evident. At higher, Fig. 3F shows that multiple clumps of cells form. From the corresponding tracks of cells shown in Fig. 3C, it is seen that the cells stay within the clumps and the clumps have almost no motion in the direction. To go beyond the visual comparison of our simulation results with our experimental observations, a quantitative description of the three modes of group cell motion described above (i. e. , uncorrelated motion, streams, and clumps) is desirable. We define the normalized mean progression, by, where the angle brackets denote an average of cells in the region between and, where (cf. , [61], [62]). We denote by the average of over, and we denote by the time average of taken over the last quarter of the simulation (). Another useful measure is the normalized averaged cell density, computed by averaging over the region to with and normalized so that. First, Fig. 4A shows the ensemble average of, denoted by, for the aca- cell experiment (in gray) and for a single model simulation (in black). The model parameters for the run are and, which correspond to the aca- mutant cells. To make a fair comparison, for the experimentally obtained we filtered out cells that move at a slower speed than what we considered in our model (i. e. ,). We calculate for a group of randomly selected cells in the region. Since our tracking algorithm cannot track all the cells available in this region, the experimentally obtained is represented by this ensemble average. To compare our experimental result to our numerical simulation results, we calculate from our simulation by sampling cells in the simulation so as to match the number of cells for which is experimentally calculated. We show in Figs. 4B and 4C how, and vary with the distance from the cell reservoir, , for the three values of used to obtain the cell tracks shown in Fig. 3 with fixed at the same value used for Fig. 3. In these plots, , and are averaged over several runs (this average is denoted by), where the error in the mean is shown by vertical error bars, which is calculated by the standard deviations of the runs divided by the square root of the number of runs. In the low regime (solid curves in Figs. 4B and 4C), corresponding to Figs. 3A and 3D, Fig. 4A shows that, saturates to in the upper half of the gradient chamber, , while Fig. 4B shows that is approximately uniform. The density profiles measured from the time lapse images (a rough estimate calculated from the image intensity) fairly agree with those obtained from our simulations. For PDE1- cells, our model suggests that the cAMP secretion levels are small compared to the wild-type cells exposed to the same imposed gradient. The density profiles measured from the time lapse images (a rough estimate calculated from the image intensity) fairly agree with those obtained from our simulations. For PDE1- mutant cells, our model suggests that the cAMP secretion levels are small compared to the wild-type cells exposed to the same imposed gradient. In determining the cAMP secretion rate we assumed same noise level compared to the wild-type cells. Therefore, in conjunction with findings from our model, our experimental observations suggest that the lack of degradation of external cAMP results in either reduced signal relay or increased noise level in gradient sensing (corresponding to receptor desensitization). The comparison and the details of the density estimate are shown in Text S1. As shown in Figs. 3B and 3E, for, streams emerge in the regime of moderate (plotted as the gray dashed curves in Figs. 4B and 4C). These streams start to aggregate in the upper half of the gradient chamber, and this results in a decrease in and a corresponding increase in. Compared to the low regime, the streams cause an increase in the cell density (the peak at). In the high regime (plotted as the black dashed curves in Figs. 4B and 4C), is even more peaked than in the moderate regime. This apparently leads to a peak in the cAMP density which leads cells to start aggregating in the lower half of the gradient chamber. Streams form close to the reservoir, where cells enter the gradient chamber. To form streams, newly entering cells acquire laterally (-directed) converging velocity components. Since the cell speeds are fixed at, this causes to decrease (see the region in Fig. 4B) and to increase. This apparently leads to a more localized secretion of cAMP, which overcomes the externally imposed cAMP concentration causing the clumping seen in Figs. 3C and 3F. In Fig. 4D the maximum in the region is plotted versus the corresponding. Each point in this figure is obtained from a single numerical run. The points are color coded with respect to the and used in the numerical run. Fig. 4D shows that points are clustered in two regions. The first region, where is large and is small [ (), ], corresponds to large clumps, while the second region, where is small and is large [ (), ], corresponds to the uncorrelated motion. The points between these two regions correspond to runs where cells form streams which either generate clumps (i. e. , points closer to the first region) or move through the region and leave the gradient chamber (i. e. , points closer to the second region). In our model there are two time scales, and (the cAMP degradation rate and the local cAMP production rate), and we explored their effects. Fig. 5 shows results for averaged over and (i. e. , the last quarter of the simulation), as well as over a large number of model simulations (). These averages are labeled in the figure. The top panel of Fig. 5A shows as a function of for. Fig. 5A shows that decreases as increases. In the region, where decreases fastest, streams occur, but clumps are rare (e. g. , Figs. 3B and 3E). The bottom panel of Fig. 5A is for a very small value of (), modeling mutant cells that cannot degrade cAMP. In this case we see that there is a sharp decrease in in the range. Below this range the simulations show roughly uniform cell density, while above this range clumps occur. Compared to the slow degradation regime, in the fast degradation regime (top panel of Fig. 5A) the streaming behavior is robust. In the slow degradation regime, the streams form for only a short period which is followed by formation of clumps. Recent experiments demonstrate that stream formation is impaired, if cells cannot degrade external cAMP [8]. Fig. 5B summarizes results for our simulations (color coded), as a function of (plotted on the horizontal axis) and (plotted on the vertical axis). The data in the top (bottom) panel of Fig. 5A corresponds to a horizontal cut through Fig. 5B at the arrow, (), on the vertical axis of Fig. 5B. Fig. 5B shows that the width of the range of, where streams occur, decreases as is lowered. Additionally, the onset of stream generation with respect to becomes smaller with decreasing.
Our model explains different observed modes of collective motion of motile cells. Our main new finding is that signal relay alone is not enough to arrange migrating cells into collectively moving streams. However, when the signal is not only relayed but also degraded, stable streams form. Our model is minimal, involving a relatively small number of potentially experimentally deducible parameters. Based on our numerical results, we suggest experiments where the transition between streaming and clumping can be experimentally tested by changing the effective values of our model parameters. One suggestion is that the value of can be effectively reduced by either mixing wild-type and PDE1- mutants or by changing the amount of PDE1 added during the PDE1- mutant cell development. The relaxation time, obtained from our experimental observations, is associated with the membrane retraction time scale. In addition, the time scale corresponding to the noise amplitude is associated with the formation time of pseudopods [63]. These parameters could be altered by adding drugs or changing the developmental procedures. For example, introducing a drug that inhibits the PI3 kinase severely reduces the pseudopod generation frequency [63] and hence both and. Additionally, recent studies show drastic change in the collective motion behavior of wild type cells when they are prepared over a longer development time [64]. In this case and are reduced in agreement with the observed reduction of stream formation [64]. Thus, we believe that our model can be utilized to quantify changes in the collective motion in response to modifications of cell characteristics. In our model, we have only focused on the extracellular cAMP dynamics given in Eq. (3) with the objective of reproducing the patterns in Fig. 1 with as few physical processes as possible. We modeled the motion of the cells according to the the dynamics of sensing the signal with the phenomenological equation Eq. (2). Models that include additional processes (not included in our model) are capable of explaining additional phenomena. E. g. , models of cAMP signal transduction including receptor desensitization [21] and adaptation [55] show the generation of experimentally observed cAMP waves including spiral waves [3], [9], [56]. In addition, the observed rotating vortex structure of the aggregates can be explained by other self-propelled particle models which allow cells to adjust their propulsive force [65]. In the future we plan to modifying our model to investigate the effect of dynamic cell-cell adhesion in stabilizing stream formation, and aggregation. Our model can be extended to include competition between the gradient steepness, , and the local cAMP concentration, , by modifying Eq. (4) and introducing a competition between the noise intensity and the concentration of the cAMP. A simple approach is to impose the following limits: For small local cAMP concentration, the noise (second term in Eq. (4) ) will have a higher effect in the directionality (i. e. independent random motion). In contrast, for high local cAMP concentration, the contribution from the noise to local cAMP concentration ratio should be small compared to the gradient steepness to local cAMP concentration ratio. When the model is extended to include this competition, we can define an organization time scale as a measure of cellular organization. Thus, we can measure the efficiency of stream formation not only with respect to signal relay but also with respect to the efficiency of directional sensing. We believe that our simplified approach, used here for D. discoideum can be extended to more complex cells exhibiting signal relay, such as neutrophils [49], [66]. For neutrophils, signal relay is less well understood [49]. However, our numerical simulations can be utilized to distinguish uncorrelated motion from weak signal relay: Using our simulations in conjunction with linear gradient experiments where cells do not converge naturally to an external signal, we can calculate the effect of signal relay on the mean progression speed, as well as the development of an inhomogeneous density due to cell-cell attraction, even in the case of very small signal relay that is not sufficient to lead to discernible clumps or streams. Moreover, our model can be potentially extended to include the dependence of signal relay on cell density, in order to compare the dynamics to those observed in Ref. [67], which proposes a quorum sensing mechanism that can quantify the persistent random walk of D. discoideum at different phases of development as well as different densities. Another potential use of our model is to model migration when subpopulations of cells have different signal sensing, and signal relay capabilities. A prominent example of such collective migration is the motion of neural crest cells, a collective process during embryonic development. Recent experiments suggest that mathematical models of the neural crest migration require subpopulations having different chemotactic responses [68].
To examine the chemotactic dose response, cell migration was recorded at second intervals for 1 hour in the EZ-TAXIScan chamber (Effector Cell Institute, Tokyo, Japan). In the absence of wild-type cells the device establishes a well-defined, stable cAMP gradient during the course of the experiment [47]. Dictyostelium discoideum cells, wild-type cells (ax3) and its ACA null mutant cells (aca-) were prepared as described previously in Ref. [6]. PDE1- cells were prepared as described previously in Ref. [8]. There are two modules in our numerical simulation code, the first module consists of the equations of motion given in Eqs. (1) – (3) which defines the position and the direction of motion of cells based on the local gradient in the neighborhood of each cell. The second module calculates the diffusive time evolution of cAMP due to the external signal and dynamic local intercellular signals and provides the updated gradient vector field for use in the first module. Simultaneous evaluation of these two modules generates cell tracks. The diffusion equation (Eq. (3) ) for the cAMP is solved explicitly on a square grid with spacing m using a forward time and central space Euler method. In the numerical simulations the time step is seconds, which is well in the stable range of the numerical algorithm. For implementing the numerical evaluation of the diffusion equation is discretized with and. The Laplace operator can be replaced by the discretized Laplace operator and the Dirac- function is discretized in one dimension as, where is the Kronecker- function, it is zero except for. Thus, the value of the cAMP field at and, where and are integers, is updated according to (12) with. In Eq. (12), rounds its argument to the nearest integer. The same is used in evaluating the equations of motion (Eqs. (1) and (2) ). Table 1 shows the definitions and values of the parameters used in the numerical simulations.
Collective cell migration is observed in various biological processes including angiogenesis, gastrulation, fruiting body formation, and wound healing. Dictyostelium discoideum, for example, exhibits highly dynamic patterns such as streams and clumps during its early phases of collective motion and has served as a model organism for the study of collective migration. In this study, facilitated by experiments, we develop a conceptual, minimalistic, computational model to analyze the dynamical processes leading to the emergence of collective patterns and the associated dependence on the external injection of a cAMP signal, the intercellular cAMP secretion rate, and the cAMP degradation rate. We demonstrate that degradation is necessary to reproduce the experimentally observed collective migration patterns, and show how our model can be utilized to uncover basic dependences of migration modes on cell characteristics. Our numerical observations elucidate the different possible types of motion and quantify the onset of collective motion. Thus, the model allows us to distinguish noisy motion guided by the external signal from weakly correlated motion.
physics
developmental biology
cell motility
mathematics
cell migration
biology
nonlinear dynamics
morphogenesis
biophysics
Modeling and Measuring Signal Relay in Noisy Directed Migration of Cell Groups