Inference of long-range cell-cell force transmission 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 long-range cell-cell force transmission 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 relied 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, once strong enough, can form a visible fibrous band of aligned and dense fibers coupling neighboring cells mechanically, and influencing the cells' internal molecular state 15 and active response 16 . This form of long-range cell-cell force transmission through the ECM (termed here cell-ECM-cell communication) was shown to coordinate various biological processes, including tissue injury 16 , fibrosis 17 , vascular assembly, capillary sprouting 1,14,18 , tissue folding 19 , and cancer invasion and metastasis 5,9 . In-vivo, fiber alignment bands can serve as ECM 'tracks' for cell migration with potential roles in wound healing, cancer metastasis and fibrosis 20,21 .
The dictionary definition of communication is "the imparting or exchanging of information" (Oxford Languages) or "a process by which information is exchanged between individuals through a common system of symbols, signs, or behavior" (Merriam Webster). Measuring the transfer of forces through the ECM during cell-ECM-cell communication is challenging and has been typically achieved indirectly by measuring changes in the density, alignment or displacement of the remodeled fibrous ECM between the cell pairs resulting from the active contraction of cells 4, 5, 9-11, 13, 22-26  Here, we present a new computational method to quantify the transmitted ECM signal in between neighboring 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 neighboring 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 neighboring pairs 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.   1B-C). While for single cells the fiber densification faded to the background level at a characteristic distance of approximately 2 cell diameters, the bands between pairs of cells were characterized with increased fiber density that extended further away from the cells and were increased for cell pairs that were located closer to one another (smaller pair distance, Fig. 1C).
The increased fiber density between cell pairs in comparison to single cells could potentially be used to quantitatively decouple the contribution of the interaction of the cell with the ECM (cell-ECM) and the mechanical information transmitted from the communicating partner (cell-ECMcell) that together form the visible band between the cells (Fig. 1D). Specifically, we hypothesized that the information encoded by the cell-ECM-cell component that is unique to a specific communicating cell pair can be used to measure the long-range mechanical communication between cells. Fibers in between pairs of cells (red and purple circles) are remodeled by the integrated mechanical forces that both cells exert on the ECM. Colored arrows depict the magnitude of the force experienced in a specific location in the fibrous gel that are generated by the two cells. (B) Quantitative visualization of representative simulated cell-ECM interactions (left) and cell-ECM-cell communication (right) at the onset ("begin") and after ("end") 50% cell contraction. Color code encodes the ECM density in z-score, the number of standard deviations above background levels. Scale bar is one cell diameter. (C) Quantifying fiber densification as a function of the distance around simulated single cells ("single cells") and as a function of the distance between cell pairs ("pair distance"). Single cells (N = 7). Pair distance of 4 (N = 19), 5 (N = 20 pairs), 7 (N = 20 pairs), or 9 (N = 19 pairs) cell diameters. Note that window distance is bounded in lower pair distances. Mean and standard deviation are displayed for each distance. (D) A dense fibrous "band" between cell pairs is formed by the combined effects of cell-ECM interactions and cell-ECM-cell communication.

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. We clarify that by denoting "communicating" cells we refer to two cells that are within pulling distance in the same (simulated) ECM network, while "non-communicating" cell pairs are two cells that come from different networks ( Fig. 2A). 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. 2B-C, Video S1).
We expected that the temporal correlation of ECM remodeling fluctuations, measured in quantification windows adjacent to each cell, between communicating cell pairs will exceed that of non-communicating cell pairs ( Fig. 2A). 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. 2C). 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. To overcome this confounder, we removed the temporal trends (detrending) via a second temporal derivative. This analysis pipeline transformed the raw simulated local ECM remodeling fluctuations to a correlation-based measure for cell-ECM-cell communication (Fig. 2D, full details in Methods).
Another confounder was that 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. 2E). 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 make a clear distinction between communicating and non-communicating cell pairs, which improved with increasing standard deviations (Fig. 2E). This distinction was negatively correlated to the pair distance -as cells were placed further apart, it became more difficult to distinguish between communicating and noncommunicating cell pairs (Fig. 2F). Moreover, this distinction was improved when moving the ECM quantification window along the band, toward the communication partner (increasing the window distance), which increased the communication partner's influence (Fig. S1C). Finally, we verified that the correlation signal stems from the contractile activity of the two neighboring cells by demonstrating that the correlation of two contracting cell pairs exceeded the correlation of pairs of one contractile and a second non-contractile passive cell (Fig. S1D).
Altogether, these results established the notion that a correlation-based approach can capture force transmission between actively contractile cells, that temporal correlation of local ECM remodeling fluctuations can distinguish between communicating and non-communicating simulated cell pairs, and that contractile heterogeneity is required to make this distinction. (B) Representative simulation visualization at the onset (top), 25% cell contraction (middle), and 50% cell contraction (bottom). (C) Quantifying fiber densification as a function of simulated contraction steps in between cell pairs. N = 20 cell pairs. Color code (B) and y-axis (C) encode the ECM density in z-score, the number of standard deviations above background levels. The variability in fiber density in is a consequence of the randomness in the network architecture. The pair distance between simulated cell centers was set to 7 cell diameters. (D) Correlation-based pipeline to measure cell-ECM-cell communication in simulations. Raw time series refer to measurement of ECM density in the quantification window (top), normalization in respect to the mean background fiber density at regions that were not influenced by the cells and detrending by second temporal derivative lead to a fluctuating signal (bottom) that is correlated to distinguish between pairs of communicating versus non-communicating cells (left versus right). (E) Distinguishing between communicating and non-communicating 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 noncommunicating 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. (F) 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). Similar to the simulations, we clarify that "communicating pairs" are pairs of cells located next to each other in the fibrous gel and thus are within a pulling distance, whereas "non-communicating cells" are pairs that are not in proximity (beyond 10 cell diameters away, see 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 background quantification windows defined at the onset of the experiment at regions that were not influenced by the cells (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 N = 13 cell pairs at pair distances of 6-8 cell diameters (90-120 μm). Error bars in C-D indicate standard deviation.

Temporal correlation of local ECM remodeling fluctuations defines a signature for cell-ECM-cell communication
Rather than comparing communicating versus non-communicating cell pairs (Fig. 2E), definitive quantification of cell-ECM-cell communication lies in the ability to distinguish between pairs of communicating cells, i.e., whether the ECM fluctuations of one cells pair has a unique signature that can be distinguished from that of a different pair of communicating cells. We tackled this challenge by testing whether the correlation between a simulated pair of communicating cells surpassed the correlation between one cell from that pair and another cell from a different simulated pair of communicating cells located in a different fibrous network (Fig. 4A).
Specifically, we compared the correlation in the second temporal derivative of the local ECM density dynamics between quantification windows adjacent to each cell pair, located in the 'same' pair versus a cell in a 'different' pair, and accordingly coined the term same-versusdifferent pair analysis (schematic in Fig. 4A). Having the "same" correlation exceeding the "different" correlation implies that the correlation between a communicating cell pair is not merely an effect of a similar ECM densification pattern that is common to any pair of communication cells, but rather is indicative of communication unique for the "same" cell pair at test. We considered all possible ordered combinations of triplets of cells that include a pair of communicating cells and a third cell that takes part in a different communicating cell pair (Fig.   4A). We then analyzed these triplets using the "same" (communicating) versus the matched "different" (non-communicating) cell pair correlations. In the corresponding plot (Fig. 4B), each data point represented the correlation of one communicating cell pair (x-axis) to multiple noncommunicating cells (y-axis), each leading to a different value (shown as a vertical line of points in Fig. 4B). Data points below the diagonal y = x (red line in Fig. 4B) indicate that correlation between the communicating pair exceeded that of the non-communicating pair. In simulated cells, a "same" pair had a higher correlation than a "different" pair in 92% of the matched correlations for pair distance of 4 cell-diameters, and this correlation gradually reduced with increased pair distance (Fig. 4B). To assess the validity of the method for fibrous ECM network with varying mechanical properties we also simulated a network with a relative linear-elastic behavior until 15% strain compared with the current fiber stiffening model (Fig. S4A). Sameversus-different analysis was improved for the stiffening model (Fig. S4B), implicating the role of strain stiffening in long-range mechanical communication, in line with our previous studies 12, 27,28 . Altogether, simulated cells communicating with one-another were more synchronized in their ECM remodeling fluctuations, and same-versus-different analysis could distinguish between different pairs of communicating cells.
We next aimed at extending our simulation results of distinguishing between pairs of communicating cells in experimental data. Our analysis pipeline was adjusted to enable analysis of 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 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). 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) Same-versus-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 "Same" -"different" correlations were distributed around a mean 0: p-value < 0.0001 for all pair distances. "Same" -"different" correlations were different for cell pairs at distances of 4 versus 9 cell diameters (T-test pvalue < 0.001). (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. 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. 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/Yaxis. The same image is also shown in Fig. 3B. Middle: Analysis pipeline: raw time series were normalized in respect to local regions close to their quantification window to reduce local spurious correlations (see panel C), and correlations were calculated using the first derivative of fiber density dynamics, which was sufficient to remove non-stationarity effects in ECM remodeling fluctuations. Full details in Methods. For the displayed cell pair (top) the Pearson correlation coefficient was 0.66 (Z offset = 0) and 0.71 (Z offset = 0.5). 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-versus-different pair analysis for cell triplets. Quantification windows were placed ~0.5 cell diameter (7.5 μm) above the connecting axis between the cells. 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  real-versus-fake pair analysis. "Fake" cell (cyan) coordinates were located on the circumference of a circle (dashed purple), with a radius r from the corresponding "real" communicating cell (dark green), and with the maximal distance (d, orange) in respect to all other real cells (light green). (B) Correlations between ECM fluctuations of pairs of communicating cell pairs ("real-real", x-axis) versus the correlation between a cell 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. (C) Schematic sketch of the matchmaking analysis. ECM remodeling dynamics of a given cell (left) was correlated with 49 other cells (with repeatssee Methods) from an experiment (right). The rank of the correlation with the true communication partner (i.e., the position in the sorted list of all correlations with other cells, marked in green) was recorded. This process repeated for all cells.

Experimental controls validate ECM remodeling correlations as a measurement for cell-ECM-cell communication
As negative controls we included a new set of experiments. Severe reduction in all three communication readouts was measured between pairs of fluorescent beads in size comparable to cells, and in experiments with dead cells (Fig. 6A-B), compared to live fibroblasts (Figs. [4][5] and cancer cells (Fig. S9). To include background ECM remodeling resulting from contraction of other live cells, we co-cultured live and dead cells within the same gel. In the presence of live cells, no communication was measured between pairs of dead cells (Fig. 6C) and between livedead cell pairs (Fig. 6D) -a validation that our measurement for cell-ECM-cell communication is robust to local ECM deformations originating from contractile activity of the two partners and not from other cells in the gel or remote sensing of physical boundaries through the ECM (e.g., 29 ) correspondingly. Summary of all measurements over all experimental conditions are reported in Table S1. "real" > "fake", p-value not significant. Matchmaking analysis: N = 50, correct matching probability = 14%. (B-D) Experimental controls with live and dead fibroblast cells. (B) Dead cells. Same-versusdifferent: 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, pvalue < 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 located far away from the connecting axis, had low same-versus-different discrimination, and this low discrimination was also apparent in denser band regions along the connecting axis. 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 simulation step (Fig. S16A). The simplest approach to quantitatively identify a temporal order is by cross-correlation with time lags, where the correlation between two timeseries 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 where the follower cell contraction is 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 leader-follower 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 [30][31][32][33] .

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

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,27,28 . 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 [42][43][44] . 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 43 . 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,43,45 . 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 46 . 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 47 .

Image analysis and quantifications Preprocessing live imaging data
We first applied Fiji's 48 "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 49 (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,45 , 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 window size was set to optimize the tradeoff between including sufficient data versus too much irrelevant data within the window (Fig. S21). 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.

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) 50 and Augmented Dickey Fuller (ADF) 51 . 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.

Same-versus-different pair analysis
To establish that cell-ECM-cell communication of one cell pair can be distinguished from a second cell pair we tested whether the correlation between a cell pair ("same" pair) surpassed the correlation between one cell from that pair and another cell from a different cell pair ("different" pair) (Fig. 4A). This comparison was termed same-versus-different pair analysis. In this analysis, we considered all combinations of triplets of cells in an experiment, that included one communicating cell pair and another cell that takes part in another communicating pair (Fig. 4A).
The quantification window of each cell in the analysis was always located in relation to its communication partner (Fig. 4A).

Real-versus-fake pair analysis
To control for misinterpreting non-communication related local ECM remodeling correlations as cell-ECM-cell communication we devised a standardized internal computational control that measures whether the correlation between a cell pair ("real" pair) surpassed the correlation between one cell from that pair and another "fake" cell (forming together a "fake" pair). The "fake" cell location was determined by following two constraints. The distance between the "fake" and its "real" partner is equal to the distance between the communicating cell pair, while maximizing the distance to other real cells. The first constraint was defined to allow comparable distance between the "real" and "fake" cells pairs. The second constraint was defined to minimize the effect of other cells in the vicinity. See depiction in Fig. 5A.

Matchmaking analysis
For each cell that takes part in a cell pair we tested our ability to identify its matched communication partner from all the other cells in that experiment (Fig. 5C) The probability of identifying the true communication partner was recorded at the matchmaking score.

Internal control to validate that our analysis is not an artifact of non-communicationrelated correlated local ECM remodeling fluctuations
The contraction of the cells in the gel lead to local ECM correlations in the fibrous network, even in cell-free areas. To further verify that our results in the same-versus-different pair and matchmaking analyses were not merely an artifact of these local ECM correlations in the fibrous network, we compared the correlation of communicating cell pairs to ECM remodeling correlations in quantification windows that were placed in cell-free, fibrous areas. The intuition behind this control experiment was to consider correlations in quantification windows that measure the proximity-component, or the local mutual ECM fluctuations, without the influence of the communicating cells. These 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.    For 'stiffening', "same" > "different" in 72% of the matched same-versus-different pair analysis (92%, 79%, 70% and 68% for pair distances of 4,5,7 and 9 cell diameters, correspondingly), whereas in for 'linear', "same" > "different" in 61% of the matched same-versus-different pair analysis (76%, 68%, 71% and 43% for pair distances of 4,5,7 and 9 cell diameters, correspondingly).  "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. Figure S9. Cell-ECM-cell communication in Murine cancer Hras mutated cells. Same-versus-different: N = 48, 94% "same" > "different, p-value < 0.0001. Real-versus-fake: N = 48, 86% "real" > "fake", p-value < 0.0001. Matchmaking analysis: N = 96, correct matching probability = 72%. Figure S10. Sensitivity of same-versus-different pair analysis to the temporal resolution. Pair distance 60-150 μm, quantification windows were placed 7.5 μm above the connecting axis between the cells (see Fig.  7E for schematics). Correlations were calculated using the first derivative of fiber density dynamics. (A) Identification of communicating cells in higher temporal resolution of 5 minutes per frame. Same-versusdifferent pair analysis. All combinations of "same"/"different" were considered. Pairs with a visible band. N = 14 cell pairs. "Same" pair had a higher correlation than "different" pair in all (100%) of the matched correlations.   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.