A critical-like collective state leads to long-range cell communication in Dictyostelium discoideum aggregation

Abstract

The transition from single-cell to multicellular behavior is important in early development but rarely studied. The starvation-induced aggregation of the social amoeba Dictyostelium discoideum into a multicellular slug is known to result from single-cell chemotaxis towards emitted pulses of cyclic adenosine monophosphate (cAMP). However, how exactly do transient, short-range chemical gradients lead to coherent collective movement at a macroscopic scale? Here, we developed a multiscale model verified by quantitative microscopy to describe behaviors ranging widely from chemotaxis and excitability of individual cells to aggregation of thousands of cells. To better understand the mechanism of long-range cell-cell communication and hence aggregation, we analyzed cell-cell correlations, showing evidence of self-organization at the onset of aggregation (as opposed to following a leader cell). Surprisingly, cell collectives, despite their finite size, show features of criticality known from phase transitions in physical systems. By comparing wild-type and mutant cells with impaired aggregation, we found the longest cell-cell communication distance in wild-type cells, suggesting that criticality provides an adaptive advantage and optimally sized aggregates for the dispersal of spores.

Fig 1. Multiscale model: From single cell—shape changes and chemotaxis to collective behavior.

(A) Schematics for membrane dynamics (left), intracellular cyclic adenosine monophosphate (cAMP) dynamics (center), and extracellular cAMP dynamics (right). (B) Single-cell “streaming" simulation in a box with periodic boundary conditions and a constant concentration of cAMP (i). Box dimensions are about 25 x 90 μm (the initial cell radius is assumed to be ∼15μm). Because of the small dimension of the box, the cell is just leaking, not pulsing, in order to avoid saturation of secreted cAMP. The simulation was repeated 12 times, and the average chemotactic index (CI) was calculated (ii). Error bars represent standard errors. Differences in CI_x and CI_y are statistically significant (p < 0.01), using a Kolmogorov—Smirnov test (KS test). (C) Cells solve “back-of-the-wave” problem. (i) A Gaussian wave (σ² ∼ 60μm) moves from right to left with a speed of about 300 μm/min [40]. At the peak of the wave, the cell emits a pulse of cAMP. After the firing, the cell enters a refractory period during which it can neither fire again nor repolarize. The cell generally moves to the right and hence does not follow the passing wave. (ii) CI in x and y as well as in left (negative x), right (positive x), up (positive y), and down (negative y) directions in order to discriminate between the directions of the incoming (right direction) and outgoing (left direction) wave. Simulations are repeated 12 times; shown are averages and standard errors. Box is about 60 x 105μm. CI in the right direction is significantly higher than CI in the other directions. (D) “Aggregation simulation.” (i) Four cells are simulated moving in a constant concentration of cAMP. At the beginning, cells are randomly distributed. (ii) Density correlation at the end of simulations is plotted for control cells without secretion (blue) and all cells leaking cAMP and one cell also emitting pulses of cAMP (red). The red line has a significant (p < 0.05, KS test) peak at a distance of about two cell radii, representing cell—cell contact. Simulations in this case are also repeated 12 times. Box dimensions are 75 x 75 μm. See Materials and methods for details on density correlation, Supporting information for a full explanation of the detailed model, and S1 and S2 Data for data and code from simulations, respectively. (E) Schematic showing cells represented as point-like objects with velocity vectors. Firing cells emit pulses of cAMP, and nonfiring cells secrete cAMP at a low constant leakage rate. Spatial cAMP profiles are derived from detailed model simulations. At every time point, cells are allowed two possible directions of movement in order to reproduce pseudopod formation at the cell front, with directions changing by ±27.5° with respect to the previous movement, corresponding to an angle between pseudopods of about 55° [41]. (F) Screenshot during streaming for n = 1,000 simulated cells (i). Red (yellow) points represent nonfiring (firing) cells. (ii) Spatial information versus time: simulations (blue) compared with experimental dataset 3 (green). Values were then normalized and shifted in time to facilitate comparison; see S3 Data for the numerical values.

Fig 2. Coarse-grained model leads to collective behavior.

