Automatic detection of spatio-temporal signaling patterns in cell collectives

Cell signaling is a finely regulated process that governs fundamental aspects of a cell’s life, its metabolism, responses to environmental cues, and cell fate determination. Time-lapse fluorescence microscopy has become a routine experimental modality that has greatly furthered our understanding of these processes and revealed the role of cell–cell heterogeneity and signaling dynamics (Purvis and Lahav, 2013; Albeck et al., 2013; Ryu et al., 2015). Live-cell imaging of fluorescent biosensors has also demonstrated that signaling propagates across cell collectives, a phenomenon that is becoming another important organizational principle in biology. Even though we are only beginning to understand the rules that govern the emergence of spatial correlations from single-cell signaling, their functional significance has been reported in various systems. Notably, the coordinated activation of extracellular signal-regulated kinase (ERK) in spatial clusters of cells is instrumental in the maintenance of epithelial homeostasis (Aikin et al., 2020; Gagliardi et al., 2021; Pond et al., 2022), acinar morphogenesis (Ender et al., 2022), osteoblast regeneration (De Simone et al., 2021), cell cycle progression (Hiratsuka et al., 2015), and in the coordination of collective cell migration (Aoki et al., 2017; Hino et al., 2020). Calcium waves observed in 2D monolayers trigger and facilitate the extrusion of oncogenically transformed cells from the epithelium, thus also contributing to its homeostasis (Takeuchi et al., 2020).

Waves rely on an active process in which individuals relay the signal from their neighbors (Talia and Vergassola, 2022). They are an important class of self-organizing phenomena that have emerged repeatedly at different biological scales, arguably because they can quickly transmit information at long distances throughout a community (Gelens et al., 2014; Wedlich-Söldner and Betz, 2018). At the cellular level, voltage-gated sodium channels mediate the action potential that propagates in neurons (Hodgkin et al., 1952). A rapid initial development phase in large animal eggs requires long-distance coordination of mitotic events, which can only be achieved by signaling waves (Chang and Ferrell, 2013). At the cell population level, waves rely on cell–cell communication as is the case of cyclic adenosine monophosphate waves in Dictyostelium discoideum (Alcantara and Monk, 1974), calcium waves in animal development, adult animal tissues, and plants (Robb-Gaspers and Thomas, 1995; Gilland et al., 1999; Choi et al., 2014), depolarization waves in the heart that control heartbeat (Noble, 1962), mechanical waves (Serra-Picamal et al., 2012), and genetic waves in the presomitic mesoderm (Tsiairis and Aulehla, 2016). Wave patterns also occur in communities of independent multicellular organisms. Giant honeybees flip their abdomens in a coordinated fashion that results in shimmering waves that repel an intruder that endangers the colony (Kastberger et al., 2014a).

Despite waves being so common in biology, a general and easy-to-use computational tool to identify them and extract useful statistics for further analysis is still missing. In parallel, other fields like geographic information science and ecology have accumulated a multitude of statistical methods to analyze patterns in weather prediction, epidemiology, animal migration, and voting preferences (Moore and Carpenter, 1999; Lawson, 2006; Pfeiffer et al., 2008; Bivand et al., 2013). Here, we build upon methods for spatio-temporal clustering and analysis of point processes and develop a computational tool for the Automatic Recognition of COllective Signaling (ARCOS) events and their characterization with various statistics. The open-source code is freely available as R and Python packages, along with extensive documentation. We also developed a plugin for napari—a fast, interactive, multidimensional image viewer for Python (Sofroniew et al., 2022). The plugin provides a graphical interface to ARCOS and facilitates an intuitive discovery environment to quantify collective phenomena in time-lapse images.

In the following sections, we explain the ARCOS algorithm and apply it to several biological systems. First, we used an optogenetic actuator that activates the mitogen-activated protein kinase (MAPK) network to generate synthetic waves of ERK activity in vitro in an MCF10A epithelial monolayer. This way we could establish a ground truth to test the algorithm in a controlled setting. We then studied mammary epithelial cells expressing different oncogenic mutations. We found that ERK activity waves are exacerbated by oncogenic mutations and that the wave propagation depends on matrix metalloproteinase (MMP)–mediated intercellular communication. Further, we quantified apoptosis-induced waves of ERK activity, calcium waves in renal cells, and collective shimmering in a colony of giant honeybees. Finally, we tackled 3D geometry and identified ERK waves in mammary acini in vitro.

To explicate the core ARCOS algorithm, let us consider a radially expanding wave, a “firework,” as shown in Fig. 1 A, where we schematically depict three consecutive frames of a developing activation cluster. Such events have been observed, for example, in epithelial monolayers where a correlated activation of ERK is typically initiated by an apoptotic cell (Aikin et al., 2020; Gagliardi et al., 2021; Pond et al., 2022). In the first step, we apply a spatial clustering algorithm, dbscan (Ester et al., 1996; Campello et al., 2013; Hahsler et al., 2019), to every frame independently to identify pockets of active cells. The clustering requires two parameters: the search radius ε and the minimum cluster size minClSz within that radius. We recommend setting ε to ≈1.5× the first nearest neighbor distance (NND), and minClSz to the dimensionality of the problem plus one, i.e., 3 for 2D, 4 for 3D, etc. (Ester et al., 1996).

In the second step, spatial clusters identified in the first step are linked between frames based on their relative distance. The first time frame contains only a single five-cell spatial cluster #1 that becomes the seed of a new collective event (CE). The second frame also contains a single but bigger cluster #2. To link clusters #1 and #2, we calculate the NND between cells that comprise both clusters and compare it to the search radius ε_prev used for tracking clusters over time. For simplicity, we set ε_prev = ε. In our case, at least one cell from cluster #2 is within ε_prev to cells in cluster #1. Hence, cells in cluster #2 inherit the identifier of the seed cluster #1. The third frame contains four spatial clusters #3–6. Based on NNDs between these cells and cells from the previous frame, we conclude that only clusters #3, #5, and #6 can be linked to the growing CE from the previous frame. Cells from cluster #4 are farther than ε_prev and therefore they become a seed of a new CE.

Thus, the algorithm first identifies spatial clusters in individual frames and then tracks the clusters between frames based on their spatial overlap (Fig. 1 B). To illustrate how this algorithm can integrate into a real-life workflow, consider a minimal pipeline to identify collective ERK activation in a 2D MCF10A epithelium (Fig. 1 C). After acquiring multichannel time-lapse movies, images are segmented and cells are tracked over time. The resulting single-cell ERK activity time series (ERKATS) form a space–time pattern (x,y,m,t), where (x,y) are the centroid coordinates of nuclei, m is the measurement of ERK activity, and t corresponds to time. The input for ARCOS should only contain “active cells,” i.e., spatial coordinates and time, (x,y,t), of cells with the measurement m greater than a threshold.

