Live-cell imaging microscopy commonly produces extensive amounts of time-lapse images of cellular systems, which can be segmented to extract morphodynamic features and interactions of individual cells under increasingly complex and physiologically relevant conditions. However, this wealth of information remains largely under-exploited due to a lack of methods and tools able to discover causal effects from spatio-temporal correlations under well-controlled experimental conditions.

CausalXtract addresses this need by integrating an advanced live-cell image feature extraction tool with a reliable and scalable causal discovery method, Fig. 1, in order to learn temporal causal networks from live-cell time-lapse imaging data, Fig. 2.

Fig. 1:

Fig. 2:

Results 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¹ and consists in three steps: detection, tracking and feature extraction of live cells within time-lapse video images. First, automatic localization/segmentation of cells (e.g. tumor and immune cells) is performed with the Circular Hough Transform (CHT) algorithm² to estimate the cell centers and radii. Second, cell trajectories along the frames are constructed by linking the positions detected at the previous time step through Munkres’ algorithm for Optimal sub-pattern Assignment Problems (OAPs)³. 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 temporal information is not available) for a broad range of biological or biomedical data, from single-cell transcriptomic and genomic alteration data^4,7 to medical records of patients^5,6. Live-cell time-lapse imaging data contain, however, information 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, reconstructs time-unfolded causal networks, where each variable is represented by several nodes at different relative time points⁹, Fig. 1c. Such a time-unfolded network framework^10–13 is required to account for the temporal correlation between successive time steps in time series data. We benchmarked tMIIC on synthetic datasets resembling the real-world data of interest analyzed 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¹⁴ for linear time series without reference to graphical models and later extended to non-linear dynamics by Schreiber^15,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 methods, Supplementary Fig. 5 (Methods, Theorem 1). In addition, Granger-Schreiber causality has long been known to infer spurious causal associations based on time delays, by excluding the presence of latent common causes a priori⁹. CausalXtract circumvents these limitations by combining graph-based and information-based approaches (Methods), while including contemporary and time-delayed effects of unobserved latent variables, 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 images 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¹ which demonstrated the effects of an anti-cancer drug (the monoclonal antibodies trastuzumab, brand name Herceptin, used to treat HER2+ breast cancers) on a reconstituted tumor microenvironment including cancer cells, immune cells, cancer-associated fibroblasts (CAF), and endothelial cells (Methods). However, a comprehensive extraction and analysis of cellular morphodynamic features and interactions remained unexplored.

To this end, cellular features such as cell geometry, velocity, division, apoptosis, cell-cell transient interactions and persistent contacts were first extracted from the raw images using CausalXtract’s feature extraction module, Fig. 1b and Supplementary Fig. 6. Then, a time-unfolded causal network, Supplementary Fig. 7, and the corresponding summary causal network, 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 biologically relevant findings, in addition to confirming known results fxrom earlier studies. In particular, CausalXtract discovers that CAFs directly inhibit cancer cell apoptosis, independently from anti-cancer treatment, Fig. 2b, while earlier studies reported that CAFs merely reduced the effect of treatment¹. 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 earlier studies. In particular, it recovers that treatment increases cancer cell apoptosis and the number of cancer-immune interactions, 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 possibly 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 eventually 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 independent computational modules, conceived to warrant interoperability with alternative live-cell segmentation and tracking methods 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 needed¹⁷ and we believe that CausalXtract can bring unique insights based on causal discovery to interpret such informationrich live-cell imaging data.

Materials and Methods

Tumor-on-chip preparation and live-cell microscopy

Videos analyzed in the present study refer to biological experiments emulating a 3D breast tumor ecosystem¹. All tumor-on-chip experiments have a central endothelium compartment containing endothelial 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 cancer BT474 cell line) and immune cells (peripheral blood mononuclear cells, PBMCs) from healthy donors, Fig. 1a. Four experimental 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 cancer 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 center 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¹ and consists in three steps: detection, tracking and feature extraction of live cells within timelapse video images. First, cell detection is based on the segmentation of circular-shaped objects using CHT² with radii set around the theoretical radii of the two cell populations (r_im = 4 px for immune cells and r_ca = 14 px for MCCs with a pixel resolution 1 px = 0.645 μm¹). Then, cell tracking is performed by linking cells detected at the i^th frame to cells located at the (i + 1)^th frame within a maximum distance 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 characteristics¹. 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³ 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 i^th 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 interaction features are extracted at successive positions along each trajectory. For each MCC, 15 descriptors were extracted, Supplementary Fig. 6, and classified into four main categories: cell shape, motility, state, and interaction descriptors.

Shape descriptors

The active contour algorithm implemented in Matlab¹⁸ was used to segment the MCC boundaries on each video crop frame. Taking as input a frame representing the i^th snapshot of the t^th MCC, it returns a binary image, where the MCC is represented by a white region. From the binary image, the shape properties of the region occupied by each MCC were extracted using the Matlab regionprops algorithm. The resulting descriptors of the extracted shape are listed below:

• area indicates the number of pixels composing the region. The equivalent diameter of the t^th MCC in the i^th frame is defined as.

• perimeter represents the distance along the MCC boundary.

• circularity is defined as 4 · area · π/perimeter², 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, , corresponding to the difference in absolute value of the equivalent diameters between the i^th and (i − 1)^th frames of the t^th MCC.

