A Critical-like Collective State Leads to Long-range Cell Communication in Dictyostelium discoideum Aggregation
AA Critical-like Collective State Leads to Long-range CellCommunication in Dictyostelium discoideum Aggregation
Giovanna De Palo , Darvin Yi ¤ , Robert G. Endres , Department of Life Sciences, Imperial College, London, United Kingdom Centre for Integrative Systems Biology and Bioinformatics, Imperial College, London,United Kingdom Joseph Henry Laboratories of Physics, Princeton University, Princeton, New Jersey,USA Lewis Siegler Institute for Integrative Genomics, Princeton University, Princeton, NewJersey, USA ¤ Current Address: Department of Radiology, School of Medicine, Stanford University,Stanford, California, USA* [email protected]
Abstract
The transition from single-cell to multicellular behavior is important in earlydevelopment but rarely studied. The starvation-induced aggregation of the socialamoeba Dictyostelium discoideum into a multicellular slug is known to result fromsingle-cell chemotaxis towards emitted pulses of cyclic adenosine monophosphate(cAMP). However, how exactly do transient short-range chemical gradients lead tocoherent collective movement at a macroscopic scale? Here, we use a multiscale modelverified by quantitative microscopy to describe wide-ranging behaviors from chemotaxisand excitability of individual cells to aggregation of thousands of cells. To betterunderstand the mechanism of long-range cell-cell communication and hence aggregation,we analyze cell-cell correlations, showing evidence for self-organization at the onset ofaggregation (as opposed to following a leader cell). Surprisingly, cell collectives, despitetheir finite size, show features of criticality known from phase transitions in physicalsystems. Application of external cAMP perturbations in our simulations near thesensitive critical point allows steering cells into early aggregation and towards certainlocations but not once an aggregation center has been chosen.
Author Summary
Cells are often coupled to each other in cell collectives, such as aggregates during earlydevelopment, tissues in the developed organism, and tumors in disease. How do cellscommunicate over macroscopic distances much larger than their typical cell-cellneighboring distance to decide how they should behave? Here, we develop a multiscalemodel of social amoeba, spanning behavior from individuals to thousands of cells. Weshow that local cell-cell coupling via secreted chemicals may be tuned to a critical value,resulting in emergent long-range communication and heightened sensitivity. Hence,these aggregates are remarkably similar to bacterial biofilms and neuronal networks, allcommunicating in a pulse-like fashion. Similar organizing principles may also aid ourunderstanding of the remarkable robustness in cancer development.
PLOS a r X i v : . [ q - b i o . CB ] J a n ntroduction Many living systems exhibit collective behavior, leading to beautiful patterns found innature. Collective behavior is most obvious in animal groups with clear advantages interms of mating, protection, foraging, and other decision-making processes [1, 2].However, how cells form collectives without visual cues is less understood [3]. There aretwo main strategies to achieve synchrony (or long-range order) among individuals:firstly, a leader, i.e. a special cell or an external field, may influence the behavior of theothers in a hierarchical fashion (top-down). An example is the developing fruit-flyembryo in maternally provided morphogen gradients [4, 5]. Secondly, all individuals areequivalent and order emerges spontaneously by self-organization (bottom-up). Examplesmay include organoids [6] and other cell clusters [7]. While order itself cannot be used todifferentiate between the two mechanisms, the response to perturbations or simply thecorrelations among fluctuations can be examined [8]. In top-down ordering fluctuationsare independent, while in bottom-up ordering fluctuations are correlated [9].To test these ideas of achieving order, we consider the well-known social amoeba
Dictyostelium discoideum , which undergoes aggregation in response to starvation [10–12].During this developmental process, cells start to secrete pulses of cAMP, a moleculethat also acts as a chemoattractant for the other cells in the vicinity. When a cell is ‘hit’by a high concentration of cAMP, it secretes a pulse of cAMP itself, relaying the signaland thus causing the formation of cAMP waves, inferred indirectly from optical densitywaves in dark field movies [13, 14]. These waves propagate through the wholepopulation [15–18]. As their development proceeds, cells pulse at higher frequencies,reaching a maximal frequency of about one pulse every six minutes [19]. Cell movementalso accompanies the secretion process: before cells start to secrete cAMP, they normallymove incoherently; when cAMP waves form, cells move towards the direction of theincoming wave by following the cells emitting the pulse in an orderly fashion (streamingphase). Interestingly, cells do not follow the passing wave in microfluidics (and hencesolve the ‘back-of-the-wave’ problem) [20, 21]. While single-cell chemotaxis [20–25] andlarge-scale pattern formation [14, 26–30] have been extensively studied, a precisecharacterization of the transition from single cells to multicellularity is still missing.Here, we develop a multiscale model to capture the mechanism of aggregation,focusing on the distinction between induced and self-organized order. Specifically, weare able to unify single-cell behavior and multicellularity at wide-rangingspatio-temporal scales. We achieve this by extending a single-cell model, which is ableto describe
Dictyostelium cell shape and behavior [22], by adding intracellular cAMPdynamics, secretion, and extracellular dynamics for cell-cell communication. Tosimulate hundreds of cells, we extract a set of minimal rules for building acoarse-grained model. Hence, our approach is able to capture all stages of aggregation,ranging from single-cell chemotaxis to the multicellular collective. For quantifying thetransition from disorder (pre-aggregate) to order (aggregate), we employ spatialinformation and directional correlations. We found that the transition occurs during thestreaming phase, which resembles a critical-like point known from phase transitions inphysical systems using finite-size scaling arguments. Criticality and other predictionsare tested by corresponding analysis of time-lapse movies from fluorescencemicroscopy [19, 31, 32], pointing towards universal behavior in cell collectives.
PLOS esults
A single-cell model fulfills criteria for aggregation
To model the transition from single cells to multicellularity, we started with cell shapeand behavior in single cells. Specifically, we considered a model capturing single-cellmembrane dynamics similar to the Meinhardt model [22, 23, 33]. This model describesmembrane protrusions and retractions, as well as resulting cell movement by means ofthree equations (see
Supporting Information ). The first and second variables are a localactivator and a global inhibitor (both are also considered in the local-excitationglobal-inhibition -LEGI- model [24, 25]). The third is a local inhibitor, important todestabilize the current pseudopod and to increase the responsiveness of the cell (Fig. 1A,left). To this we added dynamic equations representing the intracellular accumulationand secretion of cAMP from the cell rear [34, 35] based on the excitableFitzHugh-Nagumo model (Fig.1A, middle) [31, 32], as well as the extracellular dynamicsfor cell-cell communication (Fig. 1A, right; see Materials and Methods for furtherinformation and
Supporting Information for numerical implementation).Using this detailed model, we investigated the resulting behavior of the cell-cellinteractions in very small systems. First, we would like our model to capture streaming,i.e. the ability of a cell to precisely follow the cell in front of it. To reproduce that, wesimulated a single cell in a rectangular box with periodic boundary conditions (seeFig. 1B and supporting movie S1). Given the rectangular shape of the box, ahorizontally moving cell can sense its own secretion due to the periodic boundarycondition. In contrast, a vertically moving cell is too far away from its rear and thuscannot sense its secretion. We estimated the ability of the cell to stream by measuringthe chemotactic index (CI) in the x direction, calculated as the amount of movement inthe horizontal direction compared to the total length of the trajectory. In Fig. 1B, weshow that the CI in the x direction is significantly higher than the CI in the y direction.We then considered the wave response as measured in microfluidic experiments, inwhich cells are exposed to traveling waves of cAMP [20, 21]. When hit by a travelingwave, cells moved towards the direction of the incoming wave but did not follow thewave after it passed. In order to capture this behavior, our model cell undergoes arefractory period during which it cannot repolarize. This is achieved naturally as thecell spontaneously emits a pulse of cAMP when encountering the wave (see Fig. 1C, andsupporting movie S2). As a result, the CI is significantly higher in the right direction ofthe incoming wave. Finally, we considered a small number (four) of cells in a small box(with periodic boundary conditions) and tested whether they show signs of aggregation(see Fig. 1D and supporting movie S3). Specifically, we measured the density paircorrelation (see Materials and Methods), and compared the cases with and withoutsecretion of cAMP. In the absence of secretion, cells were randomly distributed in spaceat the end of the simulations. With secretion, cells tended to be much closer to eachother, with a clear peak in the density distribution at cell-contact distance (two cellradii). A coarse-grained model reproduces collective behavior
In order to reproduce aggregation as observed in experiments, e.g. [19, 31], we need tosimulate hundreds to thousands of cells. However, the detailed model introduced in theprevious section is computationally too expensive, forcing us to introduce severalsimplifications. In our coarse-grained simulations, cells are point-like objects moving incontinuous space. In particular, we took advantage of the spatial cAMP profiles fromthe detailed model by extracting the concentrations typically secreted by a single cellduring leakage or a pulse, shaped by degradation and diffusion. As in the detailed
PLOS
Supporting Information for a detailed explanation of the model).Using this minimal set of rules, we simulated thousands of cells with a densitysimilar to experimental ones (around a monolayer -ML- [19, 31]). Cells were initiallydistributed uniformly in space and allowed to move randomly. As soon as the celldensity increased spontaneously (and hence cAMP due to leakage by all cells), a cellmay sense a concentration of cAMP large enough to pulse and this excitation willpropagate throughout the whole population. Due to cell movement, streaming andaggregation into a small number of clusters can be observed (Fig. 1F and supportingmovie S4). To quantify aggregation in a different way, we also calculated the spatialinformation for capturing the order in an image based on the 2D Shannon information,relying on Fourier coefficients while not requiring tracking of individual cells (seeMaterial and Methods) [38]. In all simulations, this spatial information rises sharplyduring the streaming phase as expected for cells in an ordered aggregate (see Fig. 1G).Interestingly, the spatial information was previously used to capture the second-orderphase transition in the 2D Ising model [38]. Hence, we wondered whether aggregationmay be viewed as a critical-like point?
Collective behavior: hierarchical or self-organized?
Based on our model assumptions, all cells are treated the same. However, aggregationmay still be driven by the first random cell pulsing (hierarchical system) orspontaneously emerging as cells are coupled to each other by cAMP sensing andsecretion (self-organized system; Fig. 2A). The order of the collective process can bemeasured studying the directional correlations of pairs of cells. Specifically, the non-connected ( nc ) correlations C nc ( r ) = (cid:80) Ni (cid:54) = j −→ u i · −→ u j δ ( r − r ij ) (cid:80) Ni (cid:54) = j δ ( r − r ij ) (1)represent the average similarity of the direction of motion for every pair of cellsdepending on their distance, where N is the total number of cells, −→ u i is the vectorrepresenting the direction of cell i , and δ ( r − r ij ) is equal to 1 if r = r ij and 0 otherwise. C nc ( r ) also represents the order parameter in our system. By calculating this quantityfor every time frame, we can analyze its variation in time. During the preaggregationstage correlations are close to zero even at short distances, while they increase sharplyduring the streaming phase (Fig. 2B, top).Although order increases during the streaming phase, the origin and characteristicsof this order are yet to be determined. To achieve this, we calculated the connected ( c )directional correlations C c ( r ), measuring the similarity of the directional fluctuationswith respect to the average velocity [8, 9]. Thus, in this case direction −→ u i in Eq. 1 issubstituted by the velocity of the single cell when the average is subtracted, i.e. δ −→ u i = δ −→ v i / (cid:113) N (cid:80) Nj =1 δ −→ v j with δ −→ v i = −→ v i − N (cid:80) Nj =1 −→ v j . For this kind of collectivemovement, such a subtraction is not straightforward. If we compute a global average velocity for every time frame, we systematically overestimate the non-connectedcorrelations, because we still consider part of the “bulk” velocity vectors due to theposition of the cells in the image (see Supporting Information for a schematic
PLOS local averages . For every cell, weconsider the average velocity of all cells in its neighborhood up to a certain maximaldistance r c , and we compute the correlations between the cell in the center and all thecells belonging to its neighborhood, and repeat this procedure for every cell in ourimage.When applied to the simulations, Fig. 2B, bottom, shows significant connectedcorrelations, especially during streaming. Next, we considered the correlation length ξ as the minimum distance at which the correlation crosses zero, i.e. C ( r = ξ ) = 0 [8].We find that ξ is indeed much larger than the minimum nearest-neighbor distance,strongly suggesting self-organization (Fig. 2C). Streaming as a critical-like point
Above, we demonstrated that aggregation in
Dictyostelium is highly ordered andself-organized, with a correlation length much greater than the nearest-neighbor distance.Does the transition from disorder to order in this finite system show signs of criticality?In order to answer this question we considered that in critical systems thecorrelation length should scale with the size of the system as there is no intrinsic lengthscale [8]. To investigate this, we analyzed how the correlation length ξ changes in time.In all simulations ξ was small before aggregation and it increased markedly during thestreaming phase (Fig. 2C). The susceptibility can thus be computed as the maximumvalue reached by the integrated correlation χ = 1 N N (cid:88) i (cid:54) = j δ −→ u i · δ −→ u j θ ( ξ − r ij ) , (2)where θ ( ξ − r ij ) is equal to 1 for r ij < ξ and 0 otherwise [8]. This proxy for thesusceptibility peaks precisely during the streaming phase (Fig. 2D, inset), and thehigher the number of cells, the higher the susceptibility. Moreover, if we consider celldensity as a control parameter (similar to temperature or coupling in a ferromagneticIsing model), we can plot χ with respect to the rescaled nearest-neighbor distance(Fig. 2D and Materials and Methods). The resulting peaks do not only follow thenumber of cells in terms of height, but also shift to smaller nearest-neighbor distances asthe number of cells increases, further supporting the resemblance to a scale-free systemnear criticality [8]. Furthermore, normalizing the correlation and rescaling the distanceby the correlation length, the correlations collapse for all our simulations whenconsidering the average profile during the streaming phase (Fig. 2E). Finally, we takeadvantage of our image partition with different radii r c to examine how the correlationlength ξ scales with system size. We notice that for all movies, higher cell numbersdisplay longer correlation lengths for a given neighborhood radius, and that thecorrelation length increases with increasing radius. As a result, the correlation lengthscales with system size (Fig. 2F), indicating critical-like behavior in our simulated cells. Analysis of time-lapse fluorescent microscopy
To test the model, we analyzed five movies of
Dictyostelium aggregation with differentcell densities as described previously (see Materials and Methods and supportingmovie S5) [19, 31]. Briefly, during 15 hours of observation, individual cells become asingle, multicellular organism, going through different stages including preaggregation,streaming and aggregation (see Fig. 3A). Cell densities ranged from 1/3 ML to almost 1ML, ensuring aggregation while restricting our system to 2D. A 10% subpopulation ofcells expression TRED fluorescent marker were tracked using a custom-written software(see Materials and Methods). Based on these cells, we repeated the analysis from the
PLOS
Supporting Information ). Hence, cell numbers reported refer to the streaming phase(Fig. 3D, inset). Based on our analysis, the correlation length ξ increases during thestreaming phase as does the susceptibility χ (Fig. 3, B-E). Additionally, χ increaseswith cell number (and hence cell density), and the nearest-neighbor distance decreases,similar to the simulations. The correlation profiles, normalized and rescaled by thenearest-neighbor distance, largely superimpose for the different cell numbers, indicatingthat the slope of the resulting curves is not affected by the number of cells (see Fig. 3F).Finally, we studied how the correlation length changes for different system sizes byconsidering different neighborhood radii as performed for the simulations (see Materialsand Methods). We noticed that ξ increases for a given radius with increasing cellnumbers, and also for a fixed number of cells with increasing neighborhood radius(Fig. 3G). These observations strongly suggest that there is no intrinsic correlationlength, but that this length scales with system size. Taken together, our results suggestthat aggregation can be viewed as a critical-like point in this finite system. Additional model predictions and cell ‘steering’
Here, we consider the dynamics of aggregation more closely. When increasing thedensity of cells in our simulations, we noticed that cells aggregate faster at higher celldensities as measured by the slope of the spatial information (Fig. 4A). This increase inspeed appears to reflect the increased ability to relay the signal as nearby cells canbecome excited and pulse themselves, facilitating aggregation. To test this prediction,we attempted to quantify the speed of aggregation in our time-lapse movies as well.Although the experimental movies are much noisier and variable, we noticed a similartrend as in our simulation (Fig. 4B; only the dark blue curve violates the trend).Critical-like behavior is the tipping point between order and robustness on one sideand disorder and chaos on the other side. This point may highten the sensitivity of thecollective to detect changes in external cAMP concentration and to help cells make adecision on where to aggregate. While such experiments are beyond the scope of thiswork, we nevertheless attempted to investigate this with our simulations. For thispurpose, we applied both local and global cAMP perturbations to our system of cells;local perturbations represent a short local pulse of cAMP released on the 2D surface,while a global perturbation represents an overall pulse of cAMP on the whole area ofobservation (see
Supporting Information for details and supporting movies). Inparticular, we attempted to apply such perturbations before and during streaming, witha noticeable different effect on aggregation. A local perturbation prior to streaminginduces aggregation at the location of the pulse, although the aggregation centre mayshift later during the simulation. A global perturbation prior to streaming has a similaraffect, although the location of aggregation cannot be influenced. In contrast, applyingsuch perturbations during streaming has largely no effect. While local perturbation maysomewhat influence the location of aggregation, the overall dynamics stay the same.Indeed, comparing the time courses of the spatial information, we notice that the shapeof the curve always stays the same, and only the onset of aggregation can be shifted toearlier times with local or global cAMP perturbations (Fig. 4C). In particular, earlylocal perturbations allow the possibility of steering cells to locations of aggregation.
PLOS iscussion
Dictyostelium aggregation represents a fascinating example of synchronous collective cellbehavior, spanning ∼ ∼ µ m in size. Here, weasked how cells produce such long-range communication [39], when the transition fromsingle cells to the collective occurs, and how this transition can be characterizedquantitatively. To capture the main features of aggregation, we developed a multi-scalemodel. First, we focused on single cells using a detailed model combining sensing,cell-shape changes and movement with cAMP secretion/pulsing and hence cell-cellcommunication. Once this model resembled the behavior of a single cell or a smallgroup of cells, it allowed us to extract a minimal set of rules that could lead toaggregation. In particular, we extracted the cAMP concentration profile of a pulse fromthe detailed simulations and the refractory period after pulsing. By allowing cells toleak cAMP and to randomly move below a certain cAMP threshold concentration, wewere able to observe spontaneous random pulsing as soon as the local density increased,similar to what occurs in real cells. This minimal set was subsequently included in thecoarse-grained agent-based model, which is able to reproduce the collective behavior ofhundreds of cells in line with time-lapse microscopy [19, 31].Our major findings point towards previously uncharacterized features in aggregation,both observable in simulations and data. First, the transition to the collective is exactlypinpointed by a sharp rise in the spatial information of the cells during streaming.Second, to quantify the nature of the transition, we used fluctuations around the meanvelocity, allowing us to distinguish between a hierarchically driven, top-down (externalgradient from leader cells) and an emergent, self-organized, bottom-up (all cells areequal) process. Third, similar to second-order phase transitions in physical systems, thestreaming phase shows signatures of criticality using finite-size scaling arguments. As aresult there is no intrinsic length scale, allowing cells to communicate with each otherover large distances ‘for free’, i.e. only based on local cell-cell coupling. The controlparameter is cell density, affecting the cell-cell coupling via cAMP secretion and sensing.Our work provides further insights into the process of cell aggregation. By means ofour multi-scale model, we were able to answer why cells emit cAMP in pulses. Albeitshort-lived, a pulse creates a steeper spatial cAMP gradient than continuous secretion(assuming that the total amount of emitted cAMP is the same in both cases). Moreover,we noticed that so-called cAMP ‘waves’ are likely not actual macroscopic travelingwaves due to strong dissipation and diffusion. In contrast, cells are exposed toshort-range cAMP pulses, which need to be relayed from one to the next cell before theydissipate. Although cAMP waves from microfluidic devices were used to study thecellular response to positive (incoming wave) and negative (passing wave) gradients,they may not represent natural stimuli [20, 21]. Hence, cells may not have to solve thetraditional ‘back-of-the-wave’ problem, but instead have to decide which pulse to follow.However, this difficulty is eased as cells secrete cAMP from their rear [35].Our multi-scale model captures true emergence, generally not included in previousmodels of Dictyostelium aggregation. Models from the 1990s focused on the descriptionof cell populations and the generation of spiral waves [14, 26–28, 40]. These werefollowed by the biologically more detailed LEGI [24, 25] and Meinhardt [23] models toaddress the single-cell response to chemoattractant gradients. More recently, theFitzHugh-Nagumo model was adopted to explain the pulsing and synchronization ofmultiple cells (see
Supporting Information for a comparison) [31, 32]. Furthermore,hybrid models were proposed [41]. However, none of these models started from adetailed spatio-temporal single-cell model and was able to describe the type of orderand exact transition point for achieving collective behavior.When dealing with complex biological phenomena, there are necessarily limitationsin the deduced models and acquired data. To assess criticality via finite-size scaling,
PLOS
Dictyostelium , neurons in the brain, biofilms,embryos, tumors) [45–47] and animal groups (such as bird flocks) [8, 9, 48–50].Interestingly, neurons and bacteria pulse (spike) as well [51]. Operating at criticality, i.e.the tipping point between order and disorder, may allow cells to be maximallyresponsive, to communicate robustly over long distances, to act as a single coherentunit, and to make decisions on, e.g., when and where to aggregate. In the future, itwould be fascinating to conduct aggregation experiments in 3D environments, and tostudy the collective response to perturbations such as obstacles, changes in temperature,and exposure to toxins.
Materials and Methods
Detailed model
The intracellular cAMP dynamics are described by the FitzHugh-Nagumo model, aclassical model to reproduce neuronal spiking and previously adopted to describeexcitability in
Dictyostelium [31, 32]. Degradation of intracellular cAMP is achieved byphosphodiesterase RegA, which is negatively regulated by extracellular concentration ofcAMP (by means of ERK2 [34]). Secretion of cAMP from the cell rear [34, 35] is strictlycoupled to its intracellular concentration: if the extracellular cAMP concentration isbelow a threshold value cells exhibit a constant small leakage of cAMP, but a temporaryhigh concentration of cAMP is released during pulses of intracellular cAMP once abovethe threshold. If the extracellular cAMP concentration is kept above this threshold thecell becomes a sustained oscillator. Extracellular cAMP is degraded by thephosphodiesterase PDE [52]. This model correctly captures the relay of the signal andthe sustained pulsing observed in
Dictyostelium (see
Supporting Information for adetailed explanation).
Coarse-grained model
To reproduce the dynamics of thousands of cells, we simplified further therepresentation given by the detailed model. We assumed that cells are point-like objects,which secrete cAMP maximally at their rear. Specifically, spatial propagation of cAMPwas modeled as an exponential decay with a constant of 0.1 µm − (within a factor of 2of the value extracted from the detailed model simulations). The spatio-temporalconcentration profiles are rescaled according to the cosine of the angle with theopposite-to-motion direction; secretion becomes zero at 90 ◦ (lateral secretion) and is set PLOS c to determine if a given cell will emit a pulse or justleak cAMP, and a gradient threshold ∇ c determines if the cell will move randomly orfollow the local cAMP gradient. As for the detailed model, every cell undergoes arefractory period of 6 minutes after firing, during which it keeps the same motion it hadduring pulsing. To reproduce volume exclusion, cells cannot be closer to each other than3 µm (this rule is overwritten later in simulations, when cells are densely packed andlikely superimpose). To drastically speed up simulations, the algorithm is writtenwithout explicit modeling of diffusion of cAMP in space, instead it computes how muchcAMP every cell senses and what their spatial gradients are by considering positions ofcells with respect to each other. This implementation is able to reproduce aggregationof thousands of cells. More specifically, N =1000 cells were considered at experimentaldensity of about one monolayer (6600 cells/ µm ). For the other simulations of N =600,800 and 1200, the total area (of 389x389 µm ) was fixed and density varied accordingly.See Supporting Information for a detailed explanation.
Density pair correlation
The pair-correlation function was computed as described in [53], given by g ( r ) = AN ( N −
1) 12 πra N (cid:88) i (cid:54) = j δ ( r − r ij ) (3)where A is the total area of the image considered, N is the number of cells, r is theradius of a ring and a is the discretization constant. In case of a random distribution g ( r ) takes a value of 1 on average (similar to blue trace in Fig. 1Dii), while in case ofparticle clustering g ( r ) becomes greater for small distances (as for red trace in samepanel). Spatial information
All images were binarized (by means of MATLAB thresholding algorithms graythresh and im2bw for the case of experimental images). After that, 2D images were convertedin 3D binary matrices where the third dimension has a 1 corresponding to the pixelintensity (thus in this case, since the starting images were binary, the 3D matrix has a 1at level 0 if that pixel is black and at level 1 if it is white). This guaranteed that allimages had the same histogram, provided that they initially were of the same size. Forthe case of uncorrelated pixels, all Fourier coefficients P i are considered independentand Gaussian distributed. Image entropies were then calculated as: H kS = − N (cid:88) i P i log P i (4)where the probability density function P is Gaussian distributed with zero mean andvariance calculated from the sum of the pixel intensities. H kS is computed by dividingthe function into bins of width σ/
100 and summing P i log P i from − σ to 10 σ .Fourier transformation was then applied to the image. The real and imaginary part ofthe Fourier coefficients were then considered to compute I kS = (cid:88) i ( − log P Ri − log P Ii ) (5) PLOS P Ri and P Ii refers to the real and imaginary part of coefficient i . The sum wascalculated by considering bins of width σ/
100 around the values assumed by the Fouriercoefficients. k -space spatial information kSI was finally calculated as kSI = H kS − I kS .For a more comprehensive explanation, see [38]. Directional correlations and susceptibility
To calculate the connected correlations, local averages of the velocities are subtractedfrom cell velocities. For every cell we considered the average movement of all cells in itsneighborhood up to a certain maximal distance r c , and compute the correlationsbetween the cell in the center and all the cells belonging to its neighborhood. Werepeated this procedure for every cell in our image. In this way we are able to decreasethe “bulk” velocity component in the fluctuations, while keeping a continuous partitionof the image (which we would have lost in case of rigid partition of the image in smallersquares), and without preassigning the final position of the aggregation center. In orderto understand better the influence of this partitioning on the calculation of theconnected correlations, we repeated the same procedure for different radii. Specifically,if L is the image dimension, we set r c equal to L/ L/ L/ L/
8, and L/
10, with L/ L/ L/ L/ L/ L/
6, and L/
8, were considered, and L/ Experimental methods
Time-lapse movies were obtained similar to protocol in [19, 31]. Axenic
Dictyostelium cells expressing Epac1camps were starved for 4-5 hours, and then plated on hydrophobicagar for imaging. Sixteen fields of view from a microscope are combined (1.2 x1.2 mm ), resulting in the recording of thousands of cells in a wide field (invertedepifluorescence microscope (TE300, Nikon). To allow high-precision tracking ofindividual cells in a dense cell population, a different fluorescent marker (TRED) isexpressed and mixed with unmarked cells so a subpopulation of cells could be tracked(10% TRED cells). See Supporting Information for further details.
Segmentation and tracking
Images of TRED channels were segmented by using the MATLAB function imextendedmax , which outputs a binary image given by the computation of the localmaxima of the input image. The centroids positions were then computed from thismask by means of the regionprops function. The tracking of individual cells was done byconsidering the centroid positions for different times. For every time t the nearestneighbor centroid at time t + 1 was found, and the trajectory was accepted if thedistance between the two positions was smaller than the average cell size. Supporting Information
S1 Fig. Coarse-grained model: cell sensing and behavior.
PLOS
Single cell in a box to simulatestreaming, see Fig. 1B of the main text.
S2 Video. Single-cell model: wave response.
Single cell in a box to simulatethe response to an external wave of cAMP, see Fig. 1C of the main text.
S3 Video. Single-cell model: aggregation.
Four cells in a box to simulateaggregation, see Fig. 1D of the main text.
S4 Video. Coarse-grained simulations. N =1000 cells, as shown if Figs. 1F and2 of the main text. S5 Video. Experimental movie.
Dataset 3, shown in Fig. 3 of the main text.
S6 Video. Global perturbation during prestreaming for coarse-grainedsimulations. N =1000 cells, see Fig. 1C of the main text. S7 Video. Global perturbation during streaming for coarse-grainedsimulations. N =1000 cells, see Fig. 1C of the main text. S8 Video. Local perturbation during prestreaming for coarse-grainedsimulations. N =1000 cells, see Fig. 1C of the main text. S9 Video. Local perturbation during streaming for coarse-grainedsimulations. N =1000 cells, see Fig. 1C of the main text. S1 Appendix. Supporting information.
Including detailed explanation ofexperimental procedure, single-cell and coarse-grained model and experimental dataanalysis.
PLOS cknowledgments
We are grateful to the Gregor lab at Princeton University for sharing their data with us,and additionally Thomas Gregor, Allyson Sgro, and Monika Skoge for helpfuldiscussions. We also thank Luke Tweedy for help with the detailed model, LinusSchumacher for comments on the manuscript, and Mariam Elgabry and Suhail Islam forsupport with computational issues. GDP and RGE were supported by ERC StartingGrant 280492-PPHPI (http://erc.europa.eu/starting-grants).
References
1. David Sumpter, Jerome Buhl, Dora Biro, and Iain Couzin. Information transferin moving animal groups.
Theory in Biosciences , 127(2):177–186, 2008.2. Tam´as Vicsek and Anna Zafeiris. Collective motion.
Physics Reports , 517(3–4):71– 140, 2012. Collective motion.3. Elod Mehes and Tamas Vicsek. Collective motion of cells: from experiments tomodels.
Integrative Biology , 6:831–854, 2014.4. Hilary L. Ashe and James Briscoe. The interpretation of morphogen gradients.
Development , 133(3):385–394, 2006.5. Alexander H. Morrison, Martin Scheeler, Julien Dubuis, and Thomas Gregor.Quantifying the bicoid morphogen gradient in living fly embryos.
Cold SpringHarbor Protocols , 2012(4), 2012.6. Susanne C. van den Brink, Peter Baillie-Johnson, Tina Balayo, Anna-KaterinaHadjantonakis, Sonja Nowotschin, David A. Turner, and Alfonso Martinez Arias.Symmetry breaking, germ layer specification and axial organisation in aggregatesof mouse embryonic stem cells.
Development , 141(22):4231–4242, 2014.7. Sara Kaliman, Christina Jayachandran, Florian Rehfeldt, and Ana-SunˇcanaSmith. Novel growth regime of { MDCK } { II } model tissues on soft substrates. Biophysical Journal , 106(7):L25 – L28, 2014.8. Alessandro Attanasi, Andrea Cavagna, Lorenzo Del Castello, Irene Giardina,Stefania Melillo, Leonardo Parisi, Oliver Pohl, Bruno Rossaro, Edward Shen,Edmondo Silvestri, and Massimiliano Viale. Finite-size scaling as a way to probenear-criticality in natural swarms.
Physical Review Letters , 113:238102, Dec 2014.9. Andrea Cavagna, Alessio Cimarelli, Irene Giardina, Giorgio Parisi, RaffaeleSantagati, Fabio Stefanini, and Massimiliano Viale. Scale-free correlations instarling flocks.
Proceedings of the National Academy of Sciences ,107(26):11865–11870, 2010.10. Herbert Levine and Wouter-Jan Rappel. The physics of eukaryotic chemotaxis.
Physics Today , 66(2):24–30, 2013.11. Till Bretschneider, James Jonkman, Jana K¨ohler, Ohad Medalia, KarmelaBarisic, Igor Weber, Ernst H.K. Stelzer, Wolfgang Baumeister, and G¨untherGerisch. Dynamic organization of the actin system in the motile cells ofdictyostelium.
Journal of Muscle Research & Cell Motility , 23(7):639–649, 2002.12. P Devreotes. Dictyostelium discoideum: a model system for cell-cell interactionsin development.
Science , 245(4922):1054–1058, 1989.
PLOS
Developmental Biology , 204(2):525 – 536, 1998.14. Cornelis J Weijer. Morphogenetic cell movement in dictyostelium.
Seminars inCell & Developmental Biology , 10(6):609 – 619, 1999.15. Fernanda Alcantara and Marilyn Monk. Signal propagation during aggregation inthe slime mould dictyostelium discoideum.
Microbiology , 85(2):321–334, 1974.16. J.D. Gross, M.J. Peacey, and D.J. Trevan. Signal emission and signal propagationduring early aggregation in dictyostelium discoideum.
Journal of Cell Science ,22(3):645–656, 1976.17. KJ Tomchik and PN Devreotes. Adenosine 3’,5’-monophosphate waves indictyostelium discoideum: a demonstration by isotope dilution–fluorography.
Science , 212(4493):443–446, 1981.18. Jean-Louis Martiel and Albert Goldbeter. A model based on receptordesensitization for cyclic AMP signaling in dictyostelium cells.
BiophysicalJournal , 52(5):807 – 828, 1987.19. Thomas Gregor, Koichi Fujimoto, Noritaka Masaki, and Satoshi Sawai. The onsetof collective behavior in social amoebae.
Science , 328(5981):1021–1025, 2010.20. Monica Skoge, Haicen Yue, Michael Erickstad, Albert Bae, Herbert Levine, AlexGroisman, William F. Loomis, and Wouter-Jan Rappel. Cellular memory ineukaryotic chemotaxis.
Proceedings of the National Academy of Sciences ,111(40):14448–14453, 2014.21. Akihiko Nakajima, Shuji Ishihara, Daisuke Imoto, and Satoshi Sawai. Rectifieddirectional sensing in long-range cell migration.
Nature Communications , 5, 112014.22. Luke Tweedy, B¨orn Meier, J¨urgen Stephan, Doris Heinrich, and Robert G.Endres. Distinct cell shapes determine accurate chemotaxis.
Scientific Reports ,3:2606, 09 2013.23. H. Meinhardt. Orientation of chemotactic cells and growth cones: models andmechanisms.
Journal of Cell Science , 112(17):2867–2874, 1999.24. Carole A. Parent and Peter N. Devreotes. A cell’s sense of direction.
Science ,284(5415):765–770, 1999.25. Andre Levchenko and Pablo A. Iglesias. Models of eukaryotic gradient sensing:Application to chemotaxis of amoebae and neutrophils.
Biophysical Journal ,82(1):50 – 63, 2002.26. David A. Kessler and Herbert Levine. Pattern formation in
Dictyostelium via thedynamics of cooperative biological entities.
Physical Review E , 48:4801–4804, Dec1993.27. Herbert Levine. Modeling spatial patterns in dictyostelium.
Chaos , 4(3):563–568,1994.28. H Levine, I Aranson, L Tsimring, and T V Truong. Positive genetic feedbackgoverns camp spiral wave formation in dictyostelium.
Proceedings of the NationalAcademy of Sciences , 93(13):6382–6386, 1996.
PLOS
Interface Focus , 2(4):487–496, 2012.30. David M. Umulis and Hans G. Othmer. Mechanisms of scaling in patternformation.
Development , 140(24):4830–4843, 2013.31. Allyson E Sgro, David J Schwab, Javad Noorbakhsh, Troy Mestler, PankajMehta, and Thomas Gregor. From intracellular signaling to populationoscillations: bridging size- and time-scales in collective behavior.
MolecularSystems Biology , 11(1), 2015.32. Javad Noorbakhsh, David J. Schwab, Allyson E. Sgro, Thomas Gregor, andPankaj Mehta. Modeling oscillations and spiral waves in
Dictyostelium populations.
Physical Review E , 91:062711, Jun 2015.33. Matthew P. Neilson, Douwe M. Veltman, Peter J. M. van Haastert, Steven D.Webb, John A. Mackenzie, and Robert H. Insall. Chemotaxis: A feedback-basedcomputational model robustly predicts multiple aspects of real cell behaviour.
PLoS Biology , 9(5):e1000618, 05 2011.34. V.C. McMains, X. Liao, and A.R. Kimmel. Oscillatory signaling and networkresponses during the development of
Dictyostelium discoideum . Ageing ResearchReviews , 7(3):234–48, 2008.35. Paul W Kriebel, Valarie A Barr, Erin C Rericha, Guofeng Zhang, and Carole AParent. Collective cell migration requires vesicular trafficking for chemoattractantdelivery at the trailing edge.
The Journal of Cell Biology , 183(5):949–961, 122008.36. J.J. Tyson and J.D. Murray. Cyclic amp waves during aggregation ofdictyostelium amoebae.
Development , 106(3):421–426, 1989.37. Peter J. M. Van Haastert. A stochastic model for chemotaxis based on theordered extension of pseudopods.
Biophysical Journal , 99(10):3345–3354,2015/12/17 2010.38. WilliamF. Heinz, JeffreyL. Werbin, Eaton Lattman, and JanH. Hoh. Computingspatial information from fourier coefficient distributions.
The Journal ofMembrane Biology , 241(2):59–68, 2011.39. Luke Tweedy, David A. Knecht, Gillian M. Mackay, and Robert H. Insall.Self-generated chemoattractant gradients: Attractant depletion extends the rangeand robustness of chemotaxis.
PLoS Biology , 14(3):1–22, 03 2016.40. Eirikur Palsson and Hans G. Othmer. A model for individual and collective cellmovement in dictyostelium discoideum.
Proceedings of the National Academy ofSciences , 97(19):10448–10453, 2000.41. Arpan Bhowmik, Wouter-Jan Rappel, and Herbert Levine. Excitable waves anddirection-sensing in dictyostelium discoideum : steps towards a chemotaxis model.
Physical Biology , 13(1):016002, 2016.42. N. D. Mermin and H. Wagner. Absence of ferromagnetism or antiferromagnetismin one- or two-dimensional isotropic heisenberg models.
Physical Review Letters ,17:1133–1136, Nov 1966.
PLOS
Journal of Physics C: Solid State Physics , 6(7):1181,1973.44. John Toner and Yuhai Tu. Flocks, herds, and schools: A quantitative theory offlocking.
Physical Review E , 58:4828–4858, Oct 1998.45. Paolo Moretti and Miguel A. Mu˜noz. Griffiths phases and the stretching ofcriticality in brain networks.
Nature Communications , 4:2521, 10 2013.46. Dante R. Chialvo. Emergent complex neural dynamics.
Nature Physics ,6(10):744–750, 10 2010.47. Dmitry Krotov, Julien O. Dubuis, Thomas Gregor, and William Bialek.Morphogenesis at criticality.
Proceedings of the National Academy of Sciences ,111(10):3683–3688, 2014.48. Tam´as Vicsek, Andr´as Czir´ok, Eshel Ben-Jacob, Inon Cohen, and Ofer Shochet.Novel type of phase transition in a system of self-driven particles.
PhysicalReview Letters , 75:1226–1229, Aug 1995.49. Mate Nagy, Zsuzsa Akos, Dora Biro, and Tamas Vicsek. Hierarchical groupdynamics in pigeon flocks.
Nature , 464(7290):890–893, 04 2010.50. William Bialek, Andrea Cavagna, Irene Giardina, Thierry Mora, EdmondoSilvestri, Massimiliano Viale, and Aleksandra M. Walczak. Statistical mechanicsfor natural flocks of birds.
Proceedings of the National Academy of Sciences ,109(13):4786–4791, 2012.51. Arthur Prindle, Jintao Liu, Munehiro Asally, San Ly, Jordi Garcia-Ojalvo, andGurol M. Suel. Ion channels enable electrical communication in bacterialcommunities.
Nature , 527(7576):59–63, 11 2015.52. Sonya Bader, Arjan Kortholt, and Peter J. M. Van Haastert. Seven dictyosteliumdiscoideum phosphodiesterases degrade three pools of camp and cgmp.
Biochemical Journal , 402(1):153–161, 2007.53. Thomas Gurry, Ozan Kahramano˘gulları, and Robert G. Endres. Biophysicalmechanism for ras-nanocluster formation and signaling in plasma membrane.
PLoS ONE , 4(7):e6148, 07 2009.
PLOS lobalinhibitor Localinhibitor Localactivator Cell cAMP in RegA cAMP ex Cell x y00.40.8 * x y l r u d00.40.8 * Membrane dynamics Intracellular dynamics Extracellular dynamicsStreaming Wave response Aggregation
PDE cAMP ex Cell AB C I C I ρ (r) puls./secr.no cAMP * C D EF non-firing cell firing cell I n f o r m a t i on [ no r m ] Time [minutes] modelexp i i i iii ii ii ii
Fig 1. Multiscale model: from single cell-shape changes and chemotaxis tocollective behavior. (A) Schematics for membrane dynamics (left), intracellularcAMP dynamics (center), and extracellular cAMP dynamics (right). (B) Single-cell“streaming” simulation in a box with periodic boundary conditions and a constantconcentration of cAMP (i). Box dimensions are about 25x90 µm (the initial cell radiusis assumed ∼ µm ). Due to the small dimension of the box the cell is just leaking, notpulsing, in order to avoid saturation of secreted cAMP. The simulation was repeated 12times and the average chemotactic index (CI) calculated (ii). Errorbars representstandard errors. Differences in CI x and CI y are statistically significant (p < σ ∼ µm ) moves from right to left with a speed of about 300 µ m/min [36]. 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 directions, as well as for left, right, up and down directions in order todiscriminate between the directions of the incoming (right direction) and outgoing (leftdirection) wave. Simulations are repeated 12 times; shown are averages and standarderrors. Box is about 60x105 µm . CI in the right direction is significantly higher than CIin the other directions. (D) “Aggregation simulation”. (i) Four cells are simulatedmoving in a constant concentration of cAMP. At the beginning, cells are randomlydistributed. (ii) Density correlation at the end of simulations is plotted for control cellswithout secretion (blue) and all cells leaking cAMP and one cell also emitting pulses ofcAMP (red). The red line has a significant (p < µm . See Materials and Methods fordetails on density correlation and Supporting Information for a full explanation of thedetailed model. (E) Schematic showing cells represented as point-like objects withvelocity vectors. Firing cells emit pulses of cAMP, non-firing cells secrete cAMP at alow constant leakage rate. Spatial cAMP profiles are derived from detailed modelsimulations. At every time point cells are allowed two possible directions of movementin order to reproduce pseudopod formation at the cell front, with directions changing by ± . ◦ with respect to the previous movement, corresponding to an angle betweenpseudopods of about 55 ◦ [37]. (F) Screenshot during streaming for N =1000 simulatedcells (i). (ii) Spatial information versus time: simulations (blue) compared withexperimental dataset 3 (green). Values were then normalized and shifted in time tofacilitate comparison. PLOS
Time [min] 1000 0 64 200 D i s t an c e [ µ m ]
100 200 0 0 389 B Non-connected correlationConnected correlation
C D E F L/10L/6L/4 L/2102030 Radiusξ C o rr e l a t i on [ no r m ] Distance [ξ ]0 0.5 1 1.5 χ12001000800600N Fig 2. Coarse-grained model leads to collective behavior. (A) Screenshot ofsimulation for N =1000 cells at different time points: prestreaming (left), streaming(center) and after aggregation (right). Red (yellow) points represent non-firing (firing)cells. See Supporting Information for the full movie. (B) Kymograph of non-connected C nc and connected C c directional correlations for simulation of N =1000 cells.Directional correlations profiles C ( r ) were calculated for every time frame, anddisplayed depending on distance r . (C) Correlation length versus time for differentnumbers of cells. Data were smoothed with a moving average filter spanning 10consecutive frames. Inset: number of cells considered in simulations ( N =1200, 1000,800, 600). (D) Susceptibility χ plotted with respect to nearest-neighbor (NN) distancefor different number of cells and with respect to time (inset). The peak in susceptibilitybecomes the higher the larger the number of cells, and NN distances decreaseaccordingly. Profiles in the inset were smoothed with a moving average filter spanning10 points. (E) Comparison of correlation profiles for streaming phase (50-minute timewindow). Connected correlations were calculated for different numbers of cells and thennormalized and plotted against their respective correlation lengths. The four profilescollapse onto a single curve, independently of the number of cells. (F) Averagecorrelation length versus neighborhood radius for different number of simulated cells. L corresponds to the size of the images (389 µ m). ξ represents the average correlationlength during streaming phase (50 minutes window). Error bars represent standarderror. See Supporting Information for a full explanation of the model.
PLOS ig 3. Collective behavior in experimental data for validating model. (A)CFP images of
Dictyostelium aggregation of dataset 3. Images were taken after 4-5hours of starvation when cells were still moving randomly before initiating aggregation(left), during streaming phase (450 minutes after first image, center) and afteraggregation (800 min, right). (B) Kymograph of non-connected C nc and connected C c directional correlations for movie in dataset 3. Distance r is expressed in units ofaverage cell size (estimated after an ellipse was fitted to every cell contour, andcorresponding to the average of the minor axis, ∼ . µm ). (C) Spatial information(blue) and susceptibility χ (green) of movie in dataset 3 as a function of time. Theincrease in spatial information denoting a more ordered image corresponds to the peakin 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 forbetter visualization. Inset: comparison of cell number estimated from TRED imagesduring streaming phase for different movies. (E) Susceptibility χ as a function ofrescaled nearest-neighbor (NN) distance and as a function of time (inset). Note thatheight of peaks increases and that the corresponding rescaled NN distance decreaseswith number of cells, as for simulations. Rescaled NN distance was computed bynormalizing NN distance by the average cell size. In order to decrease noise, profiles inthe 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 everydataset was calculated as an average over 150 minutes of streaming phase. Error barsrepresent standard error. As in the simulated data curve collapse for different numbersof cells. (G) Average correlation length versus neighborhood radius. L corresponds tothe size of images (2033 pixels, ∼ ξ represents the average of 150 min duringstreaming phase. Error bars represent standard error. PLOS
B C
Time [min]1234 I n f o r m a t i on [ b i t s ] ×10
100 200 30002 [ b i t s / m i n ] x10 Time [min]1234 I n f o r m a t i on [ b i t s ] ×10 controlglobal psglobal slocal pslocal s Time [min] I n f o r m a t i on [ b i t s ]
250 500 7500510x10 [ b i t s / m i n ] Fig 4. Predictions from coarse-grained model. (A) Spatial information changeswith N in in silico data. Inset: Time derivative of spatial information profiles. Thechange in spatial information is larger for higher cell numbers. (B) Spatial informationas a function of time computed on the experimental data (as in panel (A), derivative isshown in the inset). As for the simulated data, the derivative tends to show higherpeaks for experiments with higher cell densities. (C) Effect of perturbations on thesystem during aggregation compared to control without perturbations. A speeding up ofaggregation is seen if a localized or a global, spatially uniform pulse of cAMP is given tothe system during prestreaming (ps). No effect on aggregation speed is noticed if thesystem is perturbed during streaming (s). See Supporting Information for the fullmovies (supporting movies S6-S9).