If the input data for the first step of the analysis does not comprise active cells, a binarization of the measurement is necessary before clustering. We implemented several methods to binarize time series data in ARCOS. The most straightforward approach applies a global threshold to the measurement. However, data from time-lapse microscopy often exhibit high variability in the baseline due to different expression levels of the biosensor. Additionally, high-intensity light and long exposures can lead to photobleaching of fluorescent probes, which may cause long-term trends in fluorescence intensities. To address this, ARCOS offers detrending methods based on a running median or a fit to a polynomial (Fig. 1 D; Materials and methods).

An alternative approach to measurement binarization utilizes local indicators of spatial association (LISAs), which are a class of statistics that quantify local spatial correlation and clustering (Getis and Ord, 1992; Anselin, 1995; Naimi et al., 2019). A LISA statistic applied to individual time frames of the space–time pattern (x,y,m,t) yields a coefficient of local spatial association ρ at every coordinate (x,y). For our purposes, we are looking for pockets where the coefficient is high, which corresponds to regions with high and positively correlated single-cell ERK activity. After thresholding ρ, we obtain a space–time pattern (x,y,t) that is restricted only to cells that form pockets of high activity. We assume these clusters estimate the location of active cells and we use this subset as an input for spatial clustering with ARCOS. LISA-based binarization was performed accurately on a cluster of cells with high ERK activity in the MCF10A epithelium (Fig. 1 E).

ARCOS applied to binarized measurements yields spatio-temporal clusters with identifiers of participating cells. Such a dataset can be further visualized or processed to extract various statistics (Fig. 1 C). Our R/Python implementations offer several helper functions to calculate CEs’ duration, size, growth over time, and center of mass displacement. Visualization of results is best performed in the napari plugin, where the original data can be viewed together with image segmentation results and CEs. The plugin also offers interactive plots to display the statistics and to examine cells’ measurements over time, their binarization, and their assignment to CEs (Video 1).

The observed statistics of CEs can be further tested against a null hypothesis that they do not differ from CEs calculated in a system with an equivalent total activity but with activations occurring randomly in space and time. We provide several R/Python implementations, with varying assumptions about the null hypothesis, to randomize the binarized measurement and to simulate the null distribution, i.e., the sampling distribution of a statistic under the null hypothesis (Fig. S1; Materials and methods; Manly, 2007). To test a one-sided alternative hypothesis that, for example, the mean CE duration in the observed time-lapse is greater than the CE duration in randomly activating cells, we repeatedly apply any of the randomization methods and calculate the statistic for every iteration. The fraction of cases when the durations calculated from the randomized dataset are equal to or exceed the observed duration is the P value, i.e., the probability of obtaining the statistic at least as extreme as observed (given the null hypothesis is true; Davison and Hinkley, 1997; Good, 2005; Phipson and Smyth, 2010). These randomization tests of CEs’ statistics quantify the “unusualness” of observed spatio-temporal correlations and thus form the validation part of our workflow.

Aside from synthetic datasets used during the algorithm development, we also wanted to test the code on a biological system, with the actual variability in protein activity, but in a controlled setting. To this end, we utilized an MCF10A cell line with a stably expressed optogenetic actuator–biosensor system described earlier (Aoki et al., 2017; Fig. 2 A). Using an OptoRAF actuator, we could induce ERK activity via the RAF-MEK-ERK pathway in individual cells with light pulses and measure ERK activity using the cytoplasmic/nuclear translocation of the ERK–KTR biosensor. Additionally, we used a digital micromirror device (DMD) to stimulate desired regions of the epithelium over time, thus generating synthetic collective ERK activity patterns that would mimic various aspects of spontaneous waves (Fig. 2 B). Notably, the ERK activity induced by OptoRAF does not propagate to neighboring cells, which allowed us to precisely control the extent of synthetic ERK activity patterns and to benchmark ARCOS.

In the first scenario, we repeatedly stimulated a growing circular region of the epithelium at 5-min intervals, starting at 11 min, with 100-ms light pulses. We observed a clear ERK activation in the stimulated regions (Fig. 2 C). ARCOS applied to detrended and binarized ERKATS identified a single, growing activity cluster induced by optogenetic stimulation (Fig. 2 C and Video 2). Additionally, several smaller CEs that lasted <5 min or comprised <20 unique cells were detected. They stemmed from spontaneous ERK activations in the epithelium and were discarded after processing.

To verify that the optogenetic stimulation was confined to desired regions, we inspected single-cell ERKATS (Fig. 2 D). Cells from the innermost region received all stimulation pulses after the 10-min mark. ERK activity increased sharply, although a few responses were delayed, and ensued only after the second or even the third pulse. We attribute this delay to biological variability and the fact that the cells are at the border of an already small activation region. Cells farther away from the center responded accordingly to pulses that they were exposed to, although a few cells were activated earlier, likely by the illumination of regions with smaller radii. The ERK activity in cells not exposed to DMD illumination remained at the baseline.

We concluded that our protocol was sufficiently selective to control the stimulation of a growing ERK activation cluster and that ARCOS correctly identified such a synthetic wave. An intuitive visualization of that result is shown in a “spaghetti” plot (Fig. 2 E). Therein, the spatio-temporal evolution of collective ERK activity was projected onto the 2D (x,t) plane. Each trace corresponds to a single cell involved in a CE identified by the algorithm. We confirmed that the collective ERK activation grew over time and within the bounds of the DMD illumination. As a cross-check, we compared the result obtained from ARCOS with an independent quantification of the global spatial correlation per frame using Moran’s I statistic (Gittleman and Kot, 1990; Bivand et al., 2013). Said statistic is commonly used in the field of geographic information science to quantify spatial phenomena. The Moran’s I coefficient increased with the growing collective ERK activity region and peaked when that region was the biggest (Fig. 2 E). This result confirmed a spatially correlated cluster as detected by ARCOS.

We also tried other illumination scenarios to further challenge the algorithm. To test the tracking, we illuminated a pattern with a circle that moved around the field of view (FOV; Fig. 2, F and G; and Video 3). Another common scenario observed in the epithelium prompted us to apply a circular pattern that was slowly splitting into two detached regions (Fig. 2, H and I; and Video 4). In both cases, ARCOS correctly identified the spatio-temporal progression of the optogenetically induced ERK activity region. Due to a longer, 2 h 40 min, acquisition time, we detected additional ERK activity clusters that stemmed from apoptotic events and/or a spontaneous activation. The results from the optogenetic experiments assured us that our image analysis pipeline performed well in a controlled setting with biological variability and that we could apply it to more challenging scenarios.

