Continuity scaling: A rigorous framework for detecting and quantifying causality accurately

Identifying and ascertaining causal relations is a problem of paramount importance to science and engineering with broad applications bunge2017causality ; pearl2009causality ; Runge2019inferring . For example, accurate detection of causation is key to identifying the origin of diseases in precision medicine CV:2015 and is important to fields such as psychiatry SSFRLPWBBA:2016 . Traditionally, associational concepts are often misinterpreted as causation CoxHinkley1979 ; CoverThomas2012 , while causal analysis in fact goes one step further beyond association in a sense that, instead of using static conditions, causation is induced under changing conditions JPearl2009 . The principle of Granger causality formalizes a paradigmatic framework Wiener:1956 ; Granger:1969 ; Haufe2013 , quantifying causality in terms of prediction improvements, but, because of its linear, multivariate, and statistical regression nature, the various derived methods require extensive data ding2006granger . Entropy based methods Schreiber:2000 ; frenzel2007partial ; vicente2011transfer ; runge2012escaping ; SCB:2014 ; CSB:2015 ; STB:2015 ; research2021 face a similar difficulty. Another issue with the Granger causality is the fundamental requirement of separability of the underlying dynamical variables, which usually cannot be met in the real world systems. To overcome these difficulties, the cross-map-based techniques, paradigms tailored to dynamical systems, have been developed and have gained widespread attention in the past decade PRE_Hirata2010 ; PRE_Quiroga2000 ; PhyD_Arnhold1999 ; harnack2017topological ; SMYHDFM:2012 ; DFHKMMPYS:2013 ; WPCFMCHMPW:2014 ; MAC:2014 ; MW:2014 ; MAC:2014 ; YDGS:2015 ; clark2015spatial ; JHHLL:2016 ; ma2017detection ; ma2017detection ; Amigo2018 ; wang2018detecting ; pcm .

The cross map is originated from nonlinear time series analysis Takens:1981 ; PCFS:1980 ; sauer1991embedology ; stark1999delay ; stark2003delay ; muldoon1998delay . A brief understanding of such a map is as follows. Consider two subsystems: $X$ and $Y$. In the reconstructed phase space of $X$, if for any state vector at a time a set of neighboring vectors can be found, the set of the cross mapped vectors, which are the partners with equal time of $X$, could be available in $Y$. The cross map underlying the reconstructed spaces can be written as $Y_{t}=\Phi(X_{t})$ (where $X_{t}$ and $Y_{t}$ are delay coordinates with sufficiently large dimensions) for the case of $Y$ unidirectionally causing $X$ while, mathematically, its inverse map does not exist Amigo2018 . In practice, using the prior knowledge on the true causality in toy models or/and the assumption on the expanding property of $\Phi$ (representing by its Jacobian’s singular value larger than one in the topological causality framework harnack2017topological ), scientists developed many practically useful techniques based on the cross map for causality detection. For instance, the “activity” method, originally designed to measure the continuity of the inverse of the cross map, compares the divergence of the cross mapped vectors to the state vector in $X$ with the divergence of the independently-selected neighboring vectors to the same state vector PRE_Quiroga2000 ; PhyD_Arnhold1999 . The topological causality measures the divergence rate of the cross mapped vectors from the state vectors in $Y$ harnack2017topological , and the convergent cross mapping (CCM), increasing the length of time series, compares the true state vector $Y$ with the average of the cross mapped vectors, as the estimation of $Y$ SMYHDFM:2012 ; DFHKMMPYS:2013 ; WPCFMCHMPW:2014 ; MAC:2014 ; MW:2014 ; YDGS:2015 ; clark2015spatial ; JHHLL:2016 ; ma2017detection ; Amigo2018 ; wang2018detecting ; PRE_Hirata2010 ; pcm . Then, the change of the divergence or the accuracy of the estimation is statistically evaluated for determining the causation from $Y$ to $X$. Inversely, the causation from $X$ to $Y$ can be evaluated in an analogous manner. The above evaluations DFHKMMPYS:2013 ; WPCFMCHMPW:2014 ; MAC:2014 ; MW:2014 ; YDGS:2015 ; clark2015spatial ; JHHLL:2016 ; ma2017detection ; harnack2017topological ; Amigo2018 ; wang2018detecting ; PRE_Hirata2010 ; pcm can be understood at a conceptional and qualitative level and perform well in many demonstrations.

