CausalXtract: a flexible pipeline to extract causal effects from live-cell time-lapse imaging data

and Discussion CausalXtract feature extraction and causal discovery modules CausalXtract’s live-cell image feature extraction module (Cell- Hunter+), Fig. 1b, is based on CellHunter software 1 and con- sists in three steps: detection, tracking and feature extraction of live cells within time-lapse video images. First, automatic lo- calization/segmentation of cells (e.g. tumor and immune cells) is performed with the Circular Hough Transform (CHT) algo- rithm2 to estimate the cell centers and radii. Second, cell tra- jectories along the frames are constructed by linking the posi- tions detected at the previous time step through Munkres’ algo- rithm for Optimal sub-pattern Assignment Problems (OAPs) 3. Finally, relevant descriptors related to the shape, motility, and state of the cells, as well as cell-cell interactions are quantified from each cell trajectory (Methods). CausalXtract’s temporal causal discovery module (tMIIC), Fig. 1c, is adapted from the causal discovery method, MIIC 4–6, which learns contemporaneous causal networks (i.e. when tem- poral information is not available) for a broad range of bio- logical or biomedical data, from single-cell transcriptomic and genomic alteration data 4,7 to medical records of patients 5,6,8. Live-cell time-lapse imaging data contain, however, informa- tion about cellular dynamics, which can in principle facilitate the discovery of novel cause-effect functional processes, based on the assumption that future events cannot cause past ones. To this end, CausalXtract’s discovery module, tMIIC, recon- structs time-unfolded causal networks, where each variable is represented by several nodes at different relative time points 9, Fig. 1c. Such a time-unfolded network framework 10–13 is re- quired to account for the temporal correlation between succes- sive time steps in time series data. We benchmarked tMIIC on synthetic datasets resembling the real-world data of interest an- alyzed in this study (i.e. number of time steps, network size and degree distribution) and found that it matches or outperforms state-of-the-art methods, while running order of magnitudes faster on datasets of biologically relevant size including tens to hundreds of thousands time steps, Supplementary Figs. 1-4. CausalXtract’s temporal network framework goes beyond the seminal concept of temporal causality originally proposed by Granger 14 for linear time series without reference to graph- ical models and later extended to non-linear dynamics by Schreiber15,16. In particular, Granger-Schreiber causality is in fact too restrictive and may overlook actual causal effects, that can be uncovered by graph-based causal discovery meth- 1 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 8, 2024. ; https://doi.org/10.1101/2024.02.06.579177doi: bioRxiv preprint c Tumor−on−chip preparationa b CausalXtract’s temporal causal discovery module (tMIIC) CausalXtract’s live−cell image feature extraction module (CellHunter+) Fibroblast (CAF) cells Her2+ cancer cells Her2+ cancer cells Immune primary cells µ500 m Endothelial primary cells Figure 1: CausalXtract pipeline. a , Live-cell tumor ecosystem reconstituted ex vivo 1 using the tumor-on-chip technology (Methods). b, CausalXtract’s live-cell image feature extraction module (CellHunter+). The tracking of cancer and immune cells and of their mutual inter- actions is illustrated in Supplementary Movies 1-3, in absence or presence of cell division and apoptosis event. c, CausalXtract’s temporal causal discovery module (tMIIC) learns a temporal causal network from the features extracted in (b). See Methods for CausalXtract’s implementation details and theoretical foundations. A step-by-step notebook of CausalXtract pipeline is provided with the source code. 2 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 8, 2024. ; https://doi.org/10.1101/2024.02.06.579177doi: bioRxiv preprint ods, Supplementary Fig. 5 (Methods, Theorem 1). In addition, Granger-Schreiber causality has long been known to infer spu- rious causal associations based on time delays, by excluding the presence of latent common causes a priori 9. CausalXtract circumvents these limitations by combining graph-based and information-based approaches (Methods), while including con- temporary and time-delayed effects of unobserved latent vari- ables, that are ubiquitous in cell biology data ( e.g. the latent effects of cell cycle phases on cellular features and responses). Application to tumor-on-chip cellular ecosystems We showcase CausalXtract with the analysis of time-lapse im- ages of a tumor ecosystem reconstituted ex vivo using the tumor-on-chip technology, Fig. 1a. These live-cell time-lapse images come from a proof-of-concept study 1 which demon- strated the effects of an anti-cancer drug (the monoclonal an- tibodies trastuzumab, brand name Herceptin, used to treat HER2+ breast cancers) on a reconstituted tumor microenvi- ronment including cancer cells, immune cells, cancer-associated fibroblasts (CAF), and endothelial cells (Methods). However, a comprehensive extraction and analysis of cellular morphody- namic features and interactions remained unexplored. To this end, cellular features such as cell geometry, velocity, division, apoptosis, cell-cell transient interactions and persis- tent contacts were first extracted from the raw images using CausalXtract’s feature extraction module, Fig. 1b and Supple- mentary Fig. 6. Then, a time-unfolded causal network, Sup- plementary Fig. 7, and the corresponding summary causal net- work, Fig. 2a, were reconstructed between extracted cellular features, cell-cell interactions and therapeutic conditions using CausalXtract’s temporal causal discovery module, Fig. 1c. CausalXtract inferred network, Fig. 2a, uncovers novel bi- ologically relevant findings, in addition to confirming known