One of our earlier studies showed that apoptosis of single cells within epithelium triggers waves of ERK activity pulses (Gagliardi et al., 2021). These waves act as a survival signal that protects two to three layers of neighboring cells from further apoptosis, thus contributing to epithelial homeostasis. The ERK waves are trigger waves (Gelens et al., 2014) in which the extruding apoptotic cell activates MMP-mediated cleavage of pro-EGFR ligands. Subsequently, the EGFR activates ERK in adjacent cells (Aoki et al., 2017; Aikin et al., 2020). ERK-mediated activation of myosin contractility might then provide a mechanical input that activates ERK in farther cells through MMP/EGFR signaling. This explains how the ERK wave emerges through a mechanochemical feedback loop (Hino et al., 2020). Given those past results, we were curious how oncogenic mutations that affect the MAPK network might influence the spatio-temporal ERK activation and the molecular mechanism of CEs propagation. Therefore, we explored ERK waves in starved MCF10A cells with oncogenic heterozygous mutations, KRAS^G12V (Konishi et al., 2007) or PIK3CA^H1047R (Gustin et al., 2009), that activate the MAPK/ERK and frequently occur in cancer. In the KRAS^G12V mutant, loss of GTP hydrolysis leads to increased activation of the rapidly accelerated fibrosarcoma (RAF) kinase, thus activating the MAPK/ERK kinase (MEK) and ERK. The H1047R mutation in PIK3CA, the α subunit of PI3K, leads to amphiregulin expression/secretion amplification, which activates the non-mutated MAPK/ERK pathway via EGFR (Gustin et al., 2009; Young et al., 2015).

As expected, at the population level, in starved conditions, without external stimulation, the average ERK activity was higher in the mutants than in the WT cells (Fig. 3 A). Notably, the increased basal activity occurs without ERK activity waves due to apoptosis, which is negligible in the mutants. At the single-cell level, ERK activity still exhibited non-periodic pulses, although their number and shape markedly differed between the two mutant cell lines (Fig. 3 B). WT cells displayed sharp, isolated, and rare ERK pulses as observed earlier (Albeck et al., 2013). In the KRAS^G12V cells, individual ERK pulses were still discernible but often lasted longer and formed “massifs” of ERK activity that resulted in higher pulse frequency (Fig. 3 C). In the PIK3CA^H1047R cells, as already reported (Ender et al., 2022), ERK pulses retained the same shape, duration, and amplitude as in WT cells, but the frequency of these pulses was the highest (Fig. 3 C).

We then wondered how different single-cell ERK activity dynamics translated into spatial correlations in the epithelium. ERK activity was least spatially correlated in WT cells, where CEs (comprising minimum three cells and lasting at least 15 min) were sparse, short-lasting, and usually triggered by apoptotic events. In comparison, KRAS^G12V and PIK3CA^H1047R mutants exhibited a more spatially correlated ERK activity with longer and more frequent CEs (Fig. 3 D and Video 5). Given the quantification of CEs obtained from ARCOS, we quantified their duration and size, i.e., the total number of unique cells in a CE (Fig. 3 E). Even though the median statistics were similar across the cell lines, the distributions varied and confirmed our observations from spaghetti plots: small and short-lasting ERK activity waves in WT cells, and bigger and slightly longer-lasting events in the mutants (Fig. 3 F). Furthermore, the mutant cells triggered more CEs compared with the WT (Fig. 3 G).

Next, we explored the molecular mechanism of CE propagation in the presence of oncogenic mutations. We hypothesized that they were mediated by the same MMPs/EGFR intercellular communication mechanism that we identified for apoptotic waves in WT cells (Gagliardi et al., 2021). Treatment with 10 µM of the MMPs’ inhibitor batimastat reduced the average ERK activity in the mutant cell lines (Fig. 3 A). WT and PIK3CA^H1047R ERKATS became nearly flat, with a strong increase in the fraction of cells without any ERK activity pulses (WT: 24%→65%; PIK3CA: 9%→51%). In contrast, 58% of KRAS^G12V cells still displayed residual pulsatile activity (Fig. 3 C). This result suggests that WT and PIK3CA^H1047R completely depend on intercellular communication for the regulation of ERK activity dynamics, whereas the dependency of KRAS^G12V is only partial. The hampered ERK activity dynamics after MMP inhibition coincided with a strong reduction of CEs’ number, size, and duration in all three cell lines. However, we observed that in KRAS^G12V there are more residual CEs than in WT and PIK3CA^H1047R (Fig. 3, D–G).

An independent quantification with the global Moran’s I coefficient calculated directly from segmented time-series data showed a reduced spatial correlation of ERK activity in all cell lines after batimastat treatment (Fig. 3 H). This confirms less coordinated ERK activation, as quantified with ARCOS, and shows that not only apoptotic ERK waves but also aberrant CEs caused by oncogenic mutations are at least partly mediated by MMPs.

We then compared the observed statistics of CEs to statistics calculated in a system with the same dynamics but without spatial correlations. We performed a randomization test by shuffling positions of entire time series according to Fig. S1 B. This only affected spatial correlations between neighboring cells, while leaving the dynamics intact. After performing 1,000 iterations, we detected CEs that were shorter-lived and involved fewer cells than the observed CEs; thus the one-sided P value for CE duration and the size was <0.001 (Figs. S2 and S3). We concluded that at the significance level of 0.05, we can reject the null hypothesis that the observed statistics of CEs in three cell lines with and without batimastat are no different from those calculated from random events.

To further demonstrate ARCOS’s versatility, we quantified ERK waves in the NRK-52E renal epithelial cells that express the EKAREV-NLS ERK sensor. The majority of the waves in this cell line are triggered by apoptosis, as reported earlier (Aoki et al., 2013). Since waves are bigger and isolated from each other, we used a LISA-based binarization of ERK activity as the pre-processing step for ARCOS. To identify pockets of active cells, we calculated the local G*, a LISA statistic that identifies spots of activity (Naimi et al., 2019), and applied it to individual frames of the (x,y,m,t) space-time pattern from image segmentation. After thresholding the resulting statistic, we obtained the (x,y,t) pattern only with active cells, which we then used in ARCOS (Fig. 4 A and Video 6).

The spaghetti plot confirms that apoptotic events trigger ERK activity waves (Fig. 4 B). We corroborate this result with the global Moran’s I indicator of spatial correlation, which indeed peaks at the time of CEs identified with ARCOS (Fig. 4 B). We also find that the median CEs’ duration of 0.5 h in NRK-52E cells is comparable with all MCF10A cell lines analyzed before (0.58–0.67 h) but the events are bigger with their median size of 16 unique cells compared with five to six in MCF10A WT and mutants (Fig. 4 C). The one-sided P value from 1,000 iterations of the randomization test with shuffling time series’ positions was <0.001 for CE duration and size (Dobrzyński, 2023).

We then tested ARCOS on a different cell line and biosensor by analyzing calcium waves that trigger apical extrusion in the MDCK epithelium (Takeuchi et al., 2020; imaging data, courtesy of Yasuyuki Fujita). Apical cell extrusion is an important homeostatic mechanism that eradicates harmful or suboptimal cells from the epithelium. In the experiment, the Myc-RasV12 oncogene was transiently expressed mosaically within a monolayer of MDCK cells. Due to their harmful potential, RasV12-transformed cells were apically extruded from the epithelium. Notably, before the extrusion, calcium levels acutely increased in the RasV12 cells and radially propagated to neighboring cells. Such waves have been shown to trigger and facilitate the process of cell extrusion.