In this work, striving for a comprehensive understanding of causal mechanisms and inspired by the cross-map-based techniques, we develop a mathematically rigorous framework for detecting causality in nonlinear dynamical systems, turning eyes towards investigating the original systems from their cross maps, which is also logically consistent with the natural interpretation of causality as functional dependences pearl2009causality ; JPearl2009 . The skills used in cross-map-based methods are assmilated in our framework, while we directly study the original dynamical systems or the reconstructed systems instead of the cross maps. The foundation of our framework is the scaling law for the changing relation of $\varepsilon$ with $\delta$ arising from the continuity for the investigated system, henceforth the term “continuity scaling”. In addition to providing a theory, we demonstrate, using synthetic and real-world data, that our continuity scaling framework is accurate, computationally efficient, widely applicable, showing advantages over the existing methods.

To explain the mathematical idea behind the development of our framework, we use the following class of discrete time dynamical systems: $\bm{x}_{t+1}=\bm{f}(\bm{x}_{t},\bm{y}_{t})$ and $\bm{y}_{t+1}=\bm{g}(\bm{x}_{t},\bm{y}_{t})$ for $t\in\mathbb{N}$, where the state variables $\{\bm{x}_{t}\}_{t\in\mathbb{N}}$, $\{\bm{y}_{t}\}_{t\in\mathbb{N}}$ evolve in the compact manifolds $\mathcal{M}$, $\mathcal{N}$ of dimension $D_{\mathcal{M}}$, $D_{\mathcal{N}}$ under sufficiently smooth map $\bm{f}$, $\bm{g}$, respectively. We adopt the common recognition of causality in dynamical systems.

We now interpret the causal relationship stipulated by the continuity of a function. Let $\bm{f}_{\bm{x}_{\text{g}}}(\cdot)\triangleq\bm{f}(\bm{x}_{\text{g}},\cdot)$ for a given point $\bm{x}_{\text{g}}\in\mathcal{M}$. For any $\bm{y}_{P}\in\mathcal{N}$, we denote its image under the given function by $\bm{x}_{I}\triangleq\bm{f}_{\bm{x}_{\text{g}}}(\bm{y}_{P})$. Applying the logic statement of a continuous function to $\bm{f}_{\bm{x}_{\text{g}}}(\cdot)$, we have that, for any neighborhood $\mathcal{O}(\bm{x}_{I},\varepsilon)$ centered at $\bm{x}_{I}$ and of radius $\varepsilon>0$, there exists a neighborhood $\mathcal{O}(\bm{y}_{P},\delta)$ centered at $\bm{y}_{P}$ of radius $\delta>0$, such that $\bm{f}_{\bm{x}_{\text{g}}}(\mathcal{O}(\bm{y}_{P},\delta))\subset\mathcal{O}(\bm{x}_{I},\varepsilon)$. The neighborhood and its radius are defined by

where $\text{dist}_{\mathcal{S}}(\cdot,\cdot)$ represents an appropriate metric describing the distance between two given points in a specified manifold $\mathcal{S}$ with $\mathcal{S}=\mathcal{M}$ or $\mathcal{N}$. The meaning of this mathematical statement is that, if we have a neighborhood of the resulting variable $\bm{x}_{I}$ first, we can then find a neighborhood for the causal variable $\bm{y}_{P}$ to satisfy the above mapping and inclusion relation. This operation of “first-$\varepsilon$-then-$\delta$” provides a rigorous base for the principle that the information about the resulting variable can be used to estimate the information of the causal variable and therefore to ascertain causation, as indicated by the long arrow in the middle panels of Fig. 1. Note that, the existence of the $\delta>0$ neighborhood is always guaranteed for a continuous map $\bm{f}_{\bm{x}_{\text{g}}}$. In fact, due to the compactness of the manifold $\mathcal{N}$, a largest value of $\delta$ exists. However, if $\bm{y}_{P}$ does not have an explicit causal effect on the variable $\bm{x}_{I}$, i.e., $\bm{f}_{\bm{x}_{\rm g}}$ is independent of $\bm{y}_{P}$, the existence of $\delta$ is still assured but it is independent of the value of $\varepsilon$, as shown in the upper panel of Fig. 1. This means that merely determining the existence of a $\delta$-neighborhood is not enough for inferring causation - it is necessary to vary $\varepsilon$ systematically and to examine the scaling relation between $\delta$ and $\varepsilon$. In the following we discuss a number of scenarios.

