Emergence of synchronized multicellular mechanosensing from spatiotemporal integration of heterogeneous single-cell information transfer

We quantitatively characterize how noisy and heterogeneous behaviors of individual cells are integrated across a population toward multicellular synchronization by studying the calcium dynamics in mechanically stimulated monolayers of endothelial cells. We used information-theory to quantify the asymmetric information-transfer between pairs of cells and define quantitative measures of how single cells receive or transmit information in the multicellular network. We find that cells take different roles in intercellular information-transfer and that this heterogeneity is associated with synchronization. Cells tended to maintain their roles between consecutive cycles of mechanical stimuli and reinforced them over time, suggesting the existence of a cellular “memory” in intercellular information transfer. Interestingly, we identified a subpopulation of cells characterized by higher probability of both receiving and transmitting information. These “communication hub” roles were stable - once a cell switched to a “communication hub” role it was less probable to switch to other roles. This stableness property of the cells led to gradual enrichment of communication hubs that was associated with the establishment of synchronization. Our analysis demonstrated that multicellular synchronization was established by effective information spread from the (local) single cell to the (global) group scale in the multicellular network. Altogether, we suggest that multicellular synchronization is driven by single cell communication properties, including heterogeneity, functional memory and information flow.


Introduction
Synchronized multicellular dynamics is the basis of many critical physiological processes, such as the rhythmic beating of cardiomyocytes, planar cell polarity, regulation of vascular tone, and brain activities. However, a fundamental question remains elusive: how synchronization in the group emerges from the interactions of individual cells, each making stochastic decisions based on noisy cues from their local environment?
Every cell influence and respond to cells in its community through a complex interplay of chemical and physical signals 1 . Phenotypic heterogeneity, also referred to as extrinsic noise, occurs when cells respond differentially to the same external signal. Indeed, phenotypic heterogeneity is prevalent in multicellular systems, even when all the cells in the group share the same genetic background. Non-genetic heterogeneity, which may be due to stochasticity in gene expression or differences in gene expression, alternative splicing, and post translational modifications caused by external cues [2][3][4][5][6] , may have severe implications in disease 7,8 and developmental processes [9][10][11][12][13] . However, it is not known whether non-genetic heterogeneity plays a constructive role in information transfer between cells in multicellular systems. As such, we ask whether some cells within a group function as leaders or followers, promoting the spread of information through the group, while others act individually, and whether such heterogeneity is important for the synchronization of multicellular dynamics.
Previously, we demonstrated that cell-cell communication through gap junctions modulated ATP-induced calcium signaling in monolayers of fibroblast cells 14 . Tuning the levels of intercellular communications, by varying cell densities, inserting weakly-communicating cells, and by pharmacologically inhibiting gap junctions, controlled the temporal coordination of calcium signalling in neighboring cells 14,15 . To generalize these results and to elucidate how information transfer between single cells is integrated to synchronize population-level cellular responses, here we study collective mechanosensing in endothelial cells. Biologically, endothelial cells line the interior surface of blood vessels and form a monolayer that experiences varying levels of shear stress from blood flow 16,17 . Upon changing the flow rate (e.g., during acute wound), endothelial cells detect the change in shear stress, inform other cells such as smooth muscle cells, and adjust their internal signaling accordingly. Central to the cascade of events, shear stress leads to downstream ATP activation which modulates calcium signaling at the subcellular scale [18][19][20] . As a group, the endothelial cells must coordinate their signaling dynamics to achieve a coherent and collective response. Specifically, intercellular calcium levels are synchronized via gap junction-mediated cell-cell communication 14,21 . Such synchronized calcium signaling is instrumental in modulating reepithelialization, angiogenesis, and extracellular matrix remodeling, which are essential processes in wound repair [22][23][24][25][26] . Thus, we developed an integrated experimental-computational platform to quantitatively evaluate the roles that single cells take during the emergence of multicellular synchronization. Using this platform, we identified three key functions whereby single cells contribute to collective information processing that ultimately leads to multicellular synchronization.