Since single-cell time series of calcium levels had a well-defined baseline and lacked long-term trends, we used the G* LISA statistic (Naimi et al., 2019) to identify pockets of active cells (Fig. 4, D and E; and Video 7). Spatial correlation calculated for each frame using global Moran’s I coefficient peaks when calcium waves occur, which further strengthens our identification of CEs with ARCOS (Fig. 4 E). We find that a single collective calcium wave propagates radially from the extruded cell, can involve up to 300 cells, and can last up to 2 min (Fig. 4 F). The one-sided P-value from 1,000 iterations of the randomization test with shuffling time series’ positions was <0.001 for CE duration and size (Dobrzyński, 2023).

Even though we developed ARCOS to tackle collective phenomena in cell signaling, our image processing pipeline can also quantify spatio-temporal phenomena in the context of ecology. As a demonstration, we analyzed a previously published time-lapse dataset of a giant honeybee colony (Apis dorsata) that collectively responds to a computer-controlled dummy wasp hovering in front of the honeybee nest (Kastberger et al., 2011). The giant honeybees responded to the wasp by flipping their abdomens in a simultaneous and cascadic way, generating a wave-like visual signal for external observers. The visual effect achieved by this repetitive collective response aims to repel the threat (Kastberger et al., 2008, 2014b). Pulses connected to these visual patterns at the nest surface also lead to vibrations within the nest, and these are suitable for informing nest mates, even on both sides of the nest, about the threat status (Kastberger et al., 2010; Kastberger et al., 2011; Radloff et al., 2011; Kastberger et al., 2012; Kastberger et al., 2013; Kastberger et al., 2014a).

As described above, shimmering waves propagate across the surface of the giant honeybee nest in response to repeated exposure to the dummy wasp presented to the colony on the right side of the FOV (Fig. 4 G). Notably, our quantification confirms an interesting aspect of the shimmering phenomenon. Each vertical streak in the spaghetti plot (Fig. 4 H) corresponds to a wave that propagates from right to left of the FOV. However, we note that each of such streaks consists of “sub-waves” as indicated by different colors within the streak. Such a shimmering wave is typically initiated by “leader” bees on the nest surface, assembled in so-called trigger centers across the nest surface (Schmelzer and Kastberger, 2009; Kastberger et al., 2010; Kastberger et al., 2014a). Video 8 confirms our observation that a single exposure to the dummy wasp induces several smaller waves consecutively triggered by the leader bees scattered across the colony. These trigger centers were primarily arranged in the close periphery of the mouth zone of the nest, around those parts where the main locomotory activity occurs throughout daytime.

Previously, we reported waves of ERK activity during the morphogenetic program of 3D MCF10A acini in vitro (Ender et al., 2022). One crucial stage of this process is the formation of a lumen by apoptosis of inner cells. We found that ERK waves coordinate spatial separation into two domains: the outer cell population, which survives, and the inner cell population, which undergoes apoptosis. Since the ARCOS algorithm can handle any integer dimensionality, we were curious to analyze ERK activity waves in this system. The top row of Fig. 5 A shows a single wave in a maximum projection of the raw ERK-KTR channel. This example CE starts at the surface of the acinus and propagates over time to other cells in the outer layer (Fig. 5, B and C; and Video 9). Only toward the end, some inner cells participate in the event. This observation reflects a more general trend that we reported earlier based on pooled data from 11 acini (Ender et al., 2022). Most CEs initiate and propagate at the outer layer of the acinus, thus providing a survival signal to those cells. The inner domain exhibited fewer spontaneous ERK pulses and received fewer pulses due to waves that remained predominantly in the outer layer. Such a spatio-temporal distribution of waves creates a lower survival signal in the inner domain, which contributes to a higher apoptosis rate and clearing of the lumen. The formation of hollow acini recapitulates key features of in vivo breast alveoli and makes this system an attractive model to study the effect of MAPK signaling on lumen morphogenesis.

ARCOS is designed to identify spatial clusters and to track them over time, thus making it potentially applicable to diverse spatio-temporal phenomena. The computational cost lies mainly in clustering with dbscan and in the calculation of NNDs between objects in clusters, which is necessary to track CEs. For cases with many objects per cluster, such as calcium waves in Fig. 4 D, the calculation lasts longer and memory requirements increase. However, the analysis of even the longest time lapses presented here lasted several minutes on a contemporary PC. In our opinion, the main hurdle lies in the preparation of data from microscopy images. Modern machine learning–based approaches for image segmentation are robust even in the presence of a fair amount of noise (Fig. S4, A–D and G; and Video 10; Materials and methods). However, quality time series is crucial for the subsequent binarization and CE detection (Fig. S4, E, F, and H; Materials and methods).

Signal binarization may require several trials with either detrending/thresholding or LISA-based methods. In our case, the latter performed better when activity time series had well-defined baselines without long-term trends, as was the case for all datasets except MCF10A, which required detrending due to varying baselines. Furthermore, some understanding of the acquired CEs is necessary to set the search radius ε and the minimum for CE duration and size. If ε is too large, then independent CEs are merged into a single cluster; if too small, then a single CE might be split into smaller events (Fig. S5, A–D; Materials and methods).

A time-lapse with randomly active cells can also lead ARCOS to identify spurious space-time correlations as CEs, especially when the total binarized activity (TBA) is high. To identify this regime, we simulated a 2D epithelium with randomly activating cells for a range of TBA. Above TBA ≈ 10%, due to an increased probability of neighboring cells remaining active, the algorithm identified several CEs comprising at least three cells and lasting at least three frames (Fig. S5 G). Therefore, validation with randomization tests is crucial to assess whether statistics calculated from observed data differ from the statistics obtained in a system with random activations.

Another regime when ARCOS may yield imprecise results is when CEs start overlapping. Simulation of concentrically growing events (Fig. S5 H) revealed that the statistics of collective activations agree with the ground truth up to TBA ≈ 10%, which is the threshold when individual events start colliding (Fig. S5, I–K). Even though ARCOS is designed to handle collisions of CEs, such situations are indiscernible above a certain threshold, making automated detection difficult. The actual percentages above which these problems emerge are only an approximate guideline and may further depend on detection thresholds for CE size and duration and on the nature of the propagating wave. In our experiments, TBA was well below the estimated overlap threshold (0.6% for calcium waves; 2% for bee shimmering, 3 and 7% for ERK activity waves in NRK-52E and MCF10A WT), although it amounted to ≈23% in untreated MCF10A mutant cell lines, which may partly explain the higher duration and size of detected CEs.

We introduced ARCOS, a novel tool for detecting and quantifying collective spatio-temporal phenomena with a focus on cell signaling, and validated it on several biological systems with drastically different spatio-temporal collective behaviors. We analyzed ERK activity waves in an epithelial monolayer that propagate radially on the timescale of 10–30 min, radial calcium waves on the timescale of 10–30 s, directional waves in a shimmering bee colony, and a 3D culture model of mammary acini. The open-source code provided as R and Python packages can be used in Jupyter notebooks or batch processing pipelines. Supplementary R notebooks contain full analysis of all experimental modalities presented here and provide example parameters (Dobrzyński, 2023). The napari plugin enables anyone without extensive programming knowledge to explore parameters through an intuitive graphical user interface on a platform that emerges as a de facto standard for viewing multidimensional images (Video 1).

We unravel a new layer of complexity in cancer signaling by showing that two different oncogenic mutations known to increase ERK activity led to different, distinctive single-cell and collective responses. We focused our experiments on two widely studied and clinically relevant oncogenic mutations: KRAS^G12V and PIK3CA^H1047R. The former directly activates ERK by switching on RAF and MEK (Huang et al., 2021), while the latter activates ERK through a crosstalk mechanism that involves MMP-dependent amphiregulin-mediated EGFR activation (Ender et al., 2022). As observed before (Konishi et al., 2007; Gustin et al., 2009), we show that both mutations lead to increased population-averaged ERK activity in MCF10A epithelial cells. Evaluation of single-cell ERK dynamics reveals two additional levels of complexity—the oncogenic mutations modify: (1) ERK dynamics at the single-cell level; (2) the collective behavior of ERK signaling, as can be quantified using ARCOS.

Regarding ERK dynamics, we find that KRAS^G12V cells display more ERK pulses that also switch off slower than ERK pulses in WT. In marked contrast, PIK3CA^H1047R leads to increased ERK pulse frequency without modifying ERK pulse shape. The KRAS^G12V-induced characteristic ERK dynamics can be explained as follows. Active Ras triggers the tripartite RAF/MEK/ERK network with a negative feedback loop from ERK to RAF (Sturm et al., 2010; Fritsche-Guenther et al., 2011); however, it bypasses the ERK-RSK2-SOS negative feedback loop that acts upstream of the mutated Ras (Dessauges et al., 2022). This partly abrogated negative regulation explains why ERK pulses remain adaptive, albeit with slower dynamics compared with WT, rather than exhibiting sustained ERK activation. The ability of the MAPK pathway to rescale the ERK response in the presence of an activating mutation (Gillies et al., 2017) explains why cells remain sensitive to signaling inputs. In marked contrast, the PIK3CA^H1047R mutation activates the MMP/EGFR/amphiregulin-based crosstalk to activate the MAPK signaling network with intact feedback structures as in WT cells (Ender et al., 2022). This produces ERK pulses of identical shapes as in WT cells, albeit with a higher frequency.

Using ARCOS to quantify population-level ERK dynamics, we report that both KRAS^G12V and PIK3CA^H1047R lead to more, longer, and bigger CEs compared with WT. Because ERK waves rely on MMP-mediated cleavage of EGF ligands to propagate between cells, we used the MMP inhibitor batimastat to inhibit the propagation. As previously quantified, ERK waves vanish in WT cells upon inhibition of cell–cell communication (Gagliardi et al., 2021). We observe the same behavior for PIK3CA^H1047R cells. We explain this by the fact that in this cell line, the crosstalk from PI3K to ERK involves the EGFR-dependent activation by the EGF-like ligand amphiregulin, which also requires MMPs (Young et al., 2015; Ender et al., 2022). Interestingly, the number, length, and duration of CEs are only partly reduced in the KRAS^G12V cells treated with batimastat. We speculate that in addition to increased secretion of amphiregulin in this mutant (Minjgee et al., 2011), mechanical input is also responsible for cell–cell communication. ERK-dependent activation of myosin contractility has been implicated in ERK wave generation during collective epithelial motility, together with the MMP/EGF/EGFR system (Aoki et al., 2017; Hino et al., 2020). The increased periods of ERK activity in KRAS^G12V cells might exacerbate myosin activity and thus increase cell–cell communication independently of the MMP signaling system, leading to residual ERK waves upon MMP inhibition.

Our data demonstrate that the single-cell signaling dynamic does not merely result from network wiring within individual cells but also reflects emergent phenomena at the level of cell community. Disentangling these two contributions is necessary to further interpret oncogenic signaling. We showed that intercellular propagation of dynamic collective signaling plays a major role in the pathological hyperactivation of the MAPK pathway due to oncogenic mutations. Future studies of spatio-temporal dynamics in other oncogenic MAPK/Akt pathway mutants will be essential to better characterize the dependency of individual cells on intercellular communication. Quantification of these complex emergent processes might then allow for modeling signaling networks in cell communities using agent-based approaches, which will enable true understanding of signaling processes at previously inaccessible lengths and time scales.

WT human mammary MCF10A cells were a gift of J.S. Brugge, Harvard Medical School, Boston, MA, USA (Debnath et al., 2003). The isogenic variants of MCF10A harboring the heterozygous mutations PIK3CA^H1047R (Gustin et al., 2009) and KRAS^G12V (Konishi et al., 2007) were a gift of Ben Ho Park, Vanderbilt University Medical Center, TN, USA. NRK-52E expressing EKAREV-NLS were a gift from K. Aoki, National Institute for Basic Biology, Okazaki, Japan (Aoki et al., 2013). Cultivation of WT MCF10A cells and the isogenic derivative was carried out in growth medium composed of DMEM:F12, 5% horse serum, 20 ng/ml recombinant human EGF (Peprotech), 10 μg/ml insulin (Sigma-Aldrich/Merck), 0.5 μg/ml hydrocortisone (Sigma-Aldrich/Merck), 200 U/ml penicillin, and 200 μg/ml streptomycin. NRK-52E cells were cultured in DMEM supplemented with 10% FBS, 200 U/ml penicillin, and 200 mg/ml streptomycin.

To observe nuclei and measure ERK activity, we used the PiggyBac plasmid coding for the stable nuclear marker H2B-miRFP703 (pPBbSr2-H2B-miRFP703) and the PiggyBac plasmids coding for the ERK activity biosensor ERK-KTR fused with mRuby2 (pHygro-PB-ERK-KTR-mRuby2) or mTurquoise2 (pHygro-PB-ERK-KTR-mTurquoise2). To express the CRY2-based optogenetic actuator OptoRAF and the plasma membrane anchor CIBN-KrasCT, we used a PiggyBac plasmid expressing both the proteins separated by the self-cleaving peptide P2A (pPB3.0-PuroCRY2-cRAF-mCitrine-P2A-CIBN-KrasCT). All plasmids were generated for a previous study from our laboratory (Gagliardi et al., 2021). We chose OptoRAF to locally stimulate ERK activity because, in MCF10A cells, the activation of RAF is downstream of the MMPs/EGFR-mediated cell–cell communication that propagates ERK waves. This property makes the OptoRAF system an ideal tool to generate artificial waves by precisely activating only the cells in the illuminated area and not outside.

To generate cells stably expressing the biosensors or the optogenetic actuator, we transfected the PiggyBac plasmids together with the helper plasmid expressing the transposase. The transfection of WT MCF10A cells and the mutated derivatives was carried out with FuGene (Promega) according to the manufacturer’s protocol. Antibiotic selection and image-based screening were used to generate stable clones with the desired biosensors/optogenetic tools.