Case I: Dynamical variables $\{(\bm{x}_{t},\bm{y}_{t})\}_{t\in\mathbb{N}}$ are fully measurable. For any given constant $\varepsilon_{\bm{x}}>0$, the set $\{\bm{x}_{\tau}\in\mathcal{M}~{}|~{}\tau\in I^{t}_{\bm{x}}(\varepsilon_{\bm{x}})\}$ can be used to approximate the neighborhood $\mathcal{O}(\bm{x}_{t+1},\varepsilon_{\bm{x}})$, where the time index set is

The radius $\delta^{t}_{\bm{y}}=\delta^{t}_{\bm{y}}(\varepsilon_{\bm{x}})$ of the neighborhood $\mathcal{O}(\bm{y}_{t},\delta^{t}_{\bm{y}})$ satisfying $\bm{f}_{\bm{x}_{\rm g}=\bm{x}_{t}}(\mathcal{O}(\bm{y}_{t},\delta^{t}_{\bm{y}}))\subset\mathcal{O}(\bm{x}_{t+1},\varepsilon_{\bm{x}})$ can be estimated as

where $\#[\cdot]$ is the cardinality of the given set and the index set is $\bar{I}^{t}_{\bm{x}}(\varepsilon_{\bm{x}})\triangleq\{\tau\in I^{t}_{\bm{x}}(\varepsilon_{\bm{x}})~{}|~{}~{}\text{dist}_{\mathcal{M}}(\bm{x}_{t},\bm{x}_{\tau-1})<\varepsilon_{\bm{x}}\}$.

The strict mathematical steps for estimating $\delta^{t}_{\bm{y}}$ are given in Section II of Supplementary Information (SI). We emphasize that here correspondence between $\bm{x}_{t+1}$ and $\bm{y}_{t}$ is investigated, differing from the cross-map-based methods, with one-step time difference naturally arising. This consideration yields a key condition [DD], which is only need when considering the original iteration/flow and whose detailed description and universality are demonstrated in SI. We reveal a linear scaling law between $\langle\delta^{t}_{\bm{y}}\rangle_{t\in\mathbb{N}}$ and $\ln{\varepsilon_{\bm{x}}}$, as shown in the lower panels of Fig. 1, whose slope $s_{\bm{y}\hookrightarrow\bm{x}}$ is an indicator of the correspondent relation between $\varepsilon$ and $\delta$ and hence the causal relation $\bm{y}\hookrightarrow\bm{x}$. Here, $\langle\cdot\rangle_{t\in\mathbb{N}}$ denotes the average over time. In particular, a larger slope value implies a stronger causation in the direction from $\bm{y}$ to $\bm{x}$ as represented by the map functions $\bm{f}(\bm{x}_{t},\bm{y}_{t})$ [Fig. 1], while a near zero slope indicates null causation in this direction [Fig. 1]. Likewise, possible causation in the reversed direction, $\bm{x}\hookrightarrow\bm{y}$, as represented by the function $\bm{g}(\bm{x}_{t},\bm{y}_{t})$, can be assessed analogously. And the unidirectional case when $\bm{f}(\bm{x},\bm{y})=\bm{f}_{0}(\bm{x})$ independent of $\bm{y}$ is uniformly considered in Case II. We summarize the consideration below and an argument for the generic existence of the scaling law is provided in Section II of SI.

