CausalXtract’s live-cell image feature extraction module (CellHunter+) (Figure 1b) is based on CellHunter software (Nguyen et al., 2018) and consists of 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 (Davies, 2004) 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) (Munkres, 1957). 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 (‘Materials and methods’).

CausalXtract’s temporal causal discovery module (tMIIC) (Figure 1c) is adapted from the causal discovery method (MIIC) (Verny et al., 2017; Cabeli et al., 2020; Cabeli and Li, 2021; Ribeiro-Dantas et al., 2024), 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 (Verny et al., 2017; Desterke et al., 2020) to medical records of patients (Cabeli et al., 2020; Sella et al., 2022; Ribeiro-Dantas et al., 2024). 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 (Assaad et al., 2022; Figure 1c). Such a time-unfolded network framework (Entner and Granger, 2010; Malinsky and Spirtes, 2018; Runge et al., 2019) is required to account for the temporal correlation between successive time steps in time-series data. This graph-based causal framework goes beyond the seminal concept of temporal causality originally proposed by Granger, 1969 for linear time series without reference to graphical models and later extended to nonlinear dynamics by Schreiber, 2000; Barnett et al., 2009. 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 (Figure 2a and b; ‘Materials and 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 (Assaad et al., 2022). tMIIC circumvents these limitations by combining graph-based and information-based approaches (‘Materials and 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).

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, PC and PCMCI+ (Runge, 2020), while running order of magnitudes faster on datasets of biologically relevant size including tens to hundreds of thousands time steps (Figure 2c and d and Figure 2—figure supplement 3; Figure 2—figure supplement 2; Figure 2—figure supplement 1; Figure 2—figure supplements 4).

We showcase CausalXtract with the analysis of time-lapse images of a tumor ecosystem reconstituted ex vivo using the tumor-on-chip technology (Figure 1a). These live-cell time-lapse images come from a proof-of-concept study (Nguyen et al., 2018), which demonstrated the effects of an anticancer 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 (CAFs), and endothelial cells (‘Materials and 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 (Figure 1b, Figure 1—figure supplement 1).

Then, summary causal network (Figure 3a) and the corresponding time-unfolded causal network (Figure 3—figure supplement 1) were reconstructed between extracted cellular features, cell–cell interactions, and therapeutic conditions using CausalXtract’s temporal causal discovery module (Figure 1c).

CausalXtract inferred network (Figure 3a) uncovers novel biologically relevant findings, in addition to confirming known results from earlier studies. In particular, CausalXtract discovers that CAFs directly inhibit cancer cell apoptosis, independently from anticancer treatment (Figure 3b), while earlier studies reported that CAFs merely reduced the effect of treatment (Nguyen et al., 2018). CausalXtract also discovers that treatment increases cancer cell perimeter (Figure 3c), 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 (Figure 3c). Likewise, CausalXtract recovers that CAFs stimulate cancer cell migration and increase their area (Figure 3b).

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 (Figure 3d). 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, 2–3 hr 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) (Figure 3d). Likewise, the area change upon apoptosis is predicted to first decrease soon after apoptosis (blue edge) before eventually increasing upon cell lysis (red edge) (Figure 3d). These results are robust to variations in sampling rate (Figure 3—figure supplement 2).

The videos analyzed in the present study refer to biological experiments emulating a 3D breast tumor ecosystem (Nguyen et al., 2018). 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 (Figure 1a). Four experimental conditions were considered depending on the presence or absence of breast CAFs (CAF cell line Hs578T) and drug treatment (trastuzumab, Herceptin). The immortalized human BT474 and Hs578T cell lines were purchased from ATCC (#HTB-20, #HTB-126) and authenticated by SRT profiling (GenePrint 10 system, Promega, #B9510). The human primary HUVEC lines were purchased from Lonza (#C2517A). PBMCs were routinely isolated from the fresh blood of healthy donors by density gradient centrifugation. All cells were periodically tested to exclude mycoplasma contamination using a qPCR-based method (VenorGem Classic, BioValley, #11-1250). Videos were acquired using inverted motorized Leica microscopes with a frame rate of 2 min for up to 48 hr (1440 frames). Figure 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. Thirty-six video crops of up to 1440 frames were analyzed (46,935 frames in total) corresponding to nine video crops per experimental conditions.

The live-cell image feature extraction module (CellHunter+) (Figure 1b) extends the CellHunter software (Nguyen et al., 2018) and consists of 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 (Davies, 2004) with radii set around the theoretical radii of the two cell populations (${\displaystyle r_{\mathrm{i}\mathrm{m}}=4}$ px for immune cells and ${\displaystyle r_{\mathrm{c}\mathrm{a}}=14}$ px for MCCs with a pixel resolution ${\displaystyle 1\mathrm{p}\mathrm{x}{\displaystyle =0.645\mu \mathrm{m}}}$; Nguyen et al., 2018). Then, cell tracking is performed by linking cells detected at the ${\displaystyle i^{th}}$ frame to cells located at the ${\displaystyle (i+1{)}^{th}}$ frame within a maximum distance from the detected cell candidate. While the motions of both MCCs and immune cells resemble random walks with time-varying drift and volatility, these two cell types exhibit different motility characteristics (Nguyen et al., 2018). 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 Munkres, 1957 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 ${\displaystyle i^{th}}$ and ${\displaystyle (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 (Figure 1—figure supplement 1) and classified into four main categories: cell shape, motility, state, and interaction descriptors.

Traditional causal discovery methods (Pearl, 2009; Spirtes, 2000) 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 of 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 ${\displaystyle X\to Z\leftarrow Y}$ whenever the edge ${\displaystyle X-Y}$ has been removed without conditioning on the variable ${\displaystyle Z}$, which implies that ${\displaystyle Z}$ cannot be a cause of ${\displaystyle X}$ nor ${\displaystyle Y}$. This does not guarantee, however, that ${\displaystyle X}$ (or ${\displaystyle Y}$) is an actual cause of ${\displaystyle Z}$, which also requires to rule out the possibility that the edge between ${\displaystyle X}$ and ${\displaystyle Z}$ (or ${\displaystyle Y}$ and ${\displaystyle Z}$) might originate from a latent common cause, ${\displaystyle L}$, unobserved in the dataset, that is, ${\displaystyle 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) (Verny et al., 2017; Cabeli et al., 2020; Cabeli and Li, 2021; Ribeiro-Dantas et al., 2024), 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 (Affeldt and Isambert, 2015; Affeldt et al., 2016). In practice, MIIC’s strategy limits spurious conditional independences, which improves its edge sensitivity and orientation reliability compared to traditional constraint-based methods. In addition, MIIC can handle missing data (Cabeli et al., 2020) 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 (Cabeli et al., 2020; Cabeli and Li, 2021; Ribeiro-Dantas et al., 2024). Last, MIIC distinguishes genuine causal relations from putative and latent causal effects (Ribeiro-Dantas et al., 2024) that are ubiquitous in real-world applications.

In order to analyze time-series datasets, CausalXtract’s causal discovery module (tMIIC) aims to learn a time-unfolded graph, ${\displaystyle {\mathcal{G}}_t}$, where each variable is represented by a series of nodes associated with its value at different relative time points (Figure 1c). Such a time-unfolded network framework (Entner and Granger, 2010; Malinsky and Spirtes, 2018; Runge et al., 2019) is required to account for the temporal correlation between successive samples in time-series data. Assuming that the dynamics can be considered stationary (see section ‘Benchmarking of CausalXtract’s causal discovery module’), the time-unfolded graph, ${\displaystyle {\mathcal{G}}_t}$, should be translationally invariant over time and can be assigned a periodic structure a priori. In addition, ${\displaystyle {\mathcal{G}}_t}$ can be restricted to a few time steps from the running time, ${\displaystyle t}$, back to a maximum time lag, ${\displaystyle t-\tau }$, since nodes at future time points (${\displaystyle t^{\mathrm{\prime }}>t}$) cannot a priori influence the observed data at current or previous time points (${\displaystyle t^{\mathrm{\prime }}⩽t}$) (Figure 1c). The maximum time lag ${\displaystyle \tau }$ should be chosen so as to have little effect on the final graphical model, which can be achieved for instance by setting ${\displaystyle \tau }$ to twice the average relaxation time of the variables of the dataset. In practice, we may also limit the number of time points ${\displaystyle \nu }$ in ${\displaystyle {\mathcal{G}}_t}$ by introducing a time increment ${\displaystyle \delta \tau }$ between consecutive time points, which leads to ${\displaystyle \nu =\tau /\delta \tau }$ time-lagged layers in ${\displaystyle {\mathcal{G}}_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 ${\displaystyle t}$ (Figure 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 ${\displaystyle {\mathcal{G}}_t}$ skeleton. Time-lagged edges are then pre-oriented with a first arrowhead pointing toward 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 ${\displaystyle {\mathcal{G}}_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 (Figure 1c).

Learning such structural models including latent variables from time-series data was first proposed for time-lagged effects (Entner and Granger, 2010) and subsequently extended to contemporaneous effects (Malinsky and Spirtes, 2018) by adapting the constraint-based FCI method allowing for latent variables (Spirtes, 2000). While traditional constraint-based methods suffer from poor recall, the recent PCMCI (Runge et al., 2019)/PCMCI+ (Runge, 2020) method improves recall by introducing ad hoc conditioning rules for autocorrelated 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.

The concept of temporal causality was originally formulated by Granger, 1969 without reference to any graphical model by comparing linear autoregression with or without past values of possible causal variables. This was later extended to nonlinear relations by Schreiber, 2000; Barnett et al., 2009 using the notion of Transfer Entropy, ${\displaystyle T_{X\to Y}}$, which can be expressed in terms of multivariate conditional information:

where ${\displaystyle {\mathit{X}}_{t^{\mathrm{\prime }}<t}}$ and ${\displaystyle {\mathit{Y}}_{t^{\mathrm{\prime }}<t}}$ denote the sets of variables, ${\displaystyle {\displaystyle X_{t^{\mathrm{\prime }}}}}$ and ${\displaystyle {\displaystyle Y_{t^{\mathrm{\prime }}}}}$, taken at earlier time points ${\displaystyle {\displaystyle t^{\mathrm{\prime }}}}$ than ${\displaystyle {\displaystyle t}}$.

While Equation 1 is asymmetric upon ${\displaystyle X}$/${\displaystyle Y}$ permutation, a simple comparison of Transfer Entropy asymmetry (e.g., ${\displaystyle T_{X\to Y}>T_{Y\to 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.,${\displaystyle T_{Z\to X}\approx 0}$), which suggests the possibility of a causal relation in the opposite direction, ${\displaystyle X\to 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 this relation between temporal causality without reference to any structural model (Equation 1) and structural causality entailed by time-unfolded causal graphical models (${\displaystyle {\mathcal{G}}_t}$). This highlights the common foundations of temporal and structural causalities beyond their seemingly unrelated definitions.

Theorem 1. [${\displaystyle {\displaystyle T_{Y\to X}=0}}$ implies temporal (2-variable + time) v-structures]

If ${\displaystyle {\displaystyle X_t}}$ is adjacent to ${\displaystyle {\displaystyle Y_t}}$ in ${\displaystyle {\displaystyle {\mathcal{G}}_t}}$ and ${\displaystyle T_{Y\to X}=I(X_t;{\mathit{Y}}_{t^{\mathrm{\prime }}<t}|{\mathit{X}}_{t^{\mathrm{\prime }}<t})=0}$, then for all ${\displaystyle {\displaystyle Y_{t^{\mathrm{\prime }}}}}$ adjacent to ${\displaystyle {\displaystyle Y_t}}$ in ${\displaystyle {\displaystyle {\mathcal{G}}_t}}$ , with ${\displaystyle {\displaystyle t^{\mathrm{\prime }}<t}}$, there is a temporal (2-variable + time) v-structure, ${\displaystyle {\displaystyle Y_{t^{\mathrm{\prime }}}\to Y_t\leftarrow X_t}}$ , in ${\displaystyle {\displaystyle {\mathcal{G}}_t}}$ (Figure 2a).

Proof: If ${\displaystyle T_{Y\to X}=I(X_t;{\mathit{Y}}_{t^{\mathrm{\prime }}<t}|{\mathit{X}}_{t^{\mathrm{\prime }}<t})=0}$, then all pairs ${\displaystyle (X_t,Y_{t^{\mathrm{\prime }}})}$ should be unconnected (assuming ‘faithfulness’, i.e., no coincidental cancellation of effects) and all unshielded triples ${\displaystyle Y_{t^{\mathrm{\prime }}}-Y_t-X_t}$ should be temporal v-structures, ${\displaystyle Y_{t^{\mathrm{\prime }}}\to Y_t\leftarrow X_t}$ , as ${\displaystyle Y_t\notin {\mathit{X}}_{t^{\mathrm{\prime }}<t}}$ in ${\displaystyle T_{Y\to X}=I(X_t;{\mathit{Y}}_{t^{\mathrm{\prime }}<t}|{\mathit{X}}_{t^{\mathrm{\prime }}<t})=0}$.${\displaystyle ◻}$

Theorem 1 can be readily extended to include the presence of other observed variables, ${\displaystyle {\mathit{V}}_{t^{\mathrm{\prime }}⩽t}}$, by redefining Transfer Entropy as ${\displaystyle T_{Y\to X}=I(X_t;{\mathit{Y}}_{t^{\mathrm{\prime }}<t}|{\mathit{X}}_{t^{\mathrm{\prime }}<t},{\mathit{V}}_{t^{\mathrm{\prime }}⩽t})}$, which discards contributions from indirect paths through other observed variables, ${\displaystyle {\mathit{V}}_{t^{\mathrm{\prime }}⩽t}}$.

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 Figure 2b. As a result, the presence of a temporal v-structure, ${\displaystyle Y_{t^{\mathrm{\prime }}}\to Y_t\leftarrow X_t}$ in ${\displaystyle {\mathcal{G}}_t}$, does not necessarily imply a vanishing Transfer Entropy, ${\displaystyle T_{Y\to X}=0}$, as long as there remains an edge between any ${\displaystyle Y_{t^{\mathrm{\prime }\mathrm{\prime }}}}$ and ${\displaystyle X_t}$, as in the example in Figure 2b. 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 contrast, tMIIC includes time-delayed as well as synchronous effects originating from unobserved latent variables, as discussed above.

The performance of CausalXtract’s causal discovery module (tMIIC) has been assessed using Tigramite package (Runge, 2020), 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 nonlinear relations between variables.

A first series of datasets was generated for a 15-node benchmark network (Figure 2—figure supplement 1a) with linear combinations of contributions inspired by the Tigramite package (Appendix 1). Running times and scores (Precision, Recall, Fscore) have been averaged over 10 datasets (Figure 2—figure supplement 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 (Figure 2—figure supplement 1b).

Importantly, increasing the number of time-lagged layers from ${\displaystyle \tau =2}$ (as in the actual model) to 5 or 10 layers in the inferred time-unfolded network (Figure 2—figure supplement 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 nonstationary variables, whilst the rest of the network remains relatively unaffected (Figure 2—figure supplement 3). It demonstrates that CausalXtract’s causal discovery module is robust to the presence of nonstationary variables but requires long-time range interactions, and therefore multiple time-lagged layers, to account for these nonstationary dynamics at specific variables. This striking observation on benchmark networks is also consistent with the multiple self-loops observed for a number of nonstationary variables in the real-world application on cellular ecosystems (Figure 3a, Figure 1—figure supplement 1).

A second series of more complex datasets was also generated for another 15-node benchmark network (Figure 2—figure supplement 4a) with nonlinear combinations of contributors (Appendix 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.

The source code of CausalXtract is available at https://github.com/miicTeam/CausalXtract, copy archived at miicTeam, 2024. It includes a demo R markdown notebook of CausalXtract pipeline, which reproduces step-by-step the results reported in the article (Figure 3), starting from the original live-cell time-lapse images of the tumor-on-chip ecosystem (Figure 1a). The Tigramite package used for benchmark comparison is available at https://github.com/jakobrunge/tigramite, copy archived at Runge, 2024.

• ${\displaystyle X_t^1\leftarrow -0.47f_2(X_{t-1}^1)+0.29f_3(X_{t-1}^2)\times {\eta }_1}$

• ${\displaystyle X_t^2\leftarrow 0.49f_2(X_{t-1}^2)+0.4f_1(X_{t-2}^1)+{\eta }_2}$

• ${\displaystyle X_t^3\leftarrow 0.56f_1(X_{t-1}^3)+0.44f_4(X_{t-2}^4)-0.26f_2(X_{t-2}^{10})+0.56f_2(X_t^4)+{\eta }_3}$

• ${\displaystyle X_t^4\leftarrow 0.24f_3(X_{t-1}^4)-0.24f_2(X_{t-2}^6)-0.12f_4(X_{t-1}^{14})\times {\eta }_4}$

• ${\displaystyle X_t^5\leftarrow -0.39f_3(X_{t-1}^5)-0.42f_3(X_{t-2}^5)-0.39f_3(X_t^{11})+{\eta }_5}$

• ${\displaystyle X_t^6\leftarrow -0.32f_2(X_{t-1}^6)+{\eta }_6}$

• ${\displaystyle X_t^7\leftarrow -0.17f_4(X_{t-1}^7)-0.17f_1(X_{t-2}^7)+{\eta }_7}$

• ${\displaystyle X_t^8\leftarrow 0.39f_4(X_{t-1}^8)-0.46f_4(X_{t-1}^7)-0.39f_3(X_{t-1}^1)-0.4f_3(X_{t-2}^{12})+{\eta }_8}$

• ${\displaystyle X_t^9\leftarrow -0.34f_1(X_{t-1}^9)+0.43f_3(X_{t-2}^{12})+{\eta }_9}$

• ${\displaystyle X_t^{10}\leftarrow 0.2f_1(X_{t-1}^{10})+0.18f_4(X_{t-2}^9)+0.17f_1(X_{t-1}^9)+0.48f_3(X_{t-1}^7)-0.26f_4(X_{t-1}^4)+{\eta }_{10}}$

• ${\displaystyle X_t^{11}\leftarrow 0.41f_2(X_{t-1}^{11})+0.54f_3(X_t^2)-0.55f_2(X_t^{12})+{\eta }_{11}}$

• ${\displaystyle X_t^{12}\leftarrow -0.45f_2(X_{t-1}^{12})-0.43f_4(X_{t-2}^3)-0.17f_4(X_{t-2}^9)\times {\eta }_{12}}$

• ${\displaystyle X_t^{13}\leftarrow 0.45f_3(X_{t-1}^{13})+{\eta }_{13}}$

• ${\displaystyle X_t^{14}\leftarrow 0.28f_2(X_{t-1}^{14})+0.37f_1(X_{t-2}^{12})\times {\eta }_{14}}$

• ${\displaystyle X_t^{15}\leftarrow 0.52f_3(X_{t-1}^{15})+{\eta }_{15}}$

• ${\displaystyle f_1(x)=x}$

• ${\displaystyle f_2(x)=x(1-4e^{-\frac{x^2}{2}})}$

• ${\displaystyle f_3(x)=x(1-4x^3{e}^{-\frac{x^2}{2}})}$

• ${\displaystyle f_4(x)=\cos (x)}$

The ${\displaystyle \eta }$ are white noises generated for each node or contribution using a normal distribution: ${\displaystyle \eta \sim \mathcal{N}(0,1)}$

• ${\displaystyle \begin{array}{l}{X}_t^1\leftarrow \eta -0.7f_6(u(\eta +X_{t-1}^1))-0.87f_5(u(\eta +(X_{t-1}^{14}\times X_{t-2}^1)))\end{array}}$

• ${\displaystyle \begin{array}{l}{X}_t^2\leftarrow \eta +0.65f_1(u(\eta +X_{t-1}^2))-0.63f_3(u(\eta +X_{t-2}^2))+0.79f_3(u(\eta +X_{t-1}^5))\end{array}}$

• ${\displaystyle \begin{array}{l}{X}_t^3\leftarrow \eta -0.76f_5(u(\eta +X_{t-1}^3))-0.59f_6(u(\eta +X_{t-1}^7))-0.85f_2(u(\eta +X_{t-1}^{15}))\\ -0.89f_5(u(\eta +(X_{t-2}^{13}\times X_{t-1}^7)))\end{array}}$

• ${\displaystyle \begin{array}{l}{X}_t^4\leftarrow \eta -0.7f_6(u(\eta +X_{t-1}^5))-0.86f_2(u(\eta +X_{t-2}^8))+0.53f_1(u(\eta +(X_{t-1}^4\times X_{t-2}^9)))\end{array}}$

• ${\displaystyle \begin{array}{l}{X}_t^5\leftarrow \eta +0.54f_2(u(\eta +(X_{t-1}^{14}\times X_{t-2}^6)))\end{array}}$

• ${\displaystyle \begin{array}{l}{X}_t^6\leftarrow \eta -0.85f_2(u(\eta +X_{t-1}^6))-0.79f_3(u(\eta +X_{t-2}^3))+0.59f_1(u(\eta +X_{t-1}^4))\\ +0.75f_3(u(\eta +X_t^1))+0.57f_2(u(\eta +X_{t-1}^{14}))\end{array}}$

• ${\displaystyle \begin{array}{l}{X}_t^7\leftarrow \eta +0.74f_1(u(\eta +X_{t-1}^7))+0.54f_6(u(\eta +X_{t-1}^9))-0.53f_2(u(\eta +(X_{t-1}^9\times X_{t-1}^7)))\end{array}}$

• $$\begin{gathered} {\displaystyle \begin{array}{l}{X}_t^8\leftarrow \eta \times (-0.63f_1(u(\eta +X_{t-1}^6))+0.81f_5(u(\eta +X_t^{13}))+0.53f_6(u(\eta +(X_{t-2}^6\times X_{t-1}^6))) \\ -0.69f_6(u(\eta +(X_t^{13}\times X_{t-1}^6))))\end{array}} \end{gathered}$$

• ${\displaystyle \begin{array}{l}{X}_t^9\leftarrow \eta +0.79f_3(u(\eta +X_{t-2}^4))+0.69f_6(u(\eta +(X_{t-1}^9\times X_{t-1}^{15})))\end{array}}$

• ${\displaystyle X_t^{10}\leftarrow \eta +0.54f_6(u(\eta +X_{t-1}^{10}))}$

• ${\displaystyle \begin{array}{l}{X}_t^{11}\leftarrow \eta +0.83f_6(u(\eta +X_{t-1}^{11}))-0.76f_4(u(\eta +X_{t-1}^{13}))-0.73f_3(u(\eta +X_{t-1}^2))\\ +0.74f_2(u(\eta +X_t^4))-0.87f_2(u(\eta +X_{t-2}^{10}))+0.72f_4(u(\eta +X_{t-1}^{12}))\\ -0.73f_1(u(\eta +(X_{t-2}^{10}\times X_{t-1}^{13})))\end{array}}$

• ${\displaystyle \begin{array}{l}{X}_t^{12}\leftarrow \eta +0.7f_3(u(\eta +X_{t-1}^{10}))-0.55f_5(u(\eta +X_t^9))-0.54f_5(u(\eta +(X_{t-1}^{12}\times X_{t-1}^{10})))\end{array}}$

• ${\displaystyle \begin{array}{l}{X}_t^{13}\leftarrow \eta -0.62f_3(u(\eta +X_{t-2}^{14}))-0.61f_1(u(\eta +(X_{t-1}^{13}\times X_{t-2}^{14})))\end{array}}$

• ${\displaystyle X_t^{14}\leftarrow \eta -0.78f_6(u(\eta +X_{t-1}^{14}))}$

• ${\displaystyle \begin{array}{l}{X}_t^{15}\leftarrow \eta -0.68f_4(u(\eta +X_{t-1}^{15}))+0.85f_4(u(\eta +X_{t-2}^{15}))-0.6f_5(u(\eta +X_{t-2}^{10}))\\ +0.68f_6(u(\eta +X_{t-1}^{14}))+0.81f_4(u(\eta +(X_{t-1}^{14}\times X_{t-2}^{10})))\end{array}}$

• ${\displaystyle \begin{array}{ccc}u(x)& =& max(-1,min(1,x))\end{array}}$

• ${\displaystyle \begin{array}{ccc}u(x)& =& max(-1,min(1,x))\end{array}}$

• ${\displaystyle \begin{array}{ccc}u(x)& =& max(-1,min(1,x))\end{array}}$

• ${\displaystyle \begin{array}{ccc}{f}_3(x)& =& 4x^2\end{array}}$

• ${\displaystyle \begin{array}{ccc}{f}_4(x)& =& 8x^3\end{array}}$

• ${\displaystyle \begin{array}{ccc}{f}_5(x)& =& 16x^4\end{array}}$

• ${\displaystyle \begin{array}{ccc}{f}_6(x)& =& \cos (\pi x)\end{array}}$

The ${\displaystyle \eta }$ are white noises generated for each node or contribution using a normal distribution: ${\displaystyle \eta \sim \mathcal{N}(0,0.1)}$.

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

This work was supported by ITMO Cancer (grant no. 20CM106) and the European Union’s Horizon 2020 research and innovation program 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 program.

1. Sent for peer review: February 5, 2024
2. Preprint posted: February 8, 2024
3. Reviewed Preprint version 1: September 17, 2024
4. Reviewed Preprint version 2: December 6, 2024
5. Version of Record published: January 17, 2025

You can cite all versions using the DOI https://doi.org/10.7554/eLife.95485. This DOI represents all versions, and will always resolve to the latest one.

© 2024, Simon, Comes, Tocci et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.