from earlier studies. In particular, CausalXtract discov- ers that CAFs directly inhibit cancer cell apoptosis, indepen- dently from anti-cancer treatment, Fig. 2b, while earlier stud- ies reported that CAFs merely reduced the effect of treatment 1. CausalXtract also discovers that treatment increases cancer cell perimeter, Fig. 2c, which has not been reported so far either. In addition, CausalXtract confirms known results from ear- lier studies. In particular, it recovers that treatment increases cancer cell apoptosis and the number of cancer-immune inter- actions, as well as decreases the division rate of cancer cells, Fig. 2c. Likewise, CausalXtract recovers that CAFs stimulate cancer cell migration and increase their area, Fig. 2b. Interestingly, CausalXtract identifies also multiple and possi- bly antagonistic effects with different time delays. For instance, CausalXtract recovers several antagonistic relations between morphodynamic features such as cell division and eccentricity or cell apoptosis and area, Fig. 2d. Indeed, the late phases of cell division are associated to a marked increase in eccentricity (red edge) but preceded by a net decrease in eccentricity, two to three hours before cytokinesis (blue edges), once the decision to divide has been made ( i.e. the probable latent cause) and the cell is actually duplicating its biological materials (prophase), Fig. 2d. Likewise, the area change upon apoptosis is predicted to first decrease soon after apoptosis (blue edge) before eventu- ally increasing upon cell lysis (red edge), Fig. 2d. These results are robust to variations in sampling rate, Supplementary Fig. 8. All in all, CausalXtract is a flexible pipeline which uncovers novel and possibly time-lagged causal relations between cellular features under controlled conditions ( e.g. drug). CausalXtract uniquely combines live-cell feature extraction with information theory and causal discovery approaches. It consists of two in- dependent computational modules, conceived to warrant inter- operability with alternative live-cell segmentation and tracking

or alternative temporal causal discovery methods. CausalXtract opens up new avenues to analyze live-cell imaging data for a range a fundamental and translational research applications, such as the use of tumor-on-chips to screen immunotherapy responses on patient-derived tumor samples. With the advent of virtually unlimited live-cell image data, flexible hypothesis-free interpretation methods are much needed17 and we believe that CausalXtract can bring unique in- sights based on causal discovery to interpret such information- rich live-cell imaging data.