(A) Screenshot of simulation for n = 1,000 cells at different time points: prestreaming (left), streaming (center), and after aggregation (right). Red (yellow) points represent nonfiring (firing) cells. See Supporting information for the full movie. (B) Kymograph of nonconnected C_nc and connected C_c directional correlations for simulation of n = 1,000 cells. Directional correlation profiles C(r) were calculated for every time frame and displayed depending on distance r. For cell size, we used the volume exclusion length of 3 μm. (see Materials and methods). (C) Correlation length versus time for different numbers of cells (n = 1,200, 1,000, 800, and 600). Data were smoothed with a moving average filter spanning ten consecutive frames. (D) Susceptibility χ plotted with respect to nearest-neighbor (NN) distance for different numbers of cells and with respect to time (inset). Nearest-neighbor distance was rescaled by the volume exclusion length (see Materials and methods). The peak in susceptibility becomes higher the larger the number of cells, and NN distances decrease accordingly. Profiles in the inset were smoothed with a moving average filter spanning ten points. (E) Comparison of correlation profiles for streaming phase (50-min time window). Connected correlations, calculated for different numbers of cells and normalized so that the correlation length was equal to one, were plotted as a function of distance in units of their respective correlation lengths. The four profiles collapse onto a single curve, independently of the number of cells. (F) Average correlation length versus neighborhood radius for different numbers of simulated cells. L corresponds to the size of the images (389 μm). ξ₀ represents the average correlation length during streaming phase (50-min window). Error bars represent standard errors. See Supporting information for a full explanation of the model and S4 and S5 Data for MATLAB code and data, respectively.

Fig 3. Collective behavior in experimental data for validating model.

(A) Cyan fluorescent protein (CFP) images of Dictyostelium aggregation of dataset 3. Images were taken after 4–5 h of starvation, when cells were still moving randomly, before initiating aggregation (left), during streaming phase (450 min after first image, center), and after aggregation (800 min, right). (B) Kymograph of nonconnected C_nc and connected C_c directional correlations for the movie in dataset 3. Distance r is expressed in units of average cell size (estimated after an ellipse was fitted to every cell contour and corresponding to the average of the minor axis, ∼10.7 μm). (C) Spatial information (blue) and susceptibility χ (green) of movie in dataset 3 as a function of time. The increase in spatial information denoting a more ordered image corresponds to the peak in susceptibility. (D) Correlation length ξ₀ as a function of time for the six movies. Curves were smoothed with a moving average operation spanning 20 time points for better visualization. Inset: comparison of cell number estimated from TRED images during the streaming phase for different movies. (E) Susceptibility χ as a function of rescaled nearest-neighbor (NN) distance and as a function of time (inset). Note that the height of peaks increases and that the corresponding rescaled NN distance decreases with the number of cells, as it does for simulations. Rescaled NN distance was computed by normalizing NN distance by the average cell size. In order to decrease noise, profiles in the inset were smoothed with a moving average spanning 20 time points. (F) Normalized C_c as a function of correlation lengths ξ₀ for different movies. C_c for every dataset was calculated as an average over 150 min of the streaming phase. Error bars represent standard errors. Similar to the simulated data, curves collapse for different numbers of cells when correlations are plotted as a function of distance in units of their respective correlation lengths. (G) Average correlation length versus neighborhood radius. L corresponds to the size of images (2,033 pixels,∼1.3 mm). ξ₀ represents the average of 150 min during the streaming phase. Error bars represent standard errors. Cell positions and tracking for the different experimental data are provided in S6–S11 Data, with actual data provided in S12 Data. Numerical results for the correlation analysis for panels D–F are provided in S13 Data and for panel G in S14 Data.

Fig 4. Role of criticality in aggregation of wild-type and mutant cells.

(A) Correlation length (during streaming) and spatial information (of final aggregate) for coarse-grained simulations of n = 500 cells for wild-type (WT) and modified cell types. Correlation length and spatial information are normalized with respect to WT (blue symbol). Modified-cell simulations were performed with uniform (radially symmetric) secretion of cyclic adenosine monophosphate (cAMP) (red), significantly increased sensing noise (10-fold increase in standard deviation compared to WT noise; green), enhanced cell—cell adhesion (light blue), and asynchronized secretion (random pulsing; black). (B) Corresponding correlation length and spatial information for experimental data from Fig 5 of [55], considering WT cells (blue) and protein kinase A (PKA) pathway mutants (with the asynchronous regA mutant in black and the phosphorelay intermediate protein [rdeA] mutant in green). (C) Screenshots show cell distributions at the end of the simulations from (A). Error bars represent standard errors in correlation length for an average in time of 50 min during the streaming stage. See Supporting information for a detailed explanation and S15 and S16 Data for numerical values for simulations and experimental data, respectively.