Does the causation from $\bm{v}$ to $\bm{u}$ imply causation from $\bm{y}$ to $\bm{x}$? The answer is affirmative, which can be argued, as follows. If the original map function $\bm{f}$ is independent of $\bm{y}$: $\bm{f}(\bm{x},\bm{y})=\bm{f}_{0}(\bm{x})$, there is no causation from $\bm{y}$ to $\bm{x}$. In this case, the embedding $\bm{E}_{u}(\bm{x},\bm{y})$ becomes independent of $\bm{y}$, degenerating into the form of $\bm{E}_{u}(\bm{x},\bm{y})=\bm{E}_{u0}(\bm{x})$, a diffeomorphism from $\mathcal{M}$ to $\widetilde{\mathcal{L}}_{u0}=\bm{E}_{u0}(\mathcal{M})$ only. As a result, Eq. (3) becomes: $\bm{u}_{t+1}=\tilde{\bm{f}}_{0}(\bm{u}_{t})$ and $\bm{v}_{t+1}=\tilde{\bm{g}}(\bm{u}_{t},\bm{v}_{t})$, where $\tilde{\bm{f}}_{0}(\bm{u})=\bm{E}_{u0}\circ\bm{f}\circ\bm{E}_{u0}^{-1}(\bm{u})$ and the resulting mapping $\tilde{\bm{f}}_{0}$ is independent of $\bm{v}$. The independence can be validated by computing the slope $s_{\bm{v}\hookrightarrow\bm{u}}$ associated with the scaling relation between $\langle\delta^{t}_{\bm{v}}\rangle_{t\in\mathbb{N}}$ and $\ln{\varepsilon_{\bm{u}}}$, where a zero slope indicates null causation from $\bm{v}$ to $\bm{u}$ and hence null causation from $\bm{y}$ to $\bm{x}$. Conversely, a finite slope signifies causation between the variables. Thus, any type of causal relation (unidirectional or bi-directional) detected between the reconstructed state variables $\{(\bm{u}_{t},\bm{v}_{t})\}_{t\in\mathbb{N}}$ implies the same type of causal relation between the internal but inaccessible variables $\bm{x}$ and $\bm{y}$ of the original system.

Case III: The structure of the internal variables is completely unknown. Given the observational functions $\tilde{u},\tilde{v}$: $\mathcal{M}$$\times$$\mathcal{N}$$\to$$\mathbb{R}$ with $\tilde{u}_{t}=\tilde{u}(\bm{x}_{t},\bm{y}_{t})$ and $\tilde{v}_{t}=\tilde{v}(\bm{x}_{t},\bm{y}_{t})$, we first reconstruct the state space: $\tilde{\bm{u}}_{t}=(\tilde{u}_{t},\tilde{u}_{t+\tau},\cdots,\tilde{u}_{t+(d-1)\tau})^{\top}$ and $\tilde{\bm{v}}_{t}=(\tilde{v}_{t},\tilde{v}_{t+\tau},\cdots,\tilde{v}_{t+(d-1)\tau})^{\top}$. To detect and quantify causation from $\tilde{\bm{v}}$ to $\tilde{\bm{u}}$ (or vice versa), we carry out a continuity scaling analysis with the modified indices $I^{t}_{\tilde{\bm{u}}}(\varepsilon_{\tilde{\bm{u}}})$, $\delta^{t}_{\tilde{\bm{v}}}(\varepsilon_{\tilde{\bm{u}}})$ and $s_{\tilde{\bm{v}}\hookrightarrow\tilde{\bm{u}}}$. Differing from Case II, here, due to the lack of knowledge about the correspondence structure between the internal and observational variables, a causal relation for the latter does not definitely imply the same for the former.

Case IV: Continuous-time dynamical systems possessing a sufficiently smooth flow $\{\bm{S}_{t};t\in\mathbb{R}\}$ on a compact manifold $\mathcal{H}$: $\text{d}\bm{S}_{t}(\bm{u}_{0})/\text{d}t=\bm{\chi}(\bm{S}_{t}(\bm{u}_{0}))$, where $\bm{\chi}$ is the vector field. Let $\{\hat{u}_{t=\omega n+\nu}\}_{n\in\mathbb{Z}}$ and $\{\hat{v}_{t=\omega n+\nu}\}_{n\in\mathbb{Z}}$ be two respective time series from the smooth observational functions $\hat{u},~{}\hat{v}$: $\mathcal{H}$$\rightarrow$$\mathbb{R}$ with $\hat{u}_{t}=\hat{u}(\bm{S}_{t})$ and $\hat{v}_{t}=\hat{v}(\bm{S}_{t})$, where $1/\omega$ is the sampling rate and $\nu$ is the time shift. Defining $\bm{\Xi}\triangleq\bm{S}_{\omega}$: $\mathcal{H}$$\rightarrow$$\mathcal{H}$ and $\hat{\bm{S}}_{n}\triangleq{\bm{S}_{\omega n+\nu}}(\bm{u}_{0})$, we obtain a discrete-time system as $\hat{\bm{S}}_{n+1}=\bm{\Xi}(\hat{\bm{S}}_{n})$ with the observational functions as $\hat{u}_{n}=\hat{u}(\hat{\bm{S}}_{n})$ and $\hat{v}_{n}=\hat{v}(\hat{\bm{S}}_{n})$, reducing the case to Case III and rendering applicable our continuity scaling analysis to unveil and quantify the causal relation between $\{\hat{u}_{t=\omega n+\nu}\}_{n\in\mathbb{Z}}$ and $\{\hat{v}_{t=\omega n+\nu}\}_{n\in\mathbb{Z}}$. If the domains of $\hat{u}$ and $\hat{v}$ have their own restrictions on some particular subspaces, e.g., $\hat{u}$: $\mathcal{H}_{u}$$\rightarrow$$\mathbb{R}$ and $\hat{v}$: $\mathcal{H}_{v}$$\rightarrow$$\mathbb{R}$ with $\mathcal{H}=\mathcal{H}_{u}\oplus\mathcal{H}_{v}$, the case is further reduced to Case II, so the detected causal relation between the observational variables imply causation between the internal variables belonging to their respective subspaces.