Tumor-on-chip preparation and live-cell microscopy Videos analyzed in the present study refer to biological experiments emulating a 3D breast tumor ecosystem 1. All tumor-on-chip experi- ments have a central endothelium compartment containing endothe- lial cells (primary human umbilical vein endothelial cells, HUVECs) and two lateral chambers filled with biomimetic hydrogel (collagen type I at 2.3 mg/mL) seeded with cancer cells (HER2+ breast can- cer BT474 cell line) and immune cells (peripheral blood mononu- clear cells, PBMCs) from healthy donors, Fig. 1a. Four experimen- tal conditions were considered depending on the presence or absence of breast cancer-associated fibroblasts (CAF cell line Hs578T) and drug treatment (trastuzumab, Herceptin). Videos were acquired by inverted motorized Leica microscopes with a frame rate of 2 minutes for up to 48h (1440 frames). Fig. 1b shows a crop frame with can- cer cells, PBMCs and CAFs. Each video was cropped into multiple small 300×300 pixel videos (referred to as crops in the following), each of which represented a field of view at subsequent time frames containing a “main” cancer cell (MCC) initially placed at the cen- ter of the image, some PBMC immune cells, other cancer cells and possibly CAFs within the surrounding of the MCC depending on the experimental conditions. 36 video crops of up to 1440 frames were analyzed (46,935 frames in total) corresponding to 9 video crops per experimental conditions. CausalXtract’s live-cell image feature extraction module The live-cell image feature extraction module (CellHunter+), Fig. 1b, extends the CellHunter software 1 and consists in three steps: detection, tracking and feature extraction of live cells within time- lapse video images. First, cell detection is based on the segmentation of circular-shaped objects using CHT 2 with radii set around the the- oretical radii of the two cell populations ( rim = 4 px for immune cells and rca = 14 px for MCCs with a pixel resolution 1 px = 0.645 µm1). Then, cell tracking is performed by linking cells detected at the ith frame to cells located at the (i + 1)th frame within a maximum dis- tance from the detected cell candidate. While the motions of both MCCs and immune cells ressemble random walks with time-varying drift and volatility, these two cell types exhibit different motility characteristics1. Hence, different maximum distances are considered for the two cell populations: it was set to 40 px for MCCs and to 20 px for immune cells. For each cell population, an OAP using the Munkres algorithm 3 is solved: the globally best possible pairing among located objects is based on an assignment cost equal to the inverse of the distance between pairs of cell candidates at the ith and (i + 1)th frames. Cell appearing/disappearing and cell overlaps due to projection errors of the 3D scene in the 2D domain are also handled. Finally, cellular morphodynamic features and cell-cell in- teraction features are extracted at successive positions along each 3 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 8, 2024. ; https://doi.org/10.1101/2024.02.06.579177doi: bioRxiv preprint dc ba division 500 time steps 0 1000 eccentricity area apoptosis Figure 2: Application of CausalXtract to time-lapse images of tumor ecosystems reconstituted ex vivo 1. a, Summary causal network inferred by CausalXtract. The underlying time-unfolded causal network is shown on Supplementary Fig. 7. Red (resp. blue) edges correspond to positive (resp. negative) associations. Bidirected dashed edges represent the effect of unobserved (latent) common causes. Annotations on edges correspond to time delays in time-steps (1 ts = 2 min). The inferred network is largely robust to variations in sampling rate (δτ ) and maximum lag (τ ), Supplementary Fig. 8. Here δτ = 7 ts and τ = 84 ts are chosen automatically by CausalXtract, Supplementary Fig. 8b. b, The CAF presence subnetwork highlighting the direct causal effects of CAFs on cancer cells. In particular, CausalXtract uncovers that CAFs directly inhibit cancer cell apoptosis independently from treatment, which has not been reported so far. c, The treatment subnetwork highlighting the direct causal effects of treament on cancer cells. In particular, CausalXtract uncovers that treatment increases cancer cell perimeter, which has not been reported either. d, The eccentricity-area subnetwork highlighting multiple direct and possibly antagonistic time-lagged effects, notably, between cell division and eccentricity and between cell apoptosis and area, as discussed in main text. 4 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 8, 2024. ; https://doi.org/10.1101/2024.02.06.579177doi: bioRxiv preprint trajectory. For each MCC, 15 descriptors were extracted, Supple- mentary Fig. 6, and classified into four main categories: cell shape, motility, state, and interaction descriptors. Shape descriptors. The active contour algorithm implemented in Matlab18 was used to segment the MCC boundaries on each video crop frame. Taking as input a frame representing the ith snapshot of the tth MCC, it returns a binary image, where the MCC is rep- resented by a white region. From the binary image, the shape prop- erties of the region occupied by each MCC were extracted using the Matlab regionprops algorithm. The resulting descriptors of the ex- tracted shape are listed below: • area indicates the number of pixels composing the region. The equivalent diameter of the tth MCC in the ith frame is defined as dt i = √ 4 · area/π. • perimeter represents the distance along the MCC boundary. • circularity is defined as 4 · area · π/perimeter 2, which is equal to 1 when the region is perfectly circular. • eccentricity denotes the eccentricity of the ellipse with the same second moments as the region. The value is equal to 1 when the region is a line and to 0 when the region is a circle. • instantaneous shape change is defined as, |dt i − dt i−1|, corre- sponding to the difference in absolute value of the equivalent diameters between the ith and (i − 1)th frames of the tth MCC. Motility descriptors. The positions pt i = ( xt i, yt i ) and pt i−1 of the tth MCC in the ith and (i − 1)th frames were compared using the Euclidean distance d(·) to define the following motility parameters: • instantaneous cancer velocity 19 is defined as d(pt i, pt i−1)/∆t, where ∆t is the time interval between two consecutive frames. • net displacement 19 indicates the resultant distance between the initial and current positions of the tth MCC, d(pt 1, pt i). • directionality19 is defined as the ratio of net displacement, d(pt 1, pt i), and curvilinear distance, ∑i k=2 d(pt k, pt k−1). It mea- sures the persistence of motion and ranges from 0 for confined cells to 1 for cells moving perfectly straight in one direction. State descriptors. They record apoptosis or division events: • apoptosis indicates if the MCC has died during the experiment. It is set to ‘No’ as long as the cell has not died and becomes ‘Yes’ for the remaining frames after the cell undergoes apoptosis. • division indicates if the MCC has divided during the experi- ment. It is set to ‘No’ as long as the cell has not divided and becomes ‘Yes’ for the remaining frames after the cell divides. Interaction descriptors . Interactions between MCCs and immune cells were defined with respect to two radii around each MCC, r1 = rim + rca + 2 = 20 px and r2 = 2 × (rim + rca) = 36 px1. Hence, r1 refers to MCC and immune cells in actual physical contact, while r2 refers to MCC and immune cells in close vicinity. Then, for each sample the following interaction features were defined: • number of cancer-immune interactions (r 2) corresponds to the number of immune cells within the interaction radius r2 around the MCC on that frame. • number of cancer-immune interactions (r 1) corresponds to the number of immune cells in close contact with the MCC on that frame. • minimal cancer-immune distance (r 2) is the minimum distance between the MCC and the immune cells within a radius r2. • mean immune velocity (r 2) is the mean instantaneous veloc- ity norm of the immune cells within the interaction radius r2 around the MCC. • mean immune velocity ( r1) is the mean instantaneous velocity norm of the immune cells in close contact with the MCC. Overview of causal discovery methods for non-temporal data Traditional causal discovery methods 20,21 aim to learn causal net- works from datasets of independent samples by proceeding through successive steps. They first learn structural constraints in the form of unconditional or conditional independence between variables and remove the corresponding edges from an initial fully connected net- work. The second step then consists in orienting some of the re- tained edges based on the signature of causality in observational data. This corresponds to orienting three-variable “v-structure” mo- tifs as, X → Z ← Y , whenever the edge X − Y has been removed without conditioning on the variable Z, which implies that Z can- not be a cause of X nor Y . This does not guarantee, however, that X (or Y ) is an actual cause of Z, which also requires to rule out the possibility that the edge between X and Z (or Y and Z) might originate from a latent common cause, L, unobserved in the dataset, i.e. X L99 L 99K Z. In addition, classical causal discov- ery methods are prone to spurious conditional independences, which lead to many false negative edges and limit the accuracy of inferred orientations. The recent causal discovery method, MIIC 4–6, which combines constraint-based and information-based principles, learns more robust causal graphical models by first collecting iteratively sig- nificant information contributors before assessing conditional inde- pendences. In practice, MIIC’s strategy limits spurious conditional independences which improves its edge sensitivity and orientation reliability compared to traditional constraint-based methods 4–6. In addition, MIIC can handle missing data 5 and also heterogeneous multimodal data, by analyzing continuous and categorical variables on the same footing, based on a mutual information supremum prin- ciple for finite dataset 5,6. Last, MIIC distinguishes genuine causal relations from putative and latent causal effects 6, that are ubiquitous in real-world applications. CausalXtract’s causal discovery module for time series data In order to analyze time series datasets, CausalXtract’s causal dis- covery module (tMIIC) aims to learn a time-unfolded graph, Gt, where each variable is represented by a series of nodes associated to its value at different relative time points, Fig. 1c. Such a time- unfolded network framework 10–13 is required to account for the tem- poral correlation between successive samples in time series data. As- suming that the dynamics can be considered stationary (see Bench- marking of CausalXtract’s causal discovery module section, below), the time-unfolded graph, Gt, should be translationally invariant over time and can be assigned a periodic structure a priori. In addition, Gt can be restricted to a few time steps from the running time, t, back to a maximum time lag, t − τ , since nodes at future time points (t′ > t) cannot a priori influence the observed data at current or previous time points (t ′ ⩽ t), Fig. 1c. The maximum time lag τ should be chosen so as to have little effect on the final graphical model, which can be achieved for instance by setting τ to twice the average relaxation time of the variables of the dataset. In practice, we may also limit the number of time points ν in Gt by introducing a time increment δτ between consecutive time points, which leads to ν = τ /δτ time-lagged layers in Gt. Such a compact periodic graphical representation over a sliding temporal window is learned with tMIIC, which extends MIIC causal discovery method to analyze time series data. First, tMIIC identifies all necessary edges involving at least one contemporaneous node at time t, Fig. 1c. Once these time-lagged and contemporaneous nec- essary edges have been identified, they are simply duplicated at ear- lier time points to enforce the translational invariance of Gt skeleton. Time-lagged edges are then pre-oriented with a first arrowhead point- ing towards the future, considering that current time points cannot cause earlier events. Then, contemporaneous and time-lagged edges can be further oriented using MIIC orientation probability scores ap- plied to Gt, which may also uncover a second arrowhead (backward in time) for time-lagged edges. This corresponds to time-lagged latent causal effects from unobserved common causes, Fig. 1c. Learning such structural models including latent variables from 5 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 8, 2024. ; https://doi.org/10.1101/2024.02.06.579177doi: bioRxiv preprint time series data was first proposed for time-lagged effects 10 and sub- sequently extended to contemporaneous effects 11 by adapting the constraint-based FCI method allowing for latent variables 21. While traditional constraint-based methods suffer from poor recall, the re- cent PCMCI 12 / PCMCI+ 22 method improves recall by introducing ad hoc conditioning rules for auto-correlated time series. By con- trast, tMIIC does not require any ad hoc conditioning rules, as it relies on the same robust information-theoretic strategy as MIIC to limit spurious independence and improve edge recall. tMIIC also captures time-lagged and contemporaneous effects due to latent vari- ables. Relation to Granger-Schreiber temporal causality The concept of temporal causality was originally formulated by Granger14 without reference to any graphical model by compar- ing linear autoregression with or without past values of possible causal variables. This was later extended to non-linear relations by Schreiber15,16 using the notion of Transfer Entropy, TX→Y , which can be expressed in terms of multivariate conditional information, TX→Y = I(Yt; Xt′<t|Yt′<t) (1) where Xt′<t and Yt′ T Y →X ⩾ 0) does not necessarily translate into causal direction as this asymmetry is also expected for non-causal relations. Interestingly, this is in fact the absence of Transfer Entropy in one direction ( e.g. TZ→X ≈ 0) which suggests the possibility of a causal relation in the opposite direction, X → Z, as in the case of v-structures in graph-based causal discovery methods, provided that a latent common cause can be excluded between the two variables (as discussed above). We clarify in Theorem 1 below this relation between temporal causality without reference to any structural model (Eq. 1) and structural causality entailed by time-unfolded causal graphical mod- els ( Gt). This highlights the common foundations of temporal and structural causalities beyond their seemingly unrelated definitions. Theorem 1. [TY →X = 0 implies temporal (2 var + t) v-structures] If Xt is adjacent to Yt in Gt and TY →X = I(Xt; Yt′<t|Xt′<t) = 0 , then for all Yt′ adjacent to Yt in Gt, with t′ < t, there is a temporal (2 var +t) v-structure, Yt′ → Yt ← Xt, in Gt, Supplementary Fig. 5a. Proof : if TY →X = I(Xt; Yt′<t|Xt′<t) = 0 , then all pairs (Xt, Yt′ ) should be unconnected (assuming ‘faithfulness’, i.e. no coincidental cancellation of effects) and all unshielded triples Yt′ Yt Xt should be temporal v-structures, Yt′ → Yt ← Xt, as Yt /∈ Xt′<t in TY →X □ Note, however, that the converse of Theorem 1 is not true: a temporal v-structure does not imply a vanishing Transfer Entropy, as shown with the counterexample in Supplementary Fig. 5b. As a result, the presence of a temporal v-structure, Yt′ → Yt ← Xt in Gt, does not necessarily imply a vanishing transfer entropy, TY →X = 0, as long as there remains an edge between any Yt′′ and Xt, as in the example in Supplementary Fig. 5b. Hence, Granger-Schreiber causality is in fact too restrictive and may miss actual causal effects, which can be uncovered by structural causal discovery methods like tMIIC. In addition, Granger-Schreiber causality is also known to in- fer spurious causal associations by excluding the presence of latent common causes a priori . By constrast, CausalXtract’s causal dis- covery module includes time-delayed as well as synchronous effects originating from unobserved latent variables, as discussed above. Benchmarking of CausalXtract’s causal discovery module The performance of CausalXtract’s causal discovery module (tMIIC) has been assessed using Tigramite package 22, which provides differ- ent methods to learn temporal causal networks from time series data. We compared tMIIC to two methods capable of orienting contem- poraneous edges (PC and PCMCI+) and tested three different ker- nels for estimating mutual information (Parcorr, GPDC and KNN). Benchmark networks and datasets have been chosen to resemble the real-world data analyzed in this study ( i.e. similar number of time steps, network size and degree distribution) and include a large range of linear and non-linear relations between variables. A first series of datasets was generated for a 15 node benchmark network (Supplementary Fig. 1a) with linear combinations of contri- butions inspired by the Tigramite package, Supplementary Table 1. Running times and scores (Precision, Recall, F-score) have been aver- aged over 10 datasets (Supplementary Fig. 1b) and show that tMIIC scores are at par with PC and PCMCI+ using GPDC or KNN ker- nels but that tMIIC runs orders of magnitude faster, which enables to use tMIIC on much larger datasets of biological interest including a few tens or hundreds of thousands samples. Only PC or PCMCI+ using ParCorr kernel match tMIIC running speed but with signifi- cantly lower scores, as Fscores level off around 0.6-0.7 at large sample size, while tMIIC Fscore exceeds 0.9 (Supplementary Fig. 1b). Importantly, increasing the number of time-lagged layers from τ = 2 (as in the actual model, Supplementary Fig. 1a) to 5 or 10 layers in the inferred time-unfolded network (Supplementary Fig. 2) leads to very similar network reconstructions for simulated station- ary data. This demonstrates tMIIC insensitivity to an overestimated maximum lag for the reconstituted network. Interestingly, however, when the generated data is no longer stationary, increasing the num- ber of layers leads to multiple self-loops at non-stationary variables, whilst the rest of the network remains relatively unaffected (Supple- mentary Fig. 3). It demonstrates that CausalXtract’s causal discov- ery module is robust to the presence of non-stationary variables but requires long-time range interactions, and therefore multiple time- lagged layers, to account for these non-stationary dynamics at spe- cific variables. This striking observation on benchmark networks is also consistent with the multiple self-loops observed for a number of non-stationary variables in the real-world application on cellular ecosystems, Fig. 2a and Supplementary Fig. 6. A second series of more complex datasets was also generated for another 15 node benchmark network (Supplementary Fig. 4a) with non-linear combinations of contributors, Supplementary Ta- ble 2. Here, tMIIC tends to outperform both PC and PCMCI+, in terms of Recall and Fscores, while remaining orders of magnitude faster compared to GPDC and KNN kernels. Only PC or PCMCI+ using ParCorr kernel match tMIIC running speed but with signifi- cantly lower scores ( i.e. Fscores level off around 0.4-0.5 at large sam- ple size, while tMIIC Fscore exceeds 0.8). This demonstrates that CausalXtract’s causal discovery module (tMIIC) is both a reliable and scalable method to discover complex temporal causal relations in very large time series datasets including a few hundred thousand samples. Data availability The original live-cell time-lapse image data and extracted crops are available at: https://doi.org/10.5281/zenodo.7755699. Code availability The source code of the CausalXtract pipeline is available at: https://github.com/miicTeam/CausalXtract. It includes a demo R markdown notebook of CausalXtract pipeline, which reproduces step-by-step the results reported in the manuscript, Fig. 2, start- ing from the original live-cell time-lapse images of the tumor-on-chip ecosystem, Fig. 1a. Tigramite package used for benchmark compar- ison is available at: https://github.com/jakobrunge/tigramite