MCF10A and NRK-52E cells were seeded in 5 mg/ml fibronectin (PanReac AppliChem) on 24-well #1.5 glass bottom plates (Cellvis) to reach confluence on the day of the experiment. The culture medium was replaced with imaging medium (DMEM:F12, 0.5 μg/ml hydrocortisone (Sigma-Aldrich/Merck), 200 U/ml penicillin and 200 μg/ml streptomycin, and 10 mM Hepes) 24 h before the experiments and kept for their entire duration. Image acquisition was executed with an epifluorescence Eclipse Ti inverted fluorescence microscope (Nikon) controlled by NIS-Elements (Nikon) with a Plan Apo air 20× (NA 0.8) or Plan Apo air 40× (NA 0.9) objectives and an Andor Zyla 4.2 plus camera at a 16-bit depth. Laser-based autofocus was used throughout the experiments. All the experiments were performed in a humidified chamber at 37°C and 5% CO₂. Excitation and emission of the fluorescent proteins were obtained with the following monochromatic LEDs and filters. Far red: 640 nm, ET705/72m; red: 555 nm, ET652/60m; NeonGreen: 508 nm, ET605/52; mTurquoise2: 440 nm, HQ480/40. For optogenetic experiments, cells expressing OptoRAF were kept in the dark for at least 24 h before the experiments and all preparatory steps at the microscope were carried out with red or green light (wavelength >550 nm). Stimulation of the OptoRAF was obtained with monochromatic 488 nm blue light for 100 ms at 0.3 W/cm². The stimulated areas were produced by a DMD (Andor Mosaic3) controlled by NIS Jobs based on masks prepared beforehand.

We used three image segmentation approaches. In the first, we used Ilastik (Berg et al., 2019) to classify pixels in the nuclear channel and build a nuclear probability mask. The classifier was trained on manual annotations of nuclei based on the H2B-mRFP703 fluorescent marker for MCF10A cells and EKAREV-NLS for NRK-52E. We used CellProfiler 3.0 (McQuin et al., 2018) to perform further analysis. A threshold-based segmentation identified nuclei from the probability masks. The nuclear segmentation masks were used to quantify the average nuclear pixel intensity in the ERK-KTR channel. To determine the area corresponding to the cytosol, we expanded the nuclei by a predefined number of pixels, but not further than the underlying fluorescence intensity in the ERK-KTR channel. We calculated the average cytosolic pixel intensity of ERK-KTR from the resulting cytosolic ring. By dividing the average pixel intensities from the cytosolic ring and the nucleus, we obtained the cytosolic/nuclear ratio that estimates the ERK kinase activity. For experiments with NRK-52E cells that express EKAREV-NLS biosensor, we calculated the fluorescence resonance energy transfer (FRET) ratio, where the FRET image is divided pixel-by-pixel by the donor image. Tracking of single cells was done on nuclear centroids in MATLAB using μ-track 2.2.1 (Jaqaman et al., 2008).

In the second approach, we used stardist (Schmidt et al., 2018) to segment nuclei, then shrank the nuclei by one pixel, and used this binary image together with the raw ERK-KTR image to segment the cytoplasm of corresponding cells using cellpose (Stringer et al., 2021; Pachitariu and Stringer, 2022). The cytoplasm was shrunk by a predefined number of pixels and mean intensities for nuclei and cytoplasm were extracted from the KTR images to subsequently calculate the KTR ratio as before.

In the third approach, we used a custom segmentation algorithm that extracts image features from the nuclear marker channel using filters from a modified VGG16 (Arsa and Susila, 2019) convolutional neural network pretrained on ImageNet (Deng et al., 2009). We then quantified nuclear and cytosolic fluorescence intensity following the procedure described above. To track nuclei, we used a Python implementation of the linking algorithm (Crocker and Grier, 1996) provided in the trackpy library (Allan et al., 2021).

Apoptotic events were manually annotated based on nuclear shrinkage, as reported earlier (Gagliardi et al., 2021).

Live imaging of 3D mammary acini with MCF10A cells and the corresponding image analysis were executed as described earlier (Ender et al., 2022). MCF10A cells were dissociated to single cells and embedded in growth factor–reduced Matrigel (Corning) and overlaid with DMEM/F12 supplemented with 2% horse serum, 20 ng/ml recombinant human EGF, 0.5 mg/ml hydrocortisone, 10 mg/ml insulin, 200 U/ml penicillin, and 200 mg/ml streptomycin. After 3 d, EGF, insulin, and horse serum were removed from the medium. 25 mM Hepes was added to the medium before imaging.

Image acquisition was performed with an epifluorescence Eclipse Ti2 inverted fluorescence microscope (Nikon) equipped with a CSU-W1 spinning-disk confocal system (Yokogawa) and a Plan Apo VC 60× water immersion objective (NA = 1.2). Images were acquired with a Prime 95B camera at 16-bit depth (Teledyne Photometrics). 3D segmentation of nuclei and extraction of cytosolic/nuclear ERK-KTR fluorescence intensity ratio were performed using a customized version of the LEVER software (Winter et al., 2016).

Image acquisition was performed as described earlier (Kastberger et al., 2011). To segment the act of flipping the abdomen in individual bees, we used a pixel classifier from the Ilastik image processing software (Berg et al., 2019). We obtained a probability mask that could be thresholded and processed as a conventional single-cell time-lapse movie that resulted in single-bee time series data. Motion-active bees, i.e., bees that flipped their abdomen, were identified with the local G* statistic, which we then thresholded and fed into ARCOS.

Single-cell ERK activity obtained from image segmentation yields a (x,y,m,t) spatio-temporal point pattern, where (x,y) are centroid coordinates of nuclei, m is the corresponding measurement value at these coordinates (e.g., ERK activity), and t denotes time. To find CEs, we identify “active cells,” i.e., cells with m above a certain threshold. We implemented three methods in the binMeas function of the ARCOS R and Python packages to identify activity periods in time series. In the simplest case (the biasMet parameter set to “none”), the measurement is rescaled globally to [0,1], smoothed with a short-term median filter to remove fast frequencies and/or incidental peaks, and binarized with a global threshold. Setting biasMet = “runmed” activates detrending, where a long-term median-smoothed measurement is subtracted from the initial time series. With biasMet = “lm,” the detrending step subtracts a linear fit to the measurement. The aim of this step is to remove long-term trends due to photobleaching of fluorescent probes or to remove permanently active cells. The subsequent thresholding of the detrended signal may lead to identification of longer pulses only during their peak activity, thus omitting the leading and the decaying phase of the pulse (Fig. 1 D).

For the LISA-based activity binarization, we used the statistics available in the R package elsa (Naimi et al., 2019). These include local Moran’s I and local Geary’s c (Anselin, 1995), as well as G and G* statistics (Getis and Ord, 1992). The resulting spatial correlation coefficient is binarized with a fixed threshold, as exemplified in Fig. 1 E.

The choice of the binarization algorithm depends on the dataset and should be evaluated on a case-by-case basis. The main disadvantage of detrending with the running median is that it requires long time series to apply the median filter to estimate the baseline. The recommended window for such a filter usually amounts to a quarter of the time series’ length. This condition may be difficult to satisfy when image segmentation errors produce short single-cell time series. Instead, LISA-based binarization is not contingent on long time series because frames are analyzed independently. However, the difficulty lies in the choice of a LISA statistic, e.g., local Moran’s I, Geary’s c, G, and G*, for each of them measures a different aspect of the spatial pattern. For example, Moran’s I measures the local deviation from the global mean, local G measures the clustering of hot/cold spots, and Geary’s c measures the degree of dissimilarity of a value to its neighbor (Naimi et al., 2019). Moreover, most LISAs are not numerically bounded and are sensitive to the level of global autocorrelation, thus making them difficult to compare between systems with different overall clustering.

After time series binarization, the space–time point pattern (x,y,t) of active cells is passed to the clustering and tracking algorithm (Fig. 1 B). In the first step, a dbscan algorithm (Ester et al., 1996; Campello et al., 2013; Hahsler et al., 2019) is applied independently to every time frame t of the time-lapse to spatially cluster active cells. In the ARCOS R package, we use the implementation from the dbscan package and the Scikit-learn implementation in the Python version. Two, user-provided parameters determine this clustering: ε, the radius of the neighborhood, and minClSz, the minimum number of cells that form a cluster. The radius ε must be large enough to form clusters of neighboring active cells. We recommend setting it slightly larger than the mean first nearest neighbor distance calculated for all cells in the entire time-lapse. When ε is too small, cells are grouped into small clusters instead of one. When ε is too big, multiple CEs are joined. The minClSz parameter corresponds to the minimum number of cells that are considered a CE. Setting it to 1 allows detecting the initial stage of a CE, while risking that a sudden onset of simultaneously activating cells may be assigned to separate clusters. A bigger minClSz lets the algorithm delay the clustering until a larger event develops. A sudden appearance of multiple active cells may indicate an insufficient acquisition frequency during the experiment.

In the second step, identified clusters are tracked between the frames. The number of previous frames to link objects in the current frame depends on the parameter nPrev. Usually, clusters between two consecutive frames are linked (nPrev = 1) but the number might be increased if missing frames are present in the dataset. Assuming, for simplicity, that nPrev = 1, the cluster tracking is performed as follows. For each cluster in t, we calculate the first NND between all cells in the current cluster and all cells in all clusters in the frame t−1. We use the nn2 function from the RANN package in R or the Scipy implementation of the KDTree nearest neighbor search in Python to calculate the distances. If at least a single cell in the current cluster is within a threshold distance ε_prev to a cell in the previous frame, cells in the current cluster inherit cluster IDs from first nearest neighbor cells in the previous frame. The ε_prev parameter is the threshold for detecting cluster displacements between the frames. Detection of a propagating activation wave as several consecutive CEs suggests that ε_prev should be increased.

To display the relationship between the size and duration of CEs, we used hex-bin plots (Fig. 3 E; and Fig. 4, C and F), as they avoid overplotting that can occur in a conventional scatter plot when dealing with large datasets. The hexagonal binning tiles are colored according to the number of points in those bins, and we used a logarithmic color scale.

We propose five randomization implementations in ARCOS R/Python packages to simulate the null distribution of CEs’ statistics (Fig. S1). The methods vary regarding assumptions about the null hypothesis. The first shuffles the entire time series in space between initial spatial positions of existing cells (Fig. S1 B). This is the recommended method as it only disturbs the spatial component while conserving individual cells’ measurement dynamics and the population-averaged activity over time. However, since this approach moves the entire time series to a new location, they may end up outside the initial FOV. This may happen when cells migrate, and the new initial position after shuffling is close to the FOV’s border. If that is the case, we recommend other approaches: randomly shifting the measurement for every time series back and forth in time (Fig. S1 C), shuffling whole activity blocks per time series (Fig. S1 D), or shuffling individual time points along every time series (Fig. S1 E). All three methods preserve the total activity in individual time series and leave them in their original positions but affect the onset (randomly shifting the time series) or randomize the activity dynamics (shuffling blocks or time points of activity), thus potentially disrupting correlations between the neighbors. If these methods are insufficient to rearrange the spatio-temporal pattern, e.g., for a single large collective activation, they can be combined with the last approach, where measurements are shuffled between the cells independently for every frame (Fig. S1 F). Note that such complete randomization changes the tested question, for the observed pattern is compared to a uniform distribution of activations in space and time. It is likely that the observed statistics will be sufficiently away from such a null distribution to reject the null hypothesis at a given significance level but will not be extreme enough when the time series are only shuffled in space with their dynamics intact.

The choice of the randomization method depends on the type of locality in the observed phenomenon. If cell activity is local in space but uniform in time, e.g., sustained or synchronized oscillatory activation of a cell cluster, then shuffling positions of entire trajectories or shuffling positions per time point would disrupt the spatial correlation. However, when the observed phenomenon is local in time but uniform in space, e.g., a synchronized, transient activation pulse of all cells in the FOV, position shuffling will not affect such an activity cluster. Instead, shifting trajectories in time, shuffling time points, or blocks of activity per trajectory would be required to disrupt the temporal locality.

The implementation of the ARCOS package in R offers two functions to generate a single or multiple CEs on a 2D lattice over time. genSynthSingle2D creates a single event with predetermined cell activations over time. genSynthMultiple2D generates a random number of concentrically growing circles with a random duration at random points in space and time. Additionally, we used NetLogo, a multiagent programmable modeling environment (Wilensky, 1999), to simulate coordinated ERK activation in a 2D epithelium (Figs. S1 and S5). Configuration files to reproduce the simulations are in supplementary data (Dobrzyński, 2023).

In Fig. 2, E, G, and I; Fig. 3 H; and Fig. 4, B, E, and H, a global Moran’s I autocorrelation coefficient was calculated according to the method described in Gittleman and Kot (1990). We applied the Moran.I function from the R package ape (Paradis and Schliep, 2019) to the (x,y,m) point pattern and calculated the coefficient independently for every time frame t. From the observed Moran’s I, we subtracted the expected value of I under the null hypothesis that there was no correlation. Quantification with the global correlation coefficient is useful only for cases when a single (or a few) CE occurs in a FOV.

for a given SNR iteration and averaged over 10 repeats.

To test the sensitivity of CE detection on image quality, we used a time-lapse of an ERK activity wave in MDCK epithelial cells stably expressing the EKAREV-NLS FRET biosensor (Fig. S4 A and Video 10). We report data degradation using SNR expressed in decibels or dB (Eq. 2). A SNR > 0 dB indicates that the signal power, i.e., the average of squares of all pixel intensities over time, is greater than the power of added noise; SNR < 0 dB is when the power of noise dominates the power of the original data. The higher the ratio, the cleaner the signal, with +∞ corresponding to original, measured images.