To demonstrate the efficacy of our continuity scaling framework and its superior performance, we have carried out extensive numerical tests with a large number of synthetic and empirical datasets, including those from gene regulatory networks as well as those of air pollution and hospital admission. The practical steps of the continuity scaling framework together with the significance test procedures are described in Methods. We present three representative examples here, while leaving others of significance to SI.

The first example is an ecological model of two unidirectionally interacting species: $x_{1,t+1}=x_{1,t}(3.8-3.8x_{1,t}-\mu_{12}x_{2,t})$ and $x_{2,t+1}=x_{2,t}(3.7-3.7x_{2,t}-\mu_{21}x_{1,t})$. With time series $\{(x_{1,t},x_{2,t})\}_{t\in\mathbb{N}}$ obtained from different values of the coupling parameters, our continuity scaling framework yields correct results of different degree of unidirectional causation, as shown in Figs. 3-3. In all cases, there exists a reasonable range of $\ln{\varepsilon_{x_{2}}}$ (neither too small nor too large) from which the slope $s_{x_{1}\hookrightarrow x_{2}}$ of the linear scaling can be extracted. The statistical significance of the estimated slope values and consequently the strength of causation can be assessed with the standard $p$-value test lancaster2018surrogate (Methods and SI). An ecological model with bidirectional coupling has also been tested (see Section III of SI). Figures 3-3 show the results from ecological networks of five mutually interacting species on a ring and on a tree structure, respectively, where the color-coded slope values reflect accurately the interaction patterns in both cases.

The second example is the coupled Lorenz system: $\dot{x}_{i}=\sigma_{i}(y_{i}-x_{i})+\mu_{ij}x_{j}$, $\dot{y}_{i}=x_{i}(\rho_{i}-z_{i})-y_{i}$, $\dot{z}_{i}=x_{i}y_{i}-\beta_{i}z_{i}$ with $i,j=1,2$ and $i\not=j$. We use time series $\{y_{1,t},y_{2,t}\}_{t=n\omega}$ for detecting different configurations of causation (see Section III of SI). Figure 4 presents the overall result, where the color-coded estimated values of the slope from the continuity scaling are shown for different combinations of the sampling rate $1/\omega$ and coupling strength. Even with relatively low sampling rate, our continuity scaling framework can successfully detect and quantify the strength of causation. Note that the accuracy does not vary monotonously with the sampling rate, indicating the potential of our framework to ascertain and quantify causation even with rare data. Moreover, the proposed index can accurately reflect the true causal strength (denoting by the coupling parameter), which is also evidenced by numerical tests in Sections III and IV of SI. Robustness tests against different noise perturbations are provided in Section III of SI demonstrating the practical usefulness of our framework. Additionally, analogous to the first example, we present in SI several examples on causation detection in the coupled Lorenz system with nonlinear couplings, and the Rössler-Lorenz system, etc., which further demonstrates the generic efficacy of our framework.