Endothelial cells in a monolayer synchronize their calcium dynamics in response to external shear stress
We employed a microfluidics system that can precisely control the temporal profile of the shear stress that the cells experience (Fig. 1A top). We grew confluent monolayers of HUVEC (Human Umbilical Vascular Endothelial Cell) cells on the bottom surface of the flow channels (Fig. 1A bottom). A computer interfaced flow switch regulated input pressure to induce smooth flow profiles in the microfluidic channel as verified by particle image velocimetry (Fig. 1B). The shear stress -induced calcium signal of the HUVEC cells was imaged with the fluorescent calcium indicator Calbryte-520 with single cell resolution ( Fig. 1A inset, Video S1). We manually marked each cell center ( Fig. 1A inset), recorded the intracellular calcium signal as a time-series of fluorescent intensity for every cell and verified that the magnitude of the cell's calcium signal correlated with the magnitude of the corresponding flow shear stress that the cells experienced (Fig. 1C, Methods). This setting enabled us to investigate the collective calcium response of HUVEC cells when subjected to precisely controlled shear stresses.
Upon exposing the cells to a step-like increase in shear stress (0.1-0.2 Pa), which is similar to those that an endothelial cell experience during acute bleeding 27 , the variability in the cells' temporal derivative of their calcium signal (termed calcium dynamics, annotated ̂( ) (Methods)) increased and then gradually reduced until it converged to a synchronized steady state (Fig. 1D). We defined three measures to quantify multicellular synchronization that relied on the population-level standard deviation ̂( ) ( ) of single cell calcium dynamics ̂( ) at different time points. First, asynchronization is defined as the time average of ̂( ) ( ) over the final 200 seconds. Low values of asynchronization implied improved synchronization across the entire group (Fig. 1E, marked in cyan). Second, (̂( ) ( )) − (̂( ) ( )), measures the magnitude of the transition from an unsynchronized to a synchronized state, which we termed relative variability (Methods, Fig. 1E, marked in green). Third, we defined the synchronization rate (Methods, Fig. 1E, marked in red) as the ratio between the area under the curve of ̂( ) ( ) from the peak variability in calcium dynamics over 200 seconds (Fig. 1E, the area marked in yellow) with respect to the theoretical upper bound where the relative variability is zero (Fig. 1E, the combined areas marked in yellow and orange). We could not find a clear association between asynchronization and relative variability or synchronization rate, while the latter two measures were negatively correlated across experiments (Fig. 1F, Fig. S1). The mean time-correlation between neighboring cells, which we used previously 15 , was associated with the multicellular synchronization but not associated with the relative variability or the synchronization rate (Fig.   S2). These results demonstrate that multicellular synchronization is not necessarily correlated with the magnitude and speed of the response. We also found that increased levels of shear stress decreased multicellular asynchronization but not the relative variability or the synchronization rate (Fig. S3). Altogether, these results suggested that the cells formed a communication network that gradually evolved and reinforced synchronization despite the vast variability in single cell calcium response at the onset of shear stress application. In the following, we focus on relative low shear stress levels of 0.1-0.2 Pa to further elucidate the mechanisms of collective synchronization.  To avoid spurious cause-effect relations, Granger Causality requires the time-series being analyzed to be stationary, i.e., fluctuating signals with a consistent mean and variability.
Therefore, we excluded experiments where less than 85% of the cells passed two stationary tests termed hub/individual ratio (Fig. 2E), to the relative variability and to the synchronization rate, but not to the asynchronization ( Fig. S7A-C). These results suggest that increased fraction of communication hub cells along with decreased fraction of individuals were associated with an improved synchronization process (reduced relative variability and increased synchronization rate). Importantly, the hub/individual ratio is a local cell property -dependent on the spatial organization of cells in the vicinity of cells playing the hub or individual phenotype, as verified with a spatial permutation test (Fig. S8A). Thus, these results established a link between the "local" single cell communication properties and the process leading to the "global" emergent property of multicellular synchronization.
We next focused our attention to the fraction of cells taking the "communication hub" role. The assignment of a cell to the "hub" role was determined as 0.5 standard deviations above the mean in both the receiver and transmission score. Assuming that the transmission and receiver scores of a cell are independent (or negatively correlated), the expected fraction of "hubs" in the population would be bounded by P(zscore_transmission > 0.5) * P(zscore_receiver > 0.5).
However, the observed fraction of cells in the communication hub role was enriched beyond this expectation across experiments (Fig. 2F, Fig. S7D).
The hub/individual ratio is a measure for heterogeneity, relying heavily on the extreme high and low values in the transmission and receiver scores distributions (Fig. 2E, top right and bottom left corners). More hubs increase this ratio, and more individuals decrease it. But how do we quantify the heterogeneity at the bulk of the transmission/receiver score space? We defined a second measure for heterogeneity of all cells with |transmission score + receiver score| < 1, which consists of approximately 40% of the cells in the population (Fig. 2E, central region). We calculated the distribution of the absolute sum of the transmission and receiver scores for the cells in this bulk region and defined the bulk homogeneity as a scalar value measuring the deviation from the uniform distribution (Methods, Fig. 2E). For high heterogeneity experiments, the bulk homogeneity was low, following a uniform distribution, while in more homogeneous experiments, the bulk homogeneity was higher as the data distribution deviated from uniformity.
We used the Earth Mover's Distance (EMD) 31,32 to calculate dissimilarities between the observed and the uniform distributions. The EMD of 1-dimensional distributions is defined as the minimal 'cost' to transform one distribution into the other 33  Together, these results suggest that the heterogeneity in the roles that single cells took in the multicellular network, or "division of labor", was associated with the process leading to the emergence of collective calcium synchronization, where "communication hubs" were enriched.
As a note, none of the heterogeneity or the other synchronization measures was associated with the asynchronization, which was not even significantly associated with the mean GC edge probability (Fig. S10). Density Estimate plot visualization of the normalized transmission and receiver score (blue gradient contours). Partitioning of the z-score normalized transmission-receiver space to five regions (blue dashed lines), their corresponding annotations or "roles" (red text) and single cell assignment (yellow dots). Individual: transmission and receiver z-score < -0.5, Common: transmission z-score in the range of -0.5-0.5 and receiver z-score < 0.5 or receiver z-score in the range of -0.5-0.5 and transmission z-score < 0.5. Leader: transmission z-score > 0.5 and receiver z-score < 0.5. Follower: receiver z-score > 0.5 and transmission z-score < 0.5. Leader: transmission and receiver z-score > 0.5. (E) Depiction of the two measures for heterogeneity overlayed on the normalized transmission and receiver score plot. Hub/individual ratio -the ratio between the fraction of cells in the communication hub role (magenta) and the fraction of cells in the individual role (orange). Bulk homogeneity -the distance (earth movers distance, EMD) between the experimentally observed and the uniform distribution for all cells that hold |transmission score + receiver score| < 1 (green). (F) Communication hubs were enriched beyond their expected probability. p-value = 0.0076 via Wilcoxon rank-sign test on P(hub)-P(z-transmission)*P(zreceiver). For panels F and G, N = 14 biological replicates at 0.1 or 0.2 Pa (see S7D). (G) Pairwise associations between two heterogeneity measures (in gray: hub/individual ratio, bulk homogeneity) and three synchronization measures (in yellow: asynchronization, relative variability, synchronization rate). Edges color represents the level of association, as quantified by magnitude of correlation coefficients, between a pair of measures. Note that some of the measures are positive measures for synchronization properties (relative variability, synchronization rate) and heterogeneity (hub/individual ratio), while other properties are negative measures (asynchronization and bulk homogeneity) so negative correlations between measures from these two groups imply associations between synchronization and heterogeneity (for example, the negative correlation between bulk homogeneity and synchronization rate). N = 14 biological replicates per association. Full data in Fig

