Inference of long-range cell-cell mechanical communication from ECM remodeling fluctuations

Cells sense, manipulate and respond to their mechanical microenvironment in a plethora of physiological processes, yet the understanding of how cells transmit, receive and interpret environmental cues to communicate with distant cells is severely limited due to lack of tools to quantitatively infer the complex tangle of dynamic cell-cell interactions in complicated environments. We present a computational method to systematically infer and quantify longrange mechanical cell-cell communication through the extracellular matrix (cell-ECM-cell communication) by correlating ECM remodeling fluctuations in between communicating cells and demonstrating that these fluctuations contain sufficient information to define unique signatures that robustly distinguish between different pairs of communicating cells. We demonstrate our method with finite element simulations and live 3D imaging of fibroblasts and cancer cells embedded in fibrin gels. While previous studies relayed on the formation of a visible fibrous ‘band’ extending between cells to inform on mechanical communication, our method detected mechanical propagation even in cases where visible bands never formed. We revealed that while contractility is required, band formation is not necessary, for cell-ECM-cell communication, and that mechanical signals propagate from one cell to another even upon massive reduction in their contractility. Our method sets the stage to measure the fundamental aspects of intercellular long-range mechanical communication in physiological contexts and may provide a new functional readout for high content 3D image-based screening. The ability to infer cell-ECM-cell communication using standard confocal microscopy holds the promise for wide use and democratizing the method.

These deformations can form a visible fibrous band of aligned and dense fibers coupling the communicating cells mechanically, influencing the cells' internal molecular state 15 and active response 16 . This long-range mechanical intercellular communication through the ECM (termed here cell-ECM-cell communication) was observed for cells embedded in synthetic 2, 7 and biopolymer fibrous matrices 1, 3-6, 8, 9, 11-13 and are thought to play a key role in tissue regeneration and disease progression 5,[16][17][18][19][20][21][22][23] . In-vivo, fiber alignment bands can serve as ECM 'tracks' for cell migration with potential roles in wound healing, cancer metastasis and fibrosis 24,25 .
Quantification of mechanical cell-ECM-cell communication has been typically achieved by measuring the density, alignment or displacement of the remodeled ECM between the cell pairs 4, 5, 9-11, 13, 26-29 . While current measurements characterize the formation of the visible fibrous band extending between communicating cells, they lack the sensitivity to measure the dynamic reciprocal mechanical information transfer between the cells. Thus, current methods are limited in measuring potential cell-ECM-cell communication in the absence of visible bands, hampering our ability to distinguish which cells are actually communicating from the many cells that have the potential to communicate. Bridging this gap will enable tackling long-standing open questions in how tissues develop and diseases progress by enabling us to identify which cells are communicating with each other, and to what extent, in complex environments.
Here, we present a new computational method to quantify the transmitted ECM signal in between communicating cells by correlating temporal fluctuations of the remodeled matrix.
Computational simulations and 3D live imaging of fibroblasts and cancer cells embedded in fibrin gels demonstrate the power of our method in identifying unique ECM remodeling signatures that allow to robustly distinguish between different pairs of communicating cells.
Using this method, we were able to measure communication between cell pairs that do not form a visible 'band' of densified ECM, and after partial depletion of contractility upon Myosin-II inhibition. These results imply that mechanical signals propagate from one cell to another even in low cell contractility levels.

Analysis of ECM density between pairs of communicating cells using finite element simulations
During cell-cell mechanical communication, every localized fibrous region in between the communicating cells is affected by two components: the contractile force of the adjacent cell (referred as 'cell-ECM' interaction), and the force transmitted from the second distant cell (referred as 'cell-ECM-cell' communication) (Fig. 1A). To quantitatively characterize the independent contribution of each of these two components, we simulated contracting cells embedded within 2D fibrous networks using finite element discrete modeling, quantified the change in the ECM density between pairs of communicating cells and compared it to single cells interacting with the ECM. While these simulations do not reflect the true complexity of biological systems, they capture the essence of the mechanical elements of cell contractility and force propagation in fibrous nonlinear elastic networks (Methods). Thus, these simulations serve as a minimal mechanical model that enables us to independently tune mechanical parameters (e.g., cell-cell distance and cell contractility), and quantitatively infer their effect on ECM remodeling. These simulations will be used to ask whether ECM remodeling can robustly encode unique ECM signature in communicating cell pairs, forming the basis for our method.
To enable quantitative comparison between simulations, we normalized the fiber density to its zscore -the number of standard deviations away from the mean background fiber density at regions that were not influenced by the cells (Methods). In simulations of single contractile cells, we found that the regions next to the cell's boundary became denser following the application of 50% cell contraction (Fig. 1B). Upon the presence of a second cell, the overall densification was governed by the integrated contractile activity of both cells (Fig. 1A). This dual contribution led to the formation of a band of increased density along the connecting axis between the cells (Fig.