Motility descriptors

The positions and of the t^th MCC in the i^th and (i − 1)^th frames were compared using the Euclidean distance d(·) to define the following motility parameters:

• instantaneous cancer velocity¹⁹ is defined as, where Δt is the time interval between two consecutive frames.

• net displacement¹⁹ indicates the resultant distance between the initial and current positions of the t^th MCC,.

• directionality¹⁹ is defined as the ratio of net displacement, , and curvilinear distance,. It measures 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 experiment. 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, r₁ = r_im + r_ca + 2 = 20 px and r₂ = 2 × (r_im + r_ca) = 36 px¹. Hence, r₁ refers to MCC and immune cells in actual physical contact, while r₂ 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₂) corresponds to the number of immune cells within the interaction radius r₂ around the MCC on that frame.

• number of cancer-immune interactions (r₁) corresponds to the number of immune cells in close contact with the MCC on that frame.

• minimal cancer-immune distance (r₂) is the minimum distance between the MCC and the immune cells within a radius r₂.

• mean immune velocity (r₂) is the mean instantaneous velocity norm of the immune cells within the interaction radius r₂ around the MCC.

• mean immune velocity (r₁) 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 networks 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 network. The second step then consists in orienting some of the retained edges based on the signature of causality in observational data. This corresponds to orienting three-variable “v-structure” motifs as, X → Z ← Y, whenever the edge X − Y has been removed without conditioning on the variable Z, which implies that Z cannot 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 ⤎ L ⤏ Z. In addition, classical causal discovery 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 significant information contributors before assessing conditional independences. 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⁵ and also heterogeneous multimodal data, by analyzing continuous and categorical variables on the same footing, based on a mutual information supremum principle for finite dataset^5,6. Last, MIIC distinguishes genuine causal relations from putative and latent causal effects⁶, 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 discovery module (tMIIC) aims to learn a time-unfolded graph, G_t, where each variable is represented by a series of nodes associated to its value at different relative time points, Fig. 1c. Such a timeunfolded network framework^10–13 is required to account for the temporal correlation between successive samples in time series data. Assuming that the dynamics can be considered stationary (see Benchmarking of CausalXtract’s causal discovery module section, below), the time-unfolded graph, 𝒢_t, should be translationally invariant over time and can be assigned a periodic structure a priori. In addition, G_t 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 G_t by introducing a time increment δτ between consecutive time points, which leads to ν = τ/δτ time-lagged layers in 𝒢_t.

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 necessary edges have been identified, they are simply duplicated at earlier time points to enforce the translational invariance of 𝒢_t skeleton. Time-lagged edges are then pre-oriented with a first arrowhead pointing 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 applied to 𝒢_t, 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 time series data was first proposed for time-lagged effects¹⁰ and subsequently extended to contemporaneous effects¹¹ by adapting the constraint-based FCI method allowing for latent variables²¹. While traditional constraint-based methods suffer from poor recall, the recent PCMCI¹² / PCMCI+²² method improves recall by introducing ad hoc conditioning rules for auto-correlated time series. By contrast, 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 variables.

Relation to Granger-Schreiber temporal causality

The concept of temporal causality was originally formulated by Granger¹⁴ without reference to any graphical model by comparing linear autoregression with or without past values of possible causal variables. This was later extended to non-linear relations by Schreiber^15,16 using the notion of Transfer Entropy, T_X→Y, which can be expressed in terms of multivariate conditional information,

where and denote the sets of variables, and, taken at earlier time points t^’ than t.

While Eq. 1 is asymmetric upon X/Y permutation, a simple comparison of Transfer Entropy asymmetry (e.g. T_X→Y > T_{Y → ⩾} 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. T_Z→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 models (𝒢 _t). This highlights the common foundations of temporal and structural causalities beyond their seemingly unrelated definitions.

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, in 𝒢 _t, does not necessarily imply a vanishing transfer entropy, T_{Y →X} = 0, as long as there remains an edge between any and X_t, 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 infer spurious causal associations by excluding the presence of latent common causes a priori. By constrast, CausalXtract’s causal discovery 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²², which provides different methods to learn temporal causal networks from time series data. We compared tMIIC to two methods capable of orienting contemporaneous edges (PC and PCMCI+) and tested three different kernels 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 contributions inspired by the Tigramite package, Supplementary Table 1. Running times and scores (Precision, Recall, F-score) have been averaged over 10 datasets (Supplementary Fig. 1b) and show that tMIIC scores are at par with PC and PCMCI+ using GPDC or KNN kernels 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 significantly 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 stationary 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 number of layers leads to multiple self-loops at non-stationary variables, whilst the rest of the network remains relatively unaffected (Supplementary Fig. 3). It demonstrates that CausalXtract’s causal discovery module is robust to the presence of non-stationary variables but requires long-time range interactions, and therefore multiple timelagged layers, to account for these non-stationary dynamics at specific 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 Table 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 significantly lower scores (i.e. Fscores level off around 0.4-0.5 at large sample 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, starting from the original live-cell time-lapse images of the tumor-on-chip ecosystem, Fig. 1a. Tigramite package used for benchmark comparison is available at: https://github.com/jakobrunge/tigramite

Acknowledgements

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.