This work was supported by ITMO Cancer (grant No 20CM106) and the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 847718. LD acknowledges support from AMX PhD fellowship, VC from ARC foundation and NL from CNRS-Imperial College joint PhD programme. 6 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 8, 2024. ; https://doi.org/10.1101/2024.02.06.579177doi: bioRxiv preprint

1. Nguyen, M. et al. Dissecting Effects of Anti-cancer Drugs and Cancer-Associated Fibroblasts by On-Chip Reconstitution of Im- munocompetent Tumor Microenvironments. Cell Reports 25, 3884–3893.e3 (2018). 2. Davies, E. Machine vision 3rd ed. (Morgan Kaufmann, 2004). 3. Munkres, J. Algorithms for the assignment and transportation problems. J. Soc. Ind. Appl. Math. 5, 32–38 (1957). 4. Verny, L., Sella, N., Affeldt, S., Singh, P. P. & Isambert, H. Learn- ing causal networks with latent variables from multivariate infor- mation in genomic data. PLoS Comput. Biol. 13, e1005662 (2017). 5. Cabeli, V. et al. Learning clinical networks from medical records based on information estimates in mixed-type data. PLoS Com- put. Biol. 16, e1007866 (2020). 6. Da Câmara Ribeiro-Dantas, M. et al. Learning interpretable causal networks from very large datasets, application to 400,000 medical records of breast cancer patients. Preprint at https://doi.org/10.48550/arXiv.2303.06423 (2023). 7. Desterke, C. et al. Inferring Gene Networks in Bone Marrow Hematopoietic Stem Cell-Supporting Stromal Niche Populations. iScience 23, 101222 (2020). 8. Sella, N. et al. Interactive exploration of a global clinical network from a large breast cancer cohort. npj Digital Med 5, 113 (2022). 9. Assaad, C., Devijver, E. & Gaussier, E. Survey and Evaluation of Causal Discovery Methods for Time Series. Journal of Artificial Intelligence Research 73, 767–819 (2022). 10. Entner, D. & Hoyer, P. On Causal Discovery from Time Series Data using FCI. Proceedings of the 5th European Workshop on Probabilistic Graphical Models, PGM 2010 (2010). 11. Malinsky, D. & Spirtes, P. Causal Structure Learning from Mul- tivariate Time Series in Settings with Unmeasured Confounding in Proceedings of 2018 ACM SIGKDD Workshop on Causal Dis- covery, CD@KDD 2018 (eds Le, T. D., Zhang, K., Kiciman, E., Hyvärinen, A. & Liu, L.) 92 (2018), 23–47. 12. Runge, J., Nowack, P., Kretschmer, M., Flaxman, S. & Sejdinovic, D. Detecting and quantifying causal associations in large nonlinear time series datasets. Science Advances 5 (2019). 13. Runge, J. et al. Inferring causation from time series in Earth sys- tem sciences. en. Nat. Commun. 10, 2553 (2019). 14. Granger, C. W. J. Investigating causal relations by econometric models and cross-spectral methods. Econometrica 37, 424 (1969). 15. Schreiber, T. Measuring Information Transfer. Physical Review Letters 85, 461–464 (2000). 16. Barnett, L., Barrett, A. B. & Seth, A. K. Granger Causality and Transfer Entropy Are Equivalent for Gaussian Variables. Physical Review Letters 103 (2009). 17. Driscoll, M. K. & Zaritsky, A. Data science in cell imaging. Journal of Cell Science 134, jcs254292 (2021). 18. Chan, T. F. & Vese, L. A. Active contours without edges. IEEE Trans. Image Process. 10, 266–277 (2001). 19. Masuzzo, P., Van Troys, M., Ampe, C. & Martens, L. Taking aim at moving targets in computational cell migration. Trends Cell Biol. 26, 88–110 (2016). 20. Pearl, J. Causality (Cambridge university press, 2009). 21. Spirtes, P., Glymour, C. N., Scheines, R. & Heckerman, D. Cau- sation, prediction, and search (MIT press, 2000). 22. Runge, J. Discovering contemporaneous and lagged causal rela- tions in autocorrelated nonlinear time series datasets in Proceed- ings of the 36 th Conference on Uncertainty in Artificial Intelli- gence (eds Peters, J. & Sontag, D.) 124 (2020), 1388–1397. 7 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 8, 2024. ; https://doi.org/10.1101/2024.02.06.579177doi: bioRxiv preprint Supplementary Fig. 1: Benchmark assessment of CausalXtract’s causal discovery module (tMIIC) using generated time series datasets. a, Example of a 15 node causal network to generate benchmark time series datasets based on linear combinations of contributions, Supplementary Table 1. Examples of temporal causal networks reconstructed by tMIIC based on 100, 1,000 or 10,000 simulated time steps. b, Running times and scores (Precision, Recall, Fscore) averaged over 10 datasets and compared to PC and PCMCI+ methods using different kernels (GPDC, KNN, ParCorr); tMIIC is at par with PC and PCMCI+ scores using GPDC and KNN kernels but runs orders of magnitude faster. Only ParCorr kernel matches tMIIC running speed but with significantly lower scores at large sample size, see Methods. 8 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 8, 2024. ; https://doi.org/10.1101/2024.02.06.579177doi: bioRxiv preprint Supplementary Fig. 2: CausalXtract insensitivity to an overestimated maximum lag τ . a, Example of a temporal causal network model with a maximum lag τ = 2 . Corresponding temporal causal networks inferred by CausalXtract’s causal discovery module (tMIIC), from 1,000 time step stationary time series (Supplementary Table 1), while assuming different maximum lags τ = 2 , 5 or 10. b, Running times and scores (Precision, Recall, Fscore) of tMIIC temporal causal network reconstructions for τ = 2 , 5 or 10, averaged over ten stationary time series of 10 to 105 time steps. Overestimating the maximum lag τ has little impact on the reconstructed networks, as long as the time series are stationary, as demonstrated in Supplementary Fig. 3. 9 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 8, 2024. ; https://doi.org/10.1101/2024.02.06.579177doi: bioRxiv preprint Supplementary Fig. 3: CausalXtract sensitivity to non-stationary variables. a, Example of a temporal causal network model (τ = 2) with a low frequency periodic input (T = 100 ) applied to X8 and a time-linear trend applied to X13. Corresponding temporal causal networks inferred by tMIIC from 1,000 time step time series (Supplementary Table 1) including non-stationary inputs to X8 and X13. Increasing the maximum lag from τ = 2 to τ = 5 or 10 leads to the appearence of multiple self-loops, which result from the non-stationary dynamics of X8 and X13, whilst the rest of the network remains largely unaffected. b, Running times and scores (Precision, Recall, Fscore ignoring X8 and X13 self-loops) of tMIIC causal network reconstructions for τ = 2, 5 or 10, averaged over ten time series of 10 to 105 time steps. 10 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 8, 2024. ; https://doi.org/10.1101/2024.02.06.579177doi: bioRxiv preprint Supplementary Fig. 4: Benchmark assessment of CausalXtract’s causal discovery module (tMIIC) using more complex time series datasets. a, Example of a 15 node causal network to generate more complex benchmark time series datasets based on non-linear combinations of contributions, Supplementary Table 2. Examples of temporal causal networks reconstructed by tMIIC based on 100, 1,000 or 10,000 simulated time steps. b, Running times and scores (Precision, Recall, Fscore) averaged over 10 datasets and compared to PC and PCMCI+