In addition, we present study on several real-world dataset, which brings new insights to the evolutionary mechanism of underlying systems. We study gene expression data from DREAM4 in silico Network Challenge marbach2009generating ; marbach2010revealing , whose intrinsic gene regulatory networks (GRNs) are known for verification [Fig. 5 and Fig. S17 of SI]. Applying our framework to these data, we ascertain the causations between each pair of genes by using the continuity scaling framework. The corresponding ROC curves for five different networks as well as their AUROC values are shown in Fig. 5, which indicates a high detection accuracy in dealing with real-world data.

We then test the causal relationship in a marine ecosystem consisting of Pacific sardine landings, northern anchovy landings and sea surfae temperature (SST). We reveal new findings to support the competing relationship hypothesis stated in Lasker1983Ocean which cannot be detected by CCM SMYHDFM:2012 . As pointed out in Fig. 6, while common influence from SST to both species is verified with both methods, our continuity scaling additionally illuminates notable influence from anchovy to sardine with its reverse direction being less significant. While competing relationship plays an important role in ecosystems, continuity scaling can reveal more essential interaction mechanism. See Section III.E of SI for more details.

Moreover, we study the transmission mechanism of the recent COVID-19 pandemic. Particularly, we analyze the daily new cases of COVID-19 of representative contries for two stages: day $1$ (January 22${}^{\text{nd}}$ 2020) to day $100$ (April 30${}^{\text{th}}$ 2020) and day $101$ (May 1${}^{\text{st}}$ 2020) to day $391$ (February 15${}^{\text{th}}$ 2021). Our continuity scaling is pairwisely applied to reconstruct the transmission causal network. As shown in Fig. 7, China shows a significant effect on a few countries at the first stage and this effect disappears at the second stage. However, other countries show a different situation with China, whose external effect lasts as shown in Section III.E and Fig. S18 of SI. Our results accord with that China holds stringent epidemic control strategies with sporadic domestic infections, as evidenced by official daily briefings, demonstrating the potential of continuity scaling in detecting causal networks for ongoing complex systems. Additionally, We emphasize that day $100$ is a suitable critical day to distinguish the early severe stage and the late well-under-control stage of the pandemic [see Fig. S18(a) of SI], while slight change of the critical day will not nullify our result. As shown in Fig. S18(b) of SI, when the critical day varies from day $94$ to day $106$, no significant change (less than $5\%$) of the detected causal links occurs at both stages, and the number of countries under influence of China at Stage $2$ remains zero. See more details in Section III.E of SI.

Additional real world examples including air pollutants and hospital admission record from Hong Kong are also shown in Section III of SI.

To summarize, we have developed a novel framework for data based detection and quantification of causation in complex dynamical systems. On the basis of the widely used cross-map-based techniques, our framework enjoys a rigorous foundation, focusing on the continuity scaling law of the concerned system directly instead of only investigating the continuity of its cross map. Therefore, our framework is consistent with the standard interpretation of causality, and works even in the typical cases where several existing typical methods do not perform that well or even they fail (see the comparison results in Section IV of SI). In addition, the mathematical reasoning leading to the core of our framework, the continuity scaling, helps resolve the long-standing issue associated with techniques directly using cross map that information about the resulting variables is required to project the dynamical behavior of the causal variables, whereas several works in the literature Quyen1999PhyD , which directly studied the continuity or the smoothness of the cross map, likely yielded confused detected results on causal directions.

Computational complexity. The computational complexity of the algorithm is $O(T^{2}N_{\varepsilon})$, which is relatively smaller than the CCM method, whose computational complexity is $O(T^{2}\log T)$.

Limitations and future works. Nevertheless, there are still some limitations in the presently proposed framework. For instance, currently, only bivariate detection algorithm is designed, so generalization to multivariate network inference requires further considerations, as analogous to those works presented in Refs. Peters2017elements ; Runge2018causal ; research2019 . In addition, the causal time delay has not been taken into account in the current framework, so it also could be further investigated, similar to the work reported in Ref. ma2017detection . Definitely, we will settle these questions in our future work.

Detecting causality in complex dynamical systems has broad applications not only in science and engineering, but also in many aspects of the modern society, demanding accurate, efficient, and rigorously justified and hence trustworthy methodologies. Our present work provides a vehicle along this feat, and indeed resolves the puzzles arising in the use of those influential methods.