Distinguishing between communicating versus non-communicating simulated cell pairs using temporal correlation of local ECM remodeling fluctuations
Given that local ECM regions located along a band experience forces exerted by both cells and the discrepancy between ECM remodeling by a single and a pair of cells (Fig. 1), we hypothesized that local ECM remodeling fluctuations contain sufficient information to distinguish communicating cells from cells that are not communicating with one another. To test this hypothesis, we simulated the temporal process that led to the final remodeling. We consecutively applied 1% cell contraction for 50 steps, reaching 50% cell contraction. For each step, we recorded the fiber density between communicating cells. These iterative simulations captured the temporal dynamics of force propagation between the cells during cell-ECM-cell communication. As expected, the fiber density close to the cells' edge gradually increased over time ( Fig. 2A-B, Video S1).
We expected that the temporal correlation of ECM remodeling fluctuations between communicating cell pairs will exceed that of non-communicating cell pairs (Fig. 2C). In simulations of pairs of both communicating and non-communicating cells, we found positive correlations throughout the simulation in the fiber density dynamics and in its temporal derivative, i.e., change in fiber density over time (Fig. S1A). The correlation in fiber density is attributed to the monotonic increase in ECM densification between communicating cells (Fig.   2B). The correlation in the temporal derivative of fiber density is attributed to an association between fiber density and the change in fiber density (Fig. S1B): the ECM became denser over time, leading to increased temporal derivatives. However, these general trends lead to high correlations between any two cells regardless of whether they communicate with each other, and confound the unique fluctuations that may be attributed to a specific pair of communicating cells.
Thus, we aimed at overcoming this confounder by dampening correlations between noncommunicating cells toward zero. This was achieved by detrending, i.e., removing the temporal trends via a second temporal derivative. Although the correlations between non-communicating cell pairs were lost with a second temporal derivative, the correlations between communicating pairs were also lost and we were unable to distinguish between communicating and noncommunicating cell pairs (Fig. S1A).
We suspected that the inability to distinguish between communicating and non-communicating simulated cell pairs was due to lack of variability in ECM remodeling over time. Every simulated cell contracted in the exact same magnitude as any other simulated cell leading to nearly exact time-dependent remodeling of the adjacent ECM regime and thus masking the unique ECM fluctuation patterns that may exist between communicating partners. Thus, our hypothesis was that variability in the cells' contraction activity (an inherent property of cells) will lead to variability in the local ECM remodeling around each cell, leading to a unique communication signature between partners. In other words, we thought that heterogeneity in ECM remodeling would enable us to decouple the general ECM remodeling that occurs between any random pair of cells from the specific signal that is transmitted between communicating partners. Indeed, when introducing heterogeneity in simulated cell contraction, we distinguished communicating from non-communicating simulated cell pairs while maintaining the non-communicating correlations around zero (Fig. 2D). Cell heterogeneity was included in the simulations by drawing the instantaneous contraction rate of each cell independently from a normal distribution with mean ( ) contraction of 1% and a standard deviation ( ) in the range of 0-0.75% (std.) (Methods). Even a minimal standard deviation of 0.25% in cell contraction was sufficient to Quantifying fiber densification as a function of simulated contraction steps in between cell pairs. N = 20 cell pairs. Color code (left) and y-axis (right) encodes the ECM density in z-score, the number of standard deviation above background levels. The pair distance between simulated cell centers was set to 7 cell diameters. (C) Schematic sketch of comparing ECM remodeling fluctuations of communicating versus non-communicating cell pairs. The correlation between quantification windows of communicating pairs (upper in green) is evaluated in relation to the correlation between two cells from two different communicating pairs (lower in orange and purple). (D) Distinguishing between communicating and noncommunicating simulated cell pairs correlations of the fluctuations of the second derivative of fiber density with different levels of contraction heterogeneity. Each data point records the correlation between a simulated communicating cell pair, each box plots record the corresponding distribution of correlations between (the many) pairs of non-communicating cells. Mean contraction of 1% and standard deviation (heterogeneity) of 0%, 0.25%, 0.5%, and 0.75%. Three statistical tests are performed according to the following order. (1) Wilcoxon sum-rank testing the null hypothesis that the two correlation distributions of communicating and non-communicating cell pairs are originating from the same distribution. Wilcoxon sign-rank testing the null hypothesis that the correlation distributions of (2) communicating or (3) non-communicating are distributed around a mean 0. Std. = 0%: N = 20, p-values: (1) not significant, (2) not significant, (3) not significant. Std. = 0.25%: N = 19, p-values: (1) < 0.0001, (2) < 0.001, (3) not significant. Std. = 0.5%: N = 19, p-values: (1) < 0.0001, (2) < 0.001, (3) not significant. Std. = 0.75%: N = 20, p-values: (1) < 0.0001, (2) < 0.001, (3) not significant. (E) Correlation of communicating versus non-communicating cell pairs using the second derivative of fiber density dynamics as a function of the pair distance. Pair distance: 4 (N = 20 pairs, p-value < 0.0001), 5 (N = 19 pairs, p-value < 0.0001), 7 (N = 19 pairs, p-value < 0.01) and 9 (N = 19 pairs, p-value < 0.01).