For a chosen set of segmentation parameters in stardist (Schmidt et al., 2018), we observed that nuclei segmentation was severely affected only when SNR < −10 dB (Fig. S4 B). Below that SNR, higher noise in images resulted in severely deteriorated time series (Fig. S4 C) to a point when ARCOS could no longer detect CEs (Fig. S4 D). The number of false negatives, i.e., cells that were not assigned to true CEs increased sharply for SNR < −10 dB, which corresponds to noise power being 10× the power of the original signal. The number of false positives, i.e., cells that were erroneously assigned to clusters, increased only slightly at SNR < −10 dB (Fig. S4 G). Overall, the image segmentation and cell tracking pipeline can tolerate noisy images for SNR > −10 dB before recall (Eq. 4) or the number of relevant items retrieved decreases (Fig. S4 G). The code to reproduce this analysis is available in the accompanying data repository (Dobrzyński, 2023).

To study the effect of noise on ARCOS analysis, we added Gaussian white noise with an increasing amplitude to single-cell ERK activity from the segmentation of the wave in MDCK cells (Fig. S4, E and F). ARCOS was then run on the newly generated dataset and compared with the original output as above. This was carried out for a range of SNRs and repeated 10 times.

CEs were faithfully detected only down to SNR ≈ 20 dB (Fig. S4 F), which corresponds to the power of the original time series being 100× higher than the power of the added Gaussian noise. For lower SNR, the number of false negatives soared, which caused the recall to drop sharply (Fig. S4 H). We conclude that for the parameters in our pipeline, CE detection is more sensitive to degradation of time series rather than images.

The core ARCOS algorithm clusters active cells in individual frames of the time-lapse and then tracks the clusters. The key parameter for spatial clustering is the neighborhood radius, ε. It should be greater than the first NND between active cells to let the dbscan algorithm form clusters. If too large, nearby independent CEs are included in a single cluster; too small, and a single CE is identified as several smaller events. We investigated the quality of event detection as a function of ε by simulating a synthetic time-lapse of 200 frames in NetLogo. We simulated 69 radially expanding, collective activations randomly spaced on a 2D hexagonal mesh, like ERK waves observed in the MCF10A epithelium (Fig. S5 A). The simulation established a ground truth with a known number (69 events), size (37 unique cells), and duration (9 time frames) of CEs. The events were identified correctly for a range of $ε∈[1.5,4]×NND$ (Fig. S5, B–F), bar several small errors due to overlaps between simulated events. For ε < 1.5, cells were not clustered together and formed their own clusters; for ε > 4, the algorithm identified large clusters that comprised many independent events. These results hint at a general rule for processing spatial datasets: ε ought to be slightly larger than the NND, but care must be taken to set it low enough to avoid bundling nearby events.

Fig. S1 presents various randomization techniques used to evaluate the accuracy of the ARCOS output. Fig. S2 depicts statistics derived from randomization tests of collective events detected in wild-type, KRAS^G12V, and PIK3CA^H1047R cell lines. Fig. S3 depicts the same statistics as in Fig. S2 but applied to batimastat treatment on the same cell lines. Fig. S4 explores the sensitivity of ARCOS to image or time-series degradation. Fig. S5 portrays the algorithm’s sensitivity to the search radius ε and the total binarized activity. Video 1 shows the features implemented in the ARCOS napari plugin. Video 2 displays the detection of a growing circular region generated in cells with the OptoRAF/ERK-KTR actuator/biosensor circuit. Video 3 displays a moving circle generated in cells with the same technique. Video 4 displays a splitting circular region generated in cells with the same technique. Video 5 shows ARCOS analysis of collective ERK activity waves in MCF10A WT, KRAS^G12V, and PIK3CA^H1047R epithelial cells. Video 6 shows collective ERK activity waves in starved NRK-52E epithelial cells detected with ARCOS. Video 7 displays collective calcium waves identified by ARCOS in starved MDCK epithelial cells expressing the calcium reporter GCaMP6S. Video 8 demonstrates the ability of ARCOS to analyze shimmering waves in a giant honeybee colony. Video 9 showcases collective 3D ERK activity waves, identified by ARCOS in a mammary acinus in vitro. Video 10 demonstrates ARCOS’s sensitivity to iterative image degradation with additive Gaussian white noise.

The data, R notebooks to reproduce the plots and randomization tests, configurations of simulations, and Jupyter notebooks to show the results in the napari image viewer were uploaded to Mendeley Data (Dobrzyński, 2023). ARCOS can be downloaded from https://arcos.gitbook.io as an R or a python package or a napari plugin together with extensive documentation.

We would like to thank J.S. Brugge (Department of Cell Biology, Harvard Medical School) for MCF10A WT cells, B.H. Park (The Vanderbilt-Ingram Cancer Center, Division of Hematology/Oncology, Department of Medicine, Vanderbilt University Medical Center) for the isogenic knock-in MCF10A line with PIK3CA H1047R and KRAS G12V mutations, K. Aoki (Division of Quantitative Biology, National Institute for Basic Biology, Okazaki, Japan) for NRK-52E cells expressing EKAREV-NLS, Y. Fujita (Division of Molecular Oncology, Institute for Genetic Medicine, Hokkaido University Graduate School of Chemical Sciences and Engineering, Sapporo, Japan) and Y. Takeuchi (Division of Cancer Cell Biology, Cancer Research Institute, Kanazawa Univerisity, Japan) for providing MDCK calcium imaging data, T. Höhener (Institute of Cell Biology, University of Bern, Bern, Switzerland) for discussions, and the reviewers for the careful and insightful assessment of our manuscript.

This work was supported by Swiss National Science Foundation grants 31003A-163061, 51PHPO-163583, Div3 310030_185376, IZKSZ3_62195, and by a Swiss Cancer League grant KLS-4867-08-2019 to O. Pertz, by a Human Frontier Science Program grant RGP0043/2019 to O. Pertz and A. Cohen, and by a Chan Zuckerberg Initiative napari Plugin Foundation Grant 2022-252527 to M. Dobrzyński and B. Grädel.

Author contributions: M. Dobrzyński and P.A. Gagliardi conceived and formulated the project. M Dobrzyński and M.-A. Jacques developed the R package. B. Grädel developed the Python package and the napari plugin. A. Cohen and L. Hinderling developed image analysis software. P.A. Gagliardi and P. Ender designed and performed the experiments and acquired images of 2D and 3D epithelia. G. Kastberger designed and performed the experiments and acquired images of giant honeybees. A. Cohen, P.A. Gagliardi, and L. Hinderling performed image analysis. M. Dobrzyński and B. Grädel analyzed and visualized the data. P.A. Gagliardi designed and performed NetLogo simulations. A. Cohen, M. Dobrzyński, B. Grädel, and O. Pertz secured the funding. M. Dobrzyński, P.A. Gagliardi, and O. Pertz wrote the manuscript, and all authors discussed the results and participated in revising the manuscript.