using different kernels (GPDC, KNN, ParCorr); tMIIC outperforms both PC and PCMCI+, in terms of Recall and Fscores, while running orders of magnitude faster, except for the ParCorr kernel, which leads, however, to significantly lower scores at large sample size. 11 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 8, 2024. ; https://doi.org/10.1101/2024.02.06.579177doi: bioRxiv preprint Yt X t Yt’<t X t’<tI( ; | ) = 0 X X tt’<t Y Ytt’<t Yt X t X t Yt’<t X t’<tI( ; | ) = 0 Yt’ Y Y X X t tt’<t t’<t X t Yt’ a b Supplementary Fig. 5: Time-unfolded causal network framework and relation to Granger-Schreiber temporal causality. a , A vanishing Transfer Entropy, i.e. TY →X = I(Xt; Yt′<t|Xt′<t) = 0 , implies i) the absence of (dashed) edge between Xt and any Yt′ , with t′ < t, and ii) if Xt is adjacent to Yt, the presence of temporal (2-variable + time) v-structures, Yt′ → Yt ← Xt, for all Yt′ adjacent to Yt, with t′ < t (Methods, Theorem 1). These results can be readily extended to include the presence of other observed variables, Vt′⩽t, by redefining Transfer Entropy as, TY →X = I(Xt; Yt′<t|Xt′<t, Vt′⩽t), which discards contributions from indirect paths through other observed variables, Vt′⩽t. b, By contrast, the presence of a temporal (2-variable + time) v-structure, Yt′ → Yt ← Xt does not imply a vanishing Transfer Entropy, as long as there remains an edge between any Yt′′<t and Xt. It implies that Granger-Schreiber temporal causality is in fact too restrictive and may overlook actual causal effects, which can be uncovered by graph-based causal discovery methods like CausalXtract’s causal discovery module (tMIIC). Hence, CausalXtract’s time-unfolded network framework, combining graph-based and information-based approaches, sheds light on the common foundations of the seemingly unrelated graph-based causality and Granger-Schreiber temporal causality, while clarifying their actual differences and limitations. 12 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 8, 2024. ; https://doi.org/10.1101/2024.02.06.579177doi: bioRxiv preprint Supplementary Fig. 6: Time series of cellular features extracted from the tumor ecosystems. Example of time series of cellular features extracted by CausalXtract’s feature extraction module (CellHunter+) from the tumor ecosystems analyzed in this study, Fig. 1a. It includes two experimental control parameters ( i.e. treatment and CAF presence) and 15 cellular features extracted every 2 minutes over a period of two days. Continuous features are highlighted for one trajectory (traj.18), while categorical features are shown for all trajectories. 13 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 8, 2024. ; https://doi.org/10.1101/2024.02.06.579177doi: bioRxiv preprint a b Supplementary Fig. 7: Time-unfolded causal network inferred by CausalXtract. a , Time-unfolded causal network assuming stationary dynamics of cellular ecosystems implying translational time invariance of the inferred causal network. b, Only edges involving at least one contemporaneous variables (i.e. at time t) need to be tested for conditional independence by tMIIC and the remaining edges are then duplicated at all previous time steps before assigning orientations when time-lagged latent variables are taken into account, Fig. 1c. Variables retaining multiple self-loops with different time-delays correspond to non-stationary variables in Supplementary Fig. 6, in agreement with benchmarks from simulated data including non-stationary variables, Supplementary Fig. 3. 14 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 8, 2024. ; https://doi.org/10.1101/2024.02.06.579177doi: bioRxiv preprint b a c Supplementary Fig. 8: Robustness of CausalXtract’s temporal causal networks to variations in sampling rate. Summary causal networks inferred by CausalXtract using different sampling rates (δτ ). a, δτ = 8 ts and τ = 80 ts, in time step units (1 ts = 2 min). b, δτ = 7 ts, and τ = 84 ts, as chosen automatically by CausalXtract based on the average relaxation time across the 15 monitored variables, τR = 40 ts, which defines a maximum lag τ = 2 τR = 80 ts. Given a total number of (time-lagged and -unlagged) nodes, chosen to be around 200 nodes for computational efficiency, it leads to 13 temporal layers (ν + 1 = 200 /15 ≃ 13) and a lag increment δτ = τ /ν ≃ 7 ts. This summary causal network corresponds to Fig. 2a. c, δτ = 5 ts and τ = 60 ts, corresponding to τ = ν · δτ with ν + 1 = 13 temporal layers, as in (b). 15 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 8, 2024. ; https://doi.org/10.1101/2024.02.06.579177doi: bioRxiv preprint Supplementary Table 1: 15 nodes model. Nodes X 1 t ← −0.47 f2(X 1 t−1) + 0.29 f3(X 2 t−1) × η1 X 2 t ← 0.49 f2(X 2 t−1) + 0.4 f1(X 1 t−2) + η2 X 3 t ← 0.56 f1(X 3 t−1) + 0.44 f4(X 4 t−2) − 0.26 f2(X 10 t−2) + 0.56 f2(X 4 t ) + η3 X 4 t ← 0.24 f3(X 4 t−1) − 0.24 f2(X 6 t−2) − 0.12 f4(X 14 t−1) × η4 X 5 t ← −0.39 f3(X 5 t−1) − 0.42 f3(X 5 t−2) − 0.39 f3(X 11 t ) + η5 X 6 t ← −0.32 f2(X 6 t−1) + η6 X 7 t ← −0.17 f4(X 7 t−1) − 0.17 f1(X 7 t−2) + η7 X 8 t ← 0.39 f4(X 8 t−1) − 0.46 f4(X 7 t−1) − 0.39 f3(X 1 t−1) − 0.4 f3(X 12 t−2) + η8 X 9 t ← −0.34 f1(X 9 t−1) + 0.43 f3(X 12 t−2) + η9 X 10 t ← 0.2 f1(X 10 t−1) + 0.18 f4(X 9 t−2) + 0.17 f1(X 9 t−1) + 0.48 f3(X 7 t−1) − 0.26 f4(X 4 t−1) + η10 X 11 t ← 0.41 f2(X 11 t−1) + 0.54 f3(X 2 t ) − 0.55 f2(X 12 t ) + η11 X 12 t ← −0.45 f2(X 12 t−1) − 0.43 f4(X 3 t−2) − 0.17 f4(X 9 t−2) × η12 X 13 t ← 0.45 f3(X 13 t−1) + η13 X 14 t ← 0.28 f2(X 14 t−1) + 0.37 f1(X 12 t−2) × η14 X 15 t ← 0.52 f3(X 15 t−1) + η15 F unctions f1(x) = x f2(x) = x (1 − 4 e− x2 2 ) f3(x) = x (1 − 4 x3 e− x2 2 ) f4(x) = cos(x) Noises The η are white noises generated for each node or contribution using a normal distribution: η ∼ N (0, 1) 16 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 8, 2024. ; https://doi.org/10.1101/2024.02.06.579177doi: bioRxiv preprint Supplementary Table 2: 15 nodes model with combinations. Nodes X 1 t ← η − 0.7 f6(u(η + X 1 t−1)) − 0.87 f5(u(η + (X 14 t−1 × X 1 t−2))) X 2 t ← η + 0.65 f1(u(η + X 2 t−1)) − 0.63 f3(u(η + X 2 t−2)) + 0.79 f3(u(η + X 5 t−1)) X 3 t ← η − 0.76 f5(u(η + X 3 t−1)) − 0.59 f6(u(η + X 7 t−1)) − 0.85 f2(u(η + X 15 t−1)) −0.89 f5(u(η + (X 13 t−2 × X 7 t−1))) X 4 t ← η − 0.7 f6(u(η + X 5 t−1)) − 0.86 f2(u(η + X 8 t−2)) + 0.53 f1(u(η + (X 4 t−1 × X 9 t−2))) X 5 t ← η + 0.54 f2(u(η + (X 14 t−1 × X 6 t−2))) X 6 t ← η − 0.85 f2(u(η + X 6 t−1)) − 0.79 f3(u(η + X 3 t−2)) + 0.59 f1(u(η + X 4 t−1)) +0.75 f3(u(η + X 1 t )) + 0.57 f2(u(η + X 14 t−1)) X 7 t ← η + 0.74 f1(u(η + X 7 t−1)) + 0.54 f6(u(η + X 9 t−1)) − 0.53 f2(u(η + (X 9 t−1 × X 7 t−1))) X 8 t ← η × (−0.63 f1(u(η + X 6 t−1)) + 0.81 f5(u(η + X 13 t )) + 0.53 f6(u(η + (X 6 t−2 × X 6 t−1))) −0.69 f6(u(η + (X 13 t × X 6 t−1)))) X 9 t ← η + 0.79 f3(u(η + X 4 t−2)) + 0.69 f6(u(η + (X 9 t−1 × X 15 t−1))) X 10 t ← η + 0.54 f6(u(η + X 10 t−1)) X 11 t ← η + 0.83 f6(u(η + X 11 t−1)) − 0.76 f4(u(η + X 13 t−1)) − 0.73 f3(u(η + X 2 t−1)) +0.74 f2(u(η + X 4 t )) − 0.87 f2(u(η + X 10 t−2)) + 0.72 f4(u(η + X 12 t−1)) −0.73 f1(u(η + (X 10 t−2 × X 13 t−1))) X 12 t ← η + 0.7 f3(u(η + X 10 t−1)) − 0.55 f5(u(η + X 9 t )) − 0.54 f5(u(η + (X 12 t−1 × X 10 t−1))) X 13 t ← η − 0.62 f3(u(η + X 14 t−2)) − 0.61 f1(u(η + (X 13 t−1 × X 14 t−2))) X 14 t ← η − 0.78 f6(u(η + X 14 t−1)) X 15 t ← η − 0.68 f4(u(η + X 15 t−1)) + 0.85 f4(u(η + X 15 t−2)) − 0.6 f5(u(η + X 10 t−2)) +0.68 f6(u(η + X 14 t−1)) + 0.81 f4(u(η + (X 14 t−1 × X 10 t−2))) F unctions u(x) = max(−1, min(1, x)) f1(x) = x f2(x) = x (1 − 4 e− x2 2 )/1.52387 f3(x) = 4 x2 f4(x) = 8 x3 f5(x) = 16 x4 f6(x) = cos(πx) Noises The η are white noises generated for each node or contribution using a normal distribution: η ∼ N (0, 0.1) 17 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 8, 2024. ; https://doi.org/10.1101/2024.02.06.579177doi: bioRxiv preprint