Analysis of 3D ECM density between pairs of communicating fibroblast cells
To test the simulated results of correlation in ECM remodeling as a marker of mechanical cell-ECM-cell communication, we used 3D time-lapse confocal microscopy to image live NIH/3T3 fibroblast cells (GFP-Actin) embedded in 3D fluorescently labeled fibrin gels, where fibrin intensity was used as a proxy of fiber density (Methods). To visualize and measure fiber intensity in between cell pairs in 3D, we transformed the microscopy axes to a new 3D coordinate system that is aligned around the connecting axis between the cells' centers ( Fig. 3A, Fig. S2, Video S2, Methods) and performed z-score fiber intensity normalization in respect to the local background (Methods). While single cells did not show visually apparent fiber densification beyond regions close to the cell, approximately 60% of all cell pairs formed a visible band of increased density extending along the connecting axis between the cells (Fig. 3B) that showed increased mechanical coupling for cells pairs that were closer to one another (Fig. 3C, Fig. S3). The fiber density close to the cells' edge gradually increased over time in pairs of communicating cells but remained constant in single cells (Fig. 3D, Video S3 versus Video S4). The higher fiber intensity at the onset of imaging was attributed to the time (approximately 30 minutes) that passed from setting up the experiment until the onset of imaging (Fig. 3D, z-score of approximately 3 standard deviations above the background intensity). During this time, cells have already deformed the fibers around them. These results conclude that the formation of dense fibrous bands between fibroblasts cell pairs embedded in 3D fibrin gels are indicative of a mechanical coupling in concurrence with previous studies 5, [9][10][11] . Schematic sketch of the same-versus-different pair analysis. The correlation between quantification windows of communicating pairs ("same", green) is evaluated in relation to the correlation between one cell from that pair and another cell from a different communicating pair ("different", purple). (B) Sameversus-different pair analysis in simulations. Each data point records the correlation between the "same" and "different" cell pairs using the second derivative of fiber density dynamics. Mean contraction of 1% and standard deviation of 0.5%. All combinations of "same"/"different" were considered. Left: Pair distances: 4 (N = 20), 5 (N = 19), 7 (N = 19) and 9 (N = 19) cell diameters. "Same" pair had a higher correlation than "different" pair in 78% of the matched correlations. Right: "Same" -"different" correlation distributions for different pair distances: 4 (N = 20 pairs), 5 (N = 19 pairs), 7 (N = 19 pairs), and 9 (N = 19 pairs) cell diameters. 92% (pair distance = 4), 79% (pair distance = 5), 70% (pair distance = 7) and 68% (pair distance = 9) of data points were positive, implying that the "same" correlation is higher than the "different" correlation. Wilcoxon signed-rank testing the null hypothesis that the data are distributed around a mean 0: p-value < 0.0001 for all pair distances. (C) Schematic sketch of a cell, the quantification window and the normalization border. The normalization border is composed of two parts, above and below the quantification window, each with a width of approximately a cell diameter (15 μm), height of 0.5 cell size and with a gap of 0.25 from the quantification window (full information in Methods). (D) Same-versus-different pair analysis in experiments. Quantification windows were placed ~0.5 cell diameter (7.5 μm) above the connecting axis between cell pairs with a visible band, no offset in X/Y-axis. Correlations were calculated using the first derivative of fiber density dynamics. Top: Representative pair of communication cells shown after 255 minutes of imaging. Scale bar = 15 μm. Line connects the cell centers in the XY/Z space. For the displayed pair the Pearson correlation coefficient was 0.66 (Z offset = 0) and 0.71 (Z offset = 0.5). The same image is also shown in Fig. 3B. Bottom left: Pair distance = 60-150 μm (~4-10 cell diameters). N = 48 cell pairs, all combinations of "same"/"different" were considered. "Same" pair had a higher correlation than "different" pair in 94% of the matched correlations. Bottom right: "Same" -"different" correlation distributions for different pair distances: 4-6 (N = 18 pairs), 6-8 (N = 19 pairs), and 8-10 (N = 11 pairs) cell diameters. "Same" pair had a higher correlation than "different" pair in 93% (pair distance 4-6), 96% (pair distance 6-8) and 91% (pair distance 8-10). Wilcoxon signed-rank test p-value < 0.0001 for all pair distances. (E) Same-versusdifferent pair analysis for cell triplets. Quantification windows were placed ~0.5 cell diameter (7.5 μm) above the connecting axis between the cells. Correlations were calculated using the first derivative of fiber density dynamics. Non-genetic phenotypic cell contractile heterogeneity, which we artificially introduced to simulations, and was necessary to identify cell-ECM-cell communication, is an inherent cellular property. Thus, we aimed at extending our simulation results of distinguishing between pairs of communicating cells to experimental data. The first temporal derivative was sufficient to remove non-stationarity effects in ECM remodeling fluctuations of single and pairs of communicating fibroblast cells (Fig. S5), so we could avoid the second derivative in the correlation analysis of experimental data. Unlike simulations, all cells in a single experiment were embedded in the same fibrous gel possibly inducing spatial ECM correlations in nearby regions that do not necessarily involve cell-ECM-cell communication. Thus, we had to control for spurious correlations that could originate from the proximity between the ECM regions in the same network. This was achieved by including an additional step in the quantification where we corrected the ECM remodeling fluctuations in the quantification window in respect to local regions close to that window to reduce local temporal artifacts that may lead to erroneous correlations (Fig. 4C, Methods).
We performed the same-versus-different pair analysis, but were able to quantitatively distinguish between different pairs of communicating fibroblast cells only after shifting the quantification window 7.5 μm above or below the axis connecting the communicating cells (Fig. 4D, Fig. S6).
A "same" pair had a higher correlation than its corresponding "different" pair in 94% of the matched correlations, for various pair distances (Fig. 4D, Fig. S7A-C) and even for different cell pairs within triplets of cells (Fig. 4E).
To control for potential masking of cell-ECM-cell communication by correlated noncommunication related local ECM remodeling we devised a computational control that we term real-versus-fake pair analysis. We created new pairs where each is composed of one real cell and another fake cell located in the exact same distance as the matched pair of communicating cells, and away from other cells (Fig. 5A, Methods). ECM remodeling correlation between a pair of communicating cells ("real-real") was higher than its corresponding "real-fake" pair in 86% of the matched correlations providing a standardized internal control and further validation that our method truly captures cell-ECM-cell communication (Fig. 5B).
Given our ability to distinguish one pair from a different pair of communicating cells, we wondered whether we can match a cell to its true communication partner when considering all other cells in the experiment (Fig. 5C, Video S5, Methods). For this "matchmaking" analysis, we considered the true communication partner as the one that is connected with a visible band. For each cell, in every cell pair with a visible band, we recorded the first derivative of its fiber density dynamics, correlated it to all other cells in the gel, and reported the rank of its true matching partner. With 50 potential partners, the expected random probability of identifying the communication partner is 1/50 = 0.02 (Fig. 5D). Using the correlation-based matching, the probability of identifying the true communication partner was 0.69 (Fig. 5E). This accuracy of identifying the matching partner dropped to 0.1 when considering quantification regions along the connecting axis between the cells (i.e., no offset in the z-axis, Fig. 5F) and was not attributed to correlated non-communication related local ECM remodeling as verified by careful analysis that considered ECM regions (without cells) located close to each other ("fake" pairs, Fig. 5G,  and non-cellular-related local deformations in the ECM ("real-fake", y-axis). Quantification windows were placed ~0.5 cell diameter (7.5 μm) above the connecting axis between the communicating cells, no offset in X/Y-axis. Correlations were calculated using the first derivative of fiber density dynamics. N = 48 "real-real" cells at cell pair distances ranging 4-10 cell diameters. "Real-real" pair had a higher correlation than "real-fake" pair in 86% of cases. Wilcoxon sign-rank p-value < 0.0001. (  cells. Same-versus-different: N = 35, 61% "same" > "different, p-value < 0.0001. Real-versus-fake: N = 29, 60% "real" > "fake", p-value not significant. Matchmaking analysis: N = 70, correct matching probability = 17%. (C) Dead cells with the presence of live cells. Same-versus-different: N = 17, 55% "same" > "different, p-value < 0.0001. Real-versus-fake: N = 17, 54% "real" > "fake", p-value not significant. Matchmaking analysis: N = 34, correct matching probability = 15%. (D) Live and dead cells. Same-versus-different: N = 17, 61% "same" > "different, p-value < 0.0001. Real-versus-fake: N = 11, 62% "real" > "fake", p-value not significant. Matchmaking analysis: N = 22, correct matching probability = 14%.

Decoupling band formation and communication sensing
Our finding that cell-ECM-cell communication can be identified quantitatively only when shifting the quantification window away from the connecting axis between the cell centers, We conclude that our method's ability to measure communication is improved when correlating quantification windows located slightly away from the connecting axis in-between the communicating cells, while the densest fiber band is formed directly along the connecting axis and that cell-ECM-cell communication exists even when no visible band is formed between the cell pair. These results challenge the current notion that band formation is the hallmark of cellcell communication in fibrous gels. This analysis demonstrated a sweet spot for measuring cell-ECM-cell communication (Fig.   S14A). Regions with low fiber density, that were mostly located away from the connecting axis, as well as denser band regions, had low same-versus-different discrimination. Analysis of the temporal change in fiber density revealed that while the band intensity is maximal along the connecting axis between the cells, the magnitude of the changes in ECM intensity is similar along the connecting axis and slightly above it (Fig. S14B). Because the communicating signal is measured as the change in ECM intensity relative to the background intensity (Fig. 4C), the 'signal to noise ratio' of the measurable communication signal is higher off the connecting axis, making it a more sensitive measurement for cell-ECM-cell communication. Further increasing the Blebbistatin dosage to 150 nM verified that contractility is required for this mode of mechanical communication (Fig. S15B). Cumulatively, these results suggest that mechanical signals propagate from one cell to another even upon massive reduction in their contractility and thus cells can mechanically communicate even with reduced contractility levels that are not sufficient to form a visible band.

Identifying leader and follower in pairs of simulated communicating cells
We next asked if we could use our approach to quantify asymmetric interactions between the communicating partners. More specifically, can we identify which cell in a communicating pair is more dominant or influential? Our quantitative interpretation of "influential" is that past ECMremodeling fluctuations of one cell are predictive of the future ECM-remodeling of its communicating partner. Such temporal order defines a leader-follower relation. We performed simulations where a "follower" cell "imitates" the previous contraction of its influential "leader cell". In other words, the contraction of both cells are determined by the leader cell with a time lag of one step (Fig. S16A). The simplest approach to quantitatively identify a temporal order is by cross-correlation with time lags, where the correlation between two time-series is calculated for a given lag, and the time-lag that leads to the maximal correlation determines the temporal order. This analysis successfully identified which simulated cell was the leader and which one was the follower with an accurate lag time (Fig. S16B, Fig. S17A). To evaluate the robustness of using correlations to identify leader-follower relations we simulated cell pairs where the follower contraction was composed of an independent component, and a component that is dependent on its leader: determined with equal contribution from its intrinsic "decision" and the extrinsic influence by its leader (Fig. S16D, Fig. S17B). In a second validation we tested whether we can identify leader/follower when the follower contracts more than its leader. This leads to increased ECM remodeling that might propagate to the leader cell and mask its influence on the follower. We set to 1 (follower is copycatting the leader's contraction), and introduced ≥ 1 -the fold increase of the follower's contraction: ( ) = * ( − 1). The correlation was nearly identical for ≤ 1.2 (20% increase in follower contraction), and the first mistaken prediction of the follower/follower assignment occurred for = 1.2, for 1 out of 7 simulated pairs (Fig. S16E, Fig. S17C). We further validated these results by pairing leaders or followers to cells from other simulated communicating pairs to create artificial pairs of cells that did not interact with one another (Fig. S18). However, cross-correlation analysis did not identify leader-follower relations in experimental data (Fig. S19). This inability to identify leaderfollower relations in experimental data could be attributed to lack of sufficient temporal resolutionthe response of the follower cell may occur in time scales faster than the 5 minutes temporal resolution in this study. Another alternative is that our method is not sufficiently sensitive to measure the subtle ECM fluctuations that distinguish leader from follower cells.
There is always the possibility that fibroblast cells do not form leader-follower communication patterns, perhaps due to insufficient heterogeneity that can be resolved by assessment of cells pairs from mixed genetic backgrounds. These possibilities are left to be explored in future studies. Cumulatively, these results verified that cross-correlation of ECM remodeling fluctuations can robustly identify leader and follower cells in simulated communicating cells but not in our experimental data. probably due to lower 'signal to noise ratio' (Fig. S14). Moreover, dose-dependent Blebbistatin experiments revealed that while cellular contractile force is required for communication, cells are still able to communicate even after substantial reduction of their contractility, where fibrous bands never form (Fig. S15). These results align with previous studies suggesting that myosin IImediated contractility acts as an inhibitor for cell-to-cell mechanical communication by arresting passive mechanical force transduction through the cellular cortex during collective cell migration [33][34][35][36] .

Discussion
We aimed at identifying leader-follower relationships within pairs of communicating cells. We simulated a situation where one cell influenced its communicating partner, determining its future contractions and demonstrated that our method can robustly identify the leader and follower cells from the ECM-remodeling fluctuations. This was even possible in challenging scenarios where followers were only partially influenced by their leader or contracted up to 20% more than their leader ( Fig. S16C-E). These results set theoretical limitations on our method's capacity to decouple cell contractility from leadership status due to increased ECM remodeling that can propagate to the leader cell and mask its influence on the follower. We were not able to identify leader-follower relations in our experimental data, either because of lack of sensitivity in the given experimental conditions (e.g., time resolution, subtle "leadership" signal in cell pairs from

Finite element simulations of cells contraction in fibrous networks
We used our previously described computational finite element model of one or two cells contracting in two-dimensional fibrous networks 12,30,31 . We used the finite element software

Fiber network architecture
We used Matlab R2018b to construct the network geometry and architecture as previously described 12 . The networks were designed to optimize the fiber orientations distribution toward uniformity (i.e., isotropic) and toward homogeneous fiber density. Briefly, we devised a random process to create network geometries, as we previously reported in 12 . The process starts from uniformly scattering nodes in a circular domain. The nodes were then connected by fiber elements by considering an objective cost function which controls the fiber length, fiber connectivity (i.e., the mean number of fibers intersected at each node) and the angle between fibers connected at each node. The network coordinates spanned from -2 to 2 (AU) in X-and Yaxis (Fig. S20B). Cell centers were located along the X-axis, with a cell diameter of 0.08 (Fig.   S20B). The mean fiber thickness was 0.2 AU, and mean fiber length was 20 AU, fitting typical fiber density-to-length ratio for collagen\fibrin gels [43][44][45] . The cell diameter before contraction was set to 0.08 AU, so that the cell diameter/mean fiber length ratio was 4:1, a typical ratio for fibroblast cells embedded in fibrin gel 44 . The average connectivity of the network was set to eight, to balance the tradeoff between the finite element software numerical stability and physiological relevance.

The mechanical properties of the simulated fiber networks
Fibers were connected to one another by nodes, which acted as a freely rotating hinge, allowing for a rotation of the fibers without resistance. The fibers were modeled as linear truss elements, undergoing uniaxial tension or compression, without bending. They were characterized by nonlinear behavior typical to ECM fibers (such as collagen), including buckling under compression 6,8 and stiffening under tension 4 . We represented the buckling of the fibers by an elastic modulus which was ten times smaller at compressive strains exceeding 2% relative to the elastic modulus at small strains (-2% to 2%). Stiffening was achieved by an exponential increase in the elastic modulus for tensile strains larger than 2% 12,44,46 . In all simulations, the outer boundary of the network was fixed for translations and rotations.
To assess the validity of the method for fibrous ECM networks with varying mechanical properties, we performed additional simulations in which we tested the effect of changing the stiffening behavior to account for a material with a relative linear behavior until 15% strain, representative of fibrin gels (Fig. S4A), based on, for example 47 . where > 1, thus imitating the leader with a certain augmentation.

Experiments Cell culture and chemical reagents
Swiss 3T3 fibroblasts stably expressed GFP-actin (obtained as gifts from S. Fraser, University of Southern California, Los Angeles, CA) were cultured in DMEM supplemented with 10% fetal bovine serum, nonessential amino acids, sodium pyruvate, L-glutamine, 100 units/ml penicillin, 100 µg/ml streptomycin, and 100 µg/ml Neomycin, in a 37•C humid incubator.
Murine Hras mutated-transform cells have been generated by infecting normal epithelial cells of the tongue with mutated-HRAS and shTP53 as described in 48 .

Image analysis and quantifications Preprocessing live imaging data
We first applied Fiji's 49 "Bleach Correction" on the raw fiber channel with the "Histogram Matching" option. On the actin (cell) channel we applied a median filter (radius = 2) followed by a Gaussian blurring filter (sigma = 2), before segmenting the cells over time using the "3D Objects Counter" Fiji's plugin 50 (Threshold = 15, Size filter > 400). The cell's center coordinates in 3D for each time frame was recorded and used for cell tracking. Custom Python code was used for the cell tracking, by identification of cells to track in the first time frame and simply assigning the nearest cell (in 3D Euclidean distance) in the next time frame to construct the trajectory. Shorter trajectories were recorded for cells that moved beyond the field of view during imaging. This simple approach was sufficient thanks to the sparsity of the cell seeding.

Transforming 3D images for visualization and quantification
To visualize and quantify the 3D band between a cell pair we transformed the image to a new coordinate system that is defined in relation to the spatial relation between the pair. We We implemented custom Python code to quantify the ECM density between a pair of communicating cells. First, we transformed the 3D axes to XY slices and Z, replicating the visualization without the Z-interpolation and without the slices averaging. Second, we performed another transformation, rotating the image onto the connecting axis between the cells to reach a common Z-axis. This transformation generates a 3D image where the Z-axis is perpendicular to the XY-axis that is perpendicular to the connecting axis between the cell's centers. This property allows us to move axially in relation to the 3D line connecting the cells. The whole process is depicted in Fig. S20C.

Manual filtering of defected image frames
A small fraction of frames in a few experiment locations had imaging-related artifacts that hampered our ability to accurately segment the cell and quantify ECM densities. These artifacts included incorrect cell segmentation, dark areas due to imaging malfunctioning of the microscope and "light waves" (Video S7) that may have been the result of an air or water bubble trapped in the lenese's immersion oil. To include these experiments in our analysis we manually identified and recorded defected frames that had these artifacts and considered only valid time frames (without this artifact), when computing fiber density and correlations.

Manual annotation of cell pairs with or without a visible band of increased density
For analysis we considered cell pairs with pair distance ranging at 60-150 μm (4-10 cell diameters, assuming fibroblast diameter of 15μm). Based on previous studies that focused solely on cell pairs with a visible band of increased density between them 6,9,11,46 , we partitioned our dataset to cell pairs that formed and ones that did not form a visible band of denser fibers between the cells. This partition was performed manually by visual assessment of the pixel intensity along the full length of the connecting axis between the cells at the end of imaging.
Visually apparent bands appeared in approximately 60% of the imaged cell pairs.

Quantification window size
To quantify the local ECM density we used a quantification window of the size of a cell diameter in all axes, in 2D simulations (0.08 AU) or 3D in experiments (15 μm). This size was set to optimize the tradeoff between including sufficient data versus too much irrelevant data within the window. The number of pixels defining the quantification window (cell diameter in simulations, 15 μm in experiments) were calculated according to the transformed image resolutions. The same scale was also used upon shifting the quantification windows.

Quantifying local fiber density in simulations
The local fiber density was calculated as the accumulated fiber volume within the quantification window. We assume that the fiber volume is preserved even when the fiber is deformed.
However, this property does not hold in the simulated 2D representation of the fibers where their buckling property reduces the simulated fiber lengths. This is an inherent limitation of simulating a 3D process in 2D. To overcome this limitation we normalize each fiber to its initial length before summing the fiber length in the quantification window. More specifically, we considered two scenarios (Fig. S20E). (1) For the case where the fiber was located exclusively within the quantification window, its length at the onset of the simulation was used for quantification. (2) For the case where the fiber was not located exclusively within the quantification window (i.e., crossing the window boundaries), we used only the sub-fiber within the quantification window while adjusting to the full fiber length at the onset of the simulation: where ( ) is the fiber volume at time , ℎ ( ) is the length of the sub-fiber within the quantification window at a time , ℎ ( 0 ) is the overall fiber length at the onset of simulation and ℎ ( ) is the overall fiber length at time .
The fiber density within a quantification window was defined as the accumulated fiber volume within it. For single cells the mean fiber density of four windows, above, below, to the left and to the right of the cell was recorded (Fig. S20F).

Quantifying local fiber density in experiments
The mean fluorescent fibrin channel intensity was used as a proxy of fiber density within the transformed 3D quantification windows (see earlier). Quantification windows with over 5% of pixels extending beyond the image boundaries were marked as "invalid" and were excluded from further analysis. Quantification windows for single cells were calculated similarly to simulations, but in 3D, averaging the mean intensity in 32 orientations (see earlier).

Normalizing the local fiber density in simulations and experiments
To enable quantitative comparison across experiments and between simulations and experiments, we normalized the fiber density to its z-score -the number of standard deviations away from the mean background fiber density at quantification windows that were not influenced by the cells.  (Fig. S20J).

Correlation-based analyses
For correlation-based analysis we considered the longest sub-sequences with continuous valid time frames and mutual timestamps (Fig. S20K). We considered only cell pairs with mutual subsequences of at least 15 (temporal resolution = 15 minutes) or 50 (temporal resolution = 5 minutes) time frames. Correlation was calculated for the second derivative (simulations) or first derivative (experiments) of the fiber density dynamics according to stationarity criteria to avoid high correlations stemming from the monotonic increase of the fiber density (simulations and experiments) and its derivative (simulations) (Fig. S1, Fig. S5). We determined the first/second derivative detrending according to two stationarity tests Kwiatkowski-Phillips-Schmidt-Shin (KPSS) 51 and Augmented Dickey Fuller (ADF) 52 . The null hypothesis in the KPSS test is time series stationarity, while in the ADF test is time series non-stationarity. Temporal correlations were calculated using Pearson correlation on the derived time series.
controls were performed by manually annotating "fake" cell pairs, and analyzing them while "following" the corresponding "real" cells (Fig. S8A, cyan cells). Specifically, the quantification windows' locations followed their corresponding cell pair's motion by shifting in X' and Y' while maintaining a fixed distance from the real pair.

Assessing sensitivity to temporal resolution
To examine the sensitivity of our method to the temporal resolution we performed same-versus-

Pooling data across experiments for statistical assessment
Each experiment was analyzed independently to avoid erroneous relations stemming from correlating two cells in different fibrous networks. Such correlations would be inherently lower when correlating ECM fluctuations between different networks versus in the same fibrous network, due to global network remodeling. Thus, correlating ECM fluctuations between different experiments will lead to erroneously optimistic results, which we avoided here by considering all possible combinations of "same"/"different", "real"/"fake", and matchmaking for each experiment independently. After analyzing each experiment independently we pooled all the results across experiments for statistical assessment. The non-parametric Wilcoxon signedrank test was used for statistical analysis testing the null hypothesis that a distribution, such as of "same"-"different" correlations, was distributed around zero.
Summary of all measurements over all experimental conditions are reported in Table S1.

Leader-follower analysis with cross-correlation
Finite-element simulations with a predefined leader and follower were examined using crosscorrelation analysis, measuring the correlation between two time series under different time lags.
We generated simulations where one cell in a pair was predetermined as the "leader" and its communication partner as the "follower". The "follower" cell "imitated" the "leader" cell contraction in the previous time step, thus the contraction of the "follower" cell is lagging one simulation round behind the "leader". The magnitude of influence that the leader had on the following was defined with the parameter : In a second simulation we tested whether our approach can identify "leader" and "follower" even when the "follower" contracted more than the "leader", but was still copycating the leader's contraction with a time lag of 1 simulation round. The fold increase of the "follower's" contraction was defined using a parameter : Cross-correlation analysis was performed by comparing different time-lags to the cells' time series and evaluating in relation to the simulation ground truth. Cross-correlation was also performed on experimental data, where ground truth was not available, but did not identify any leader/follower relations, either due to lacking temporal resolution or lack of leaders/followers in our dataset.

Data
See Tables S2-3 for detailed information regarding all simulated and experimental data in this work.

Software and data availability
Processed data and source code are publicly available at https://github.com/assafna/cell-ecmproject. The data include the processed ECM remodeling fluctuation time series for each cell, in the key simulations and experiments. The source code to perform all analyses presented here is included.      "Same" pair had a higher correlation than "different" pair in 56% of the matched correlations. Wilcoxon signed-rank p-value < 0.0001. Bottom panel (Z = -0.5, XY = 0): N = 36 cell pairs. "Same" pair had a higher correlation than "different" pair in 88% of the matched correlations. Wilcoxon signed-rank p-value < 0.0001. Figure S7: Identification of communicating cells in different pair distances. Same-versus-different pair analysis. Quantification windows were placed 7.5 μm above the connecting axis between the cells (see Fig. 5E for schematics). Correlations were calculated using the first derivative of fiber density dynamics. All combinations of "same"/"different" were considered. Partitioning the data to cell pairs intervals of 4-6, 6-8, or 8-10 cell diameters was performed under the assumption of a mean fibroblast diameter of 15 μm. (A) Pair distance 60-90 μm (4-6 cell diameters). N = 18 cell pairs. "Same" pair had a higher correlation than "different" pair in 93% of the matched correlations. Wilcoxon signed-rank p-value < 0.0001. (B) Pair distance 90-120 μm (6-8 cell diameters). N = 19 cell pairs. "Same" pair had a higher correlation than "different" pair in 96% of the matched correlations. Wilcoxon signed-rank p-value < 0.0001. (C) Pair distance 120-150 μm (8)(9)(10). N = 11 cell pairs. "Same" pair had a higher correlation than "different" pair in 91% of the matched correlations. Wilcoxon signed-rank p-value < 0.0001. Figure S8: Controlling for potential masking of cell-ECM-cell communication by local ECM remodeling fluctuations. (A) Schematic sketch of a "fake-following" pair (cyan) mimicking a real cell pair (green) by repeating the cells' shifts in X and Y axes (green shadow) in a fixed distance from the real pair (cyan shadow). (B) Quantification of same-versus-different analysis for "fake-following" pairs. N = 25 fake cell pairs. "Same" pair had a higher correlation than "different" pair in 46% of the matched correlations. Wilcoxon signed-rank p-value < 0.0001. (C-F) Cell-ECM-cell communication is more prominent than the masking by local ECM remodeling fluctuations. Same-vs-different pair analysis of cell pairs and their corresponding "fake-follower" cell pair. Pair distance 60-150 μm, correlations were calculated using the first derivative of fiber density dynamics. All combinations of "same"/"different" were considered for both the cell pairs and the "fake" cell pairs independently. Left panels compare the distribution of sameversus-different pair analysis of cell pairs ("real") versus "fake-follower" pairs ("fake"). Right panels compare matched same-vs-different pair analysis for each cell pair ("real") and its corresponding "fakefollower" pair ("fake"). Panels C-D and E-F are the same analysis with low (15 minutes) versus high (5 minutes) temporal resolution correspondingly. In panels C and E the quantification window is located at the focal plane between the cell pair, while in panels D and F 7.5 μm above the focal plain between the cell pair. (C) Quantification window at the focal plane between the cell pair. Temporal resolution = 15 minutes. N = 20 matched pairs. Left: "same" > "different" in 54% ("real") and 62% ("fake") of the sameversus-different pair analysis. Wilcoxon signed-rank p-values < 0.0001. Right: "real" > "fake" in 43% of the matched same-versus-different pair analysis. Wilcoxon signed-rank p-value < 0.0001. (D) Quantification window 7.5 μm above focal plane between the cell pair. Temporal resolution = 15 minutes. N = 25 matched pairs. Left: "same" > "different" in 98% ("real") and 46% ("fake") of the same-versusdifferent pair analysis. Wilcoxon signed-rank p-values < 0.0001. Right: "real" > "fake" in 90% of the matched same-versus-different pair analysis. Wilcoxon signed-rank p-value < 0.0001. (E) Quantification window at the focal plane between the cell pair. Temporal resolution = 5 minutes. N = 19 matched pairs. Left: "same" > "different" in 56% ("real") and 57% ("fake") of the same-vs-different pair analysis. Wilcoxon signed-rank p-values < 0.001. Right: "real" > "fake" in 54% of the matched same-vs-different pair analysis observations. Wilcoxon signed-rank p-value < 0.05. (F) Quantification window 7.5 μm above focal plane between the cell pair. Temporal resolution = 5 minutes. N = 14 matched pairs. (G) "Same" > "different" in 100% ("real") and 67% ("fake") of the same-versus-different pair analysis. Wilcoxon signed-rank p-values < 0.0001. (H) "Real" > "fake" in 100% of the matched same-vs-different pair analysis observations. Wilcoxon signed-rank p-value < 0.0001.     Mean fraction of higher "same", correlation between communicating pairs, versus "different", correlation between one cell from that pair and another cell from a different communicating pair for systematic offsets in Z and XY axes. Red 'x' marked that the null hypothesis that "same" -"different" correlations are distributed around zero was not rejected with p-value ≤ 0.05. Mean fraction of higher "same", correlation between communicating pairs, versus "different", correlation between one cell from that pair and another cell from a different communicating pair for systematic offsets in Z and XY axes. Red 'x' marked that the null hypothesis that "Same" -"different" correlations are distributed around zero was not rejected with p-value ≤ 0.05. Each data point represents the mean fiber density in the quantification window (x-axis) and the fraction of cell pairs with "same" correlation higher than "different" correlation (y-axis) for different locations of the quantification window. The data here was derived from matched bins in Fig. 7B (x-axis) and Fig. 7C (yaxis) in a region that spans 1 cell diameter away (15 µm) from the quantification plane defined adjacent to the cells in the connecting axis between the communicating cells: all bins in the range -1≤XY≤1 and -1≤Z≤1. Pearson correlation coefficient = -0.26, p-value < 0.0001. Note the stiff drop in cell-ECM-cell communication around fiber density of ~3-4 z-scores. (B) The band intensity is maximal along the connecting axis between the cells. The magnitude of the changes in ECM intensity is similar along the connecting axis and slightly above it, where communication is optimal. Top: imaging at temporal resolution of 15 minutes per frame (left-most and right most are duplication of Fig. 7B and Fig. 7C correspondingly), bottom: imaging at temporal resolution of 5 minutes per frame (left-most and right most are duplication of Fig. S12A and Fig. S12B correspondingly).       Table S1: Summary of experimental results, measuring cell-ECM-cell communication in different cell systems, imaging and experimental conditions. Same-versus-different, real-versus-fake and matchmaking results were recorded for each condition. Table S2: Simulation data table. These data include single cells and cell pairs, heterogeneity and leader / follower in pair distances 5, 7 and 9 cell diameters. (1) Whether the experiment involved single cell or cell pairs, (2) Which cell line (fibroblasts/cancer) or beads were used, (3) Whether cells were treated with Blebbistatin or whether the (control) experiment included dead cells, (4) Are the cells "real" or is it a "fake" cells analysis, (5) Whether the cell pairs had a dense fibrous band between them, (6) What was the temporal resolution (in minutes), (7) What was the pair distance between the cells (4-6, 6-8 and 8-10 cell diameters, where cell diameter is estimated at ~15µm), (8) What was the Blebbistatin dosage (in nM), (9) Whether the experiment consisted exclusively dead cells, (10) What was the number (N) of observations (mostly cell pairs) in the category.

Supplementary video legends
Video S1: Finite-element simulation of a cell pair contracting in a 2D fibrous network. Pair distance = 7 cell diameters. Both cell contracts 1% for 50 steps (Methods). Video S7: Time lapse imaging of a cell pair with a manually identified illumination artifact in several time frames that were excluded from analysis.