Experiments with periodic mechanical stimuli reveal that information flow is associated with multicellular calcium synchronization
After characterizing the communication networks exhibited by endothelial cell monolayers to shear stress, we asked if the network could be trained to adapt to changing external stimuli. To this end, we extended our assay to multiple rounds of repeated mechanical stimuli (Video S2).
By treating each round as an independent cycle, and comparing single cell responses across cycles, we could focus on the evolution of synchronization in the multicellular communication network (Fig. 3A, Methods). We found that the HUVEC monolayer gradually reinforced synchronization as observed by the gradual decrease in the standard deviation of the cells' calcium dynamics ̂( ) ( ) (Fig. 3B). This trend of improved synchronization coincided with a gradual increase of the cell's mean GC edge probability (Fig. 3C) and the mean receiver and transmission scores across the population along the progression in cycles over time ( Fig. 3D-E,  Kernel density estimate plot visualization of the normalized transmission and receiver score over the cycles (blue gradient contours). Left: partitioning of the z-score normalized transmission-receiver space to five regions (blue dashed lines), each cell (yellow dot) was assigned to a group or "role" (red text) according to the region they resided at.

Functional cell memory: cells maintain their roles in the communication network and reinforce them over time
We next asked to what extent the communication properties of cells were intrinsic cellular properties. To this end we correlated single cells' transmission and receiver scores across the repeated mechanical stimulus cycles while testing the null hypothesis that these scores were assigned randomly between consecutive cycles. We found that single cells' transmission and receiver scores were strongly correlated between consecutive stimulus cycles, that this correlation gradually increased as cells underwent additional stimulus cycles (Fig. 4A, Fig.   S12A), and could not be explained by the autonomous cells' response to the external mechanical stimuli (Fig. S13). Measuring single cell correlation between larger temporal gaps of 2-4 cycles did not show a dramatic diminishing pattern suggesting that the cellular memory is stable for time scales of at least 4-8 minutes (Fig. S12B). These results suggest that cells maintain and gradually reinforce memory regarding their role in the multicellular communication network at the timescales relevant for collective synchronization.
The combined effect of the increasing information flow and the functional cell memory was associated with a gradual increase in the fraction of followers, leaders and communication hubs and a decreased fraction of common and individual cells (Fig. 4B, Video S6), leading to an increased ratio between "hub" and "individual" cells that was correlated with the increased synchronization (Fig. 4C). To follow the dynamic trajectory of single cells between roles we analyzed cell state plasticity, the probability of transitioning from one state (functional role) to another cell state in consecutive cycles. In particular, we computed the enrichment factor -transition probabilities between any two states and normalized the quantity by the fully random transition probabilities (Methods, Fig. S14). As expected from our earlier observation of functional memory (Fig. 4A), we found that cells tended to maintain their roles or "similar" roles, as reflected by self-transition enrichment factors above one (Fig. S14, Fig. 4D). Generally, single cells followed a temporal trajectory from the roles characterized with less communication capacity to roles with increased communication (Fig. 4D -showing edges only for enrichment factors > 1). Intriguingly, we found symmetric transition folds between the follower-and leader roles, and the transition from communication hub to the follower/leader roles was enriched compared to the opposite transition to a communication hub (Fig. 4D, Fig. S14). The communication hub role was found to be much more stable than other roles or transitions (2.4 fold dwell probability compared with a fully random process), explaining their rapid spread in the population (Fig. 4B). The fraction of communication hubs also systematically exceeded its expected fraction if transmission and receiver scores were independent. Along with the gradually increasing information flow, the results highlight the role of functional memory in communication hub enriching (Fig. 4E).  For this analysis, we randomly selected for each cell at each topological distance at most ten neighboring cells due to computational cost and performed FDR multiple hypotheses correction (see Methods for full details).

Discussion
The emergence of robust multicellular behaviors from heterogeneous single cell dynamics is a poorly understood, but fundamentally important phenomenon in living systems 34  We showed the cells were actively communicating with one another locally, as supported by multiple lines of evidence throughout our study. First, we demonstrated that both hub/individual ratio and bulk homogeneity depended on the spatial organization of cells in their vicinity (Fig.   S8). Second, we found that the activation time, a cell's autonomous response to the external stress, was not associated with the transmission or receiver score (Fig. S13A), which would also be conflicting with the enrichment of communication hubs (which are both leaders and followers). Third, cells "remembered" and reinforced their roles in the multicellular communication network over time, as a local, spatially-dependent property (Fig. 4A), but did not "remember" their activation time in previous cycles (Fig. S13B). Fourth, neighbor pair cross correlation was a local cell property throughout the experiment (Fig. 5A). Together, our data established the decoupling of the local cell-cell communication from the global external stimuli, and excluded the possibility that the emergence of multicellular synchronization was an artifact caused by differential cell autonomous response to the external shear stress that was applied on the group.
Our data reveal that single cells take different roles in cell-cell communication (heterogeneity or "division of labor"), gradually reinforce these roles ("functional cell memory"), and increase the connectivity ("information flow") in the multicellular network. These three mechanisms work in concert to synergistically contribute to the emergence of multicellular synchronization. We hypothesize that cell heterogeneity expands the dynamic range of multicellular responses that cell memory stabilizes the dynamics against intrinsic and extrinsic noise and that information flow sustains and reinforces the multicellular dynamics.
We found that heterogeneity in cells' communication properties were associated with improved convergence to synchronization (Fig. 2G). We also observed that the fractions of cells at each functional role, excluding individuals, became more balanced through periodic cycles ( Fig. 3E and 4B), in agreement with our hypothesis that heterogeneity constructively contributes to the synchronization of a noisy multicellular system. Cell heterogeneity could arise from stochastic Previous studies have reported multiple sources of microenvironment-dependent cell memory.
For instance, cells can remember past mechanical properties of their substrate, which influence their differentiation 46 . Cells can also sense changes in their extracellular signal by remembering past extracellular stimulation via a receptor-mediated mechanism 47 . In the context of collective cell migration, a recent study showed that cells remembered their polarized state independently of cell-cell junctions 48 , and another study revealed associative memory of electric field and chemoattractant at stimuli in a unicellular organism migration patterns 49 . In our study, single cell memory of communication properties contributes to the temporal evolution of the multicellular network to its synchronized state. The dissociation between a cell's activation time and its functional role in information processing underpins the dynamic nature of memory, which is also consistent with the unidirectional evolution of the multicellular network (Figs. 4-5).
Our study reveals a self-organized multicellular network that supports information flow from local to global scales. Such information may be carried by two main signaling mechanisms, juxtacrine (contact-dependent) and autocrine (secreted-dependent) 50 . A juxtacrine channel allows a cell to establish conversation with its (physically touching) immediate neighbors without interference from extracellular space. For HUVEC cells such communication can be realized by gap junctions 51 . On the other hand, an autocrine channel allows a cell to broadcast its information through diffusive messengers in the extracellular space. For HUVEC cells stresstriggered ATP release and ATP-induced calcium dynamics constitute an autocrine pathway 52 .
While both mechanisms could contribute to the information flow within the multicellular network, we suggest contact-dependent signaling as the dominant mechanism. While a recent study suggested that positive feedback of a diffusive signaling mechanism can drive accelerated, long-range information transmission 53 , the external flow in our system is likely to rapidly dilute the diffusive messenger 54 . The contact-dependent information flow hypothesis is also supported by our previous studies where we demonstrated that blocking gap junctions, or inserting weakly communicating cells impared the information flow 14,15 .
Altogether, our results suggest the following phenomenological model for multicellular synchronization. Cells are gradually "learning" the local network structure around them

Microfluidics
The organic elastomer polydimethylsiloxane (PDMS, Sylgard 184, Dow-Corning) used to create the microfluidic devices was comprised of a two part mixture -a base and curing agent -that were mixed in a 10:1 ratio, degassed, and poured over a stainless steel mold before curing at 65•C overnight. Once cured, the microfluidic devices were cut from the mold, inlet/outlet holes were punched, and the device was affixed to a No. 1.5 coverslip via corona treatment.
The cross section of the flow chambers was rectangular (2 mm X 1 mm). See Fig. 1A for depiction.

Applying controlled shear stress on the cells
The microfluidics flow rate was controlled by a PID-regulated pressure pump and was monitored using an inline flow sensor (Elveflow). To verify the stability of the flow profile we mixed 1 micrometer fluorescent particles in the solution and used particle image velocimetry to quantify the flow rate (Fig. 1B). To calculate the shear stress, we approximated the flow profile in the flow chamber as low-reynolds number pipe flow. We considered the cells in the field of view to experience uniform shear stress calculated at the center of the flow chamber. This was possible because the imaging window was narrow (470 μm) compared with the chamber width (1 mm).
In the "step" experiments we exposed the cells to a "step"-like shear stress of 0.

Measuring single cell calcium signaling
We manually annotated every cell center (Fig. 1A inset),

Measuring synchronization
We defined three measures to quantify multicellular synchronization that relied on the standard deviation of single cell calcium dynamics at different time points (̂( )).
The first measure was asynchronization, which measured the mean standard deviation of the cells' calcium dynamic at the cells' synchronized state (Fig. 1E,  The second measure was relative variability (Fig. 1E, marked in green), which measured the magnitude of the transition from an unsynchronized to a synchronized state by calculating max(̂( ) ( ))-min(̂( ) ( )).
The third measure was synchronization rate (Fig. 1E, marked in red), which measured the ratio between the area under the curve of (̂( ) ( )) from the peak variability in calcium dynamics over 200 seconds (Fig. 1E, the area marked in yellow) with respect to the theoretical upper bound where the relative variability is zero (Fig. 1E, the combined areas marked in yellow and orange). The area under the curve was calculated using trapezoidal rule where the space between sample points is 2 sec (1 frame), ∫̂( ) ( ) .

Granger causality
Granger causality (GC) is a statistical method to quantify the information flow among multiple variables' time-series 28 . Intuitively, time-series B is said to be "Granger causal" of time-series A, if the variability of A can be better explained by previous values of B and A, compared to using only previous values of A. Granger causality is an approximation to "transfer entropy" and under the assumption of Gaussian distribution it is exactly equivalent 56 .
Formally, given two-time series ( ) and ( ), where t∈Z. The autoregressive model of is: Where, p is the lag order, the number of previous observations used for prediction, is the coefficient of and is the prediction error at time t. The autoregressive model of based also on the previous observations of is: Where, p is the lag order, is the coefficient of ( − ), is the coefficient of ( − ) and ( ) is the joint error of and predicting .

Stationarity test
To avoid spurious causality connection, and both must comply with a stationary process before applying the granger causality test. Intuitively, stationary means that the statistical

Granger causality statistical test
We applied a statistical test to infer granger causality between two ̂( ) time series and (denoted, −> ) 28 . GC tests the null hypothesis that is not contributing to the explained variance of (Equation (2)) in relation to the model derived solely from past values of (Equation (1)). This null hypothesis is rejected when at least one of the coefficients in equation (2) is different from zero. The statistic is based on the distribution: Where, is the sum of square residuals of the model which take into account only self-previous observation of the random variable (Equation (1)), and is the sum of square residuals of the other model which also takes into account the previous observation of the second random variable (Equation (2)). T is the sample size (number of observations in the time series used for prediction), p is the number of variables which was removed from the unrestricted model, in our case, the lag order, and k is the number of variables, in our case, twice the lag order. The null hypothesis is rejected when the Fvalue is larger than the F statistic (i.e., F's critical value) to conclude that −> . We derived the p-value from the F-statistic instead of directly using the F statistic to set up a global acceptance threshold.

Calculating the transmission and receiver scores
The transmission and the receiver scores were calculated as the probability for an outgoing (respectively, ingoing) edge at topological distance ≤ 2 (nearest and next-to-nearest neighbor), where the topological distance was calculated using the Delaunay triangulation. These neighborhood sizes were determined for sufficient observations for statistics, and the expected short-range communication between the cells. The transmission (denoted Tr) and receiver (denoted Re) scores were calculated independently for each cell using its cell neighbours at topological distance one and two Importantly, we treat each cell as an independent observation thus characterizing the role of each cell in the multicellular network without specifically committing on the exact network edges. To fix spurious edges due to multiple hypothesis testing we applied the strict Bonferroni correction that defines the edge significance threshold based on the number of edges considered 63 . In our case, with a significance threshold of 0.05 and n -number of potential edges we get a new significance threshold of 0.05/n. Edges passing these strict statistical tests are termed GC edges.

Partitioning the normalized transmission-receiver space
The transmission and receiver scores of each cell were normalized across the population to allow direct comparison of single cell heterogeneity. We calculated the receiver and transmission zscore for each cell , the variation from the mean in units of standard deviations: Tr_norm(ci) = (Tr(ci)-μ)/σ, where μ is the mean transmission score across the population, and σ is the standard deviation. The same normalization was applied for the receiver score. Kernel Density Estimation 64 was used for the visualization of the 2-dimensional normalized transmission and receiver score space (Fig. 2D, Fig. 3E). We partitioned the normalized transmission-receiver space to five regions, and assigned each cell to one of these regions. Individual: transmission and receiver zscore < -0.5. Common: transmission z-score in the range of -0.5-0.5 and receiver z-score < 0.5 or receiver z-score in the range of -0.5-0.5 and transmission z-score < 0.5. Leader: transmission zscore > 0.5 and receiver z-score < 0.5. Follower: receiver z-score > 0.5 and transmission z-score < 0.5. Leader: transmission and receiver z-score > 0.5. The z-score threshold of 0.5 was selected to maintain sufficient number of cells in each role for statistical analysis.

Measuring heterogeneity
We defined two measures for heterogeneity and applied it to each biological replica of the "step" experiments. First, the ratio between the number of hub cells and the number of individual cells (termed hub/individual ratio). Second, the bulk homogeneity, which was calculated as the Earth Mover's Distance (EMD) [31][32][33] between the experimentally observed and the uniform distribution for all cells that hold |transmission score + receiver score| < 1 (~40% of the cells), thus considering the combined GC outdegree and indegree scores for each cell. The bulk homogeneity was calculated based on observed and uniform distributions that were binned to 20 intervals to cover the [-1,1] range. Importantly, hub/individual ratio and bulk homogeneity were independent and complementary measures for heterogeneity because the cells considered for the bulk homogeneity calculation did not include hubs or individuals.
To verify that these heterogeneity measures are local, i.e., dependent on the spatial organization of the cells, we performed permutation tests. For hub/individual ratio we randomly switched the location of each hub and individual cell with another cell in topological distance > 2. This was followed by re-calculating the transmission and receiver scores for all the cells that were involved in this switching defining a new hub/individual ratio value. This process was repeated multiple times (Fig. S8A). For bulk homogeneity, we shuffled the cells' spatial locations and recalculated the bulk homogeneity based on the new spatial ordering. This was followed by correlating the bulk homogeneity with the synchronization measures, which were independent of the cells' spatial ordering, across all "step" experiments. This process was repeated ten times and compared to the series of bulk homogeneity scores of the step experiments with the original location of the cells (Fig. S8D-E).

Measuring information flow
GC edge probability was defined as the probability of a GC edge in the experiment. This was calculated as the ratio between the total number of GC edges and the total number of potential edges in the experiment. Because GC edge probability is a proxy for the information flow within the multicellular network, we also used the term information flow to refer to the GC edge probability in the manuscript text.

Earth Mover Distance (Wasserstein distance) between the cycles
The Earth Mover's Distance (EMD) [31][32][33] calculates dissimilarities between distributions as the minimum 'energy' required to shift one distribution to the other distribution. We used EMD to calculate the rate of change in the information flow throughout the cycles. To compare two cycles we calculated the EMD between the corresponding two-dimensional transmissionreceiver scores cell distributions. We discretized the distribution to 20x20 bins in the range of [-4, 4] in the transmission and receiver axes.

Enrichment factor of cellular state transitions
We calculated the enrichment factor, the fold change in the observed transition probabilities of To calculate the enrichment factor, the fold changes in switching from one state to another compared to the expected probability from a null model assuming random transitions drawn from the marginal state distribution, we calculated the expected transition matrix , (Fig. S14 bottom-left) and divided the Markov transition matrix bins the corresponding bins in the expected matrix (Fig. S14 right).

Measuring cell memory
To measure the cell memory we calculated the Pearson correlation of the cells transmission or receiver scores between consecutive cycles with step Δt (Δt = 1 in Fig. 4A, Δt ≥ 1 in Fig. S12B).
We evaluated the significance of our results using a permutation test by shuffling the cells' spatial locations with over 1000 permutations. The permutation test was performed by concatenating the vector scores of the cycles c, c+Δt, shuffling the values, splitting back to two vectors, and calculating the absolute Pearson correlation. The p-value is set as the fraction of permutations where the shuffled correlation surpassed the observed experimental correlation.

Activation Time
The activation time of a cell in a given cycle is the time where its calcium dynamics exceeds a threshold value of within a cycle. The threshold is parametrized by in the range of 0.1, 0.2 or 0.3 from the calcium dynamics range -the initial value subtracted from the maximal value within the cycle.
. The initial time was shifted by 60 seconds (30 frames) from the onset of the cycle to the time where the mean value of the cells' calcium dynamics is zero to ensure that the single cell calcium signal is on the rise for the vast majority of the cells.

Correlating the topological distance between pairs of cells to their GC-edge probability
In Fig. 5B we correlated the topological distance to the corresponding GC edge probabilities. For each topological distance, for each cell, we randomly selected ten (or less in topological distances with smaller numbers) cells and calculated the GC statistical test for each cell pair in both directions. We evaluated the critical value (i.e., p-value correction) using FDR, to correct for multiple hypothesis testing. Finally, we calculated the probability for GC significant edge as the total number of significant edges divided by total GC tests performed.

Software and data availability
The source code and sample data are publicly available in the following link:             The expected transition probability between two roles based on the cells' roles frequency over the cycles. Right: The enrichment transition matrix was followed by normalizing (scalar division) the Markov transition matrix with the expected transition matrix. The enrichment factor from hub to hub is higher than any other enrichment factor. (see only transitions above the expected, Fig. 4D).