Continuity scaling framework: a detailed description of algorithms. Let $\{u_{t}\}_{t=1,2,\dots,T}$ and $\{v_{t}\}_{t=1,2,\dots,T}$ be two experimentally measured time series of internal variables $\{(\bm{x}_{t},\bm{y}_{t})\}_{t\in\mathbb{N}}$. Typically, if the dynamical variables $\{(\bm{x}_{t},\bm{y}_{t})\}_{t\in\mathbb{N}}$ are accessible, $\{(u_{t},v_{t})\}$ reduce to one-dimensional coordinate of the internal system. The key computational steps of our continuity scaling framework are described, as follows.

We present the steps for causation detection using the case of $\bm{v}\hookrightarrow\bm{u}$ as an example.

As defined, $\langle\delta^{t}_{\bm{v}}(\varepsilon_{\bm{u}})\rangle_{t\in\mathbb{N}}$ is the average of $\delta^{t}_{\bm{v}}(\varepsilon_{\bm{u}})$ over all possible time $t$. We use a finite number of pairs $\{(\langle\delta^{t}_{\bm{v}}(\varepsilon_{\bm{u},j})\rangle_{t\in\mathbb{N}_{T_{0}}},\ln\varepsilon_{\bm{u},j})\}_{j=1,\dots,N_{\varepsilon}}$ to approximate the scaling relation between $\langle\delta^{t}_{\bm{v}}(\varepsilon_{\bm{u}})\rangle_{t\in\mathbb{N}}$ and $\ln\varepsilon_{\bm{u}}$, where $\mathbb{N}_{T_{0}}=\{1,2,\cdots,T_{0}\}$. Theoretically, a larger value of $N_{\epsilon}$ and a smaller value of $e$ will result in a more accurate approximation of the scaling relation. In practice, the accuracy is determined by the length of the observational time series, the sampling duration, and different types of noise perturbations. In numerical simulations, we set $e=0.001$ and $N_{\epsilon}=33$. In addition, a too large or a too small value of $\varepsilon_{\bm{u}}$ can induce insufficient data to restore the neighborhood and/or the entire manifold. We thus set $\delta^{t}_{\bm{v}}(\varepsilon_{\bm{u},j})=\delta^{t}_{\bm{v}}(\varepsilon_{\bm{u},j+1})$ as a practical technique as the number of points is limited practically in a small neighborhood. As a result, near zero slope values would appear on both sides of the scaling curve $\langle\delta^{t}_{\bm{v}}(\varepsilon_{\bm{u}})\rangle_{t\in\mathbb{N}}$-$\ln\varepsilon_{\bm{u}}$, as demonstrated in Fig. 3 and in SI. In such a case, to estimate the slope of the scaling relation, we take the following approach.

For the other causal direction from $\bm{u}$ to $\bm{v}$, these steps are equally applicable to estimating the slope $s_{\bm{u}\hookrightarrow\bm{v}}$.

To assess the statistical significance of the numerically determined causation, we devise the following surrogate test using the case of $\bm{v}$ causing $\bm{u}$ as an illustrative example.

In simulations, we set the number of segments to be $N_{G}=25$ and the number of times for random permutations to be $Q\geq 20$.

W.L. is supported by the National Key R&D Program of China (Grant No. 2018YFC0116600), by the National Natural Science Foundation of China (Grant Nos. 11925103 & 61773125), by the STCSM (Grant No. 18DZ1201000), and by the Shanghai Municipal Science and Technology Major Project (No. 2021SHZDZX0103). Y.-C.L. is supported by AFOSR (Grant No. FA9550-21-1-0438). S.-Y.L. is supported by the National Natural Science Foundation of China (No. 12101133) and “Chenguang Program” supported by Shanghai Education Development Foundation and Shanghai Municipal Education Commission (No. 20CG01). Q.N. is partially supported by NSF (Grant No. DMS1763272) and the Simons Foundation (Grant No. 594598).

The authors declare no competing interests.

The source codes for our CS framework are available at https://github.com/bianzhiyu/ContinuityScaling.

W.L. conceived idea; X.Y., S.-Y.L., and W.L. designed and performed research; X.Y., S.-Y.L., H.-F.M., and W.L. analyzed data, H.-F.M., Y.-C.L., and Q.N. contributed data and analysis tools, and all the authors wrote the paper. X.Y. and S.-Y.L. were equally contributed.