Mathematically, an SCM for variables {V ₁, …, V _d } is a collection of so-called “structural assignments” of the form V j : = f j P A j , N j , j = 1 , … , d . (1) where the right hand side function f _j determines the assignment of the value of V _j on the left-hand side based on the values of parents variables PA _j ⊂ {V ₁, …, V _d } and of the so-called exogenous random variable N _j . The SCM is associated to a directed graph, the causal graph, where each variable V _j is represented by a node, and the parent-child relations between them is indicated by a “parent → child” arrow. While SCMs do not necessary include time information, we can exploit them to study dynamic interactions between two subsystems by considering causal relations linking variables representing one subsystem at past time points, to variables representing the other subsystem at a future time point. As an example, Figure 1A shows two uni-directionally coupled brain regions whose activities are represented by time series { … , X t 1 , X t + 1 1 , … } and { … , X t 2 , X t + 1 2 , … } and the corresponding SCM links the past of X ² to the future of X ¹. Typically, such a model also includes dependencies of each regional activity on its own past, and those dependencies can involve multiple time steps, leading to causal graphs of the form exemplified in Figure 2A. We focus on information theoretic causality measures, which are typically “model free”, in the sense that they can be expressed independently from the choice of functions f _j in the assignments of Eq. 1. However, model-free estimation of information theoretic quantities is a challenging problem that we will not address in this paper. Instead, estimation of the relevant quantities will rely on the following linear time-inhomogeneous Structural Vector Autoregressive (SVAR) model. X t 1 : = a t ⊤ X p , t 1 + b t ⊤ X p , t 2 + η t 1 , η t 1 ∼ N k t 1 , σ 1 , t 2 , (2a) X t 2 : = c t ⊤ X p , t 1 + d t ⊤ X p , t 2 + η t 2 , η t 2 ∼ N k t 2 , σ 2 , t 2 . (2b) where:

• η t k is the innovation for channel k at time t, sampled independently from the past values of X and from innovations at other time points and/or channels,

Note that the constraint of independence between exogenous variables η t k associated to different channels entails the assumption of no contemporaneous effects, in contrast to, e.g., Moneta et al. (2011). Figure 2A exemplifies the associated causal graph for p = 2 and c _t = 0, such that X ² is dependent only on X ₂ itself, encodes unidirectional causation from X ² to X ¹. The methods presented in this paper can be applied to time series generated by any other Markovian time series model (e.g., non-linear models (Marinazzo et al., 2011a)). However, the choice of SVAR 1) allows expressing the link between Granger causality and Transfer Entropy, 2) facilitates the estimation of all information theoretic quantities (which is inherently hard in a model-free setting, see, e.g., McAllester and Stratos (2020), as they get an analytic expression based on second order statistics), 3) avoids issues related to the non-parametric estimation of information theoretic quantities (e.g., finite sample bias), 4) allows to build easily interpretable models of transient oscillations.

One key question in causality is estimating the effect of (possibly imaginary) manipulations (see, e.g., Woodward, 2001) of the system of interest from data, which boils down to comparing two “worlds” or scenarios (Shpitser and Pearl, 2008): the original world where no manipulation is performed, and the “post-manipulation” world.

Both original and post-manipulation worlds typically cannot be measured simultaneously (e.g., “treatment” and “no treatment” in the same patient). However, estimating their differences arguably forms the basis of causal investigations in empirical sciences, for example, by performing randomized experiments on multiple instances of a system designed with mutually exclusive treatments to infer the outcome of manipulations of this system. However, even performing carefully controlled experiments on close to identical instances of a system is often challenging in reality, as many physical and physiological phenomena cannot be easily reproduced or manipulated. This is typically the case for spontaneous transient neural events investigated in this paper, where neurophysiological experimental techniques limit the understanding and control of their conditions of occurrence, as well as the ability to precisely modify some aspects of network activity to test assumptions on the underlying mechanisms.

Under additional model assumptions, such as the absence of unobserved confounders, the framework of SCMs (as briefly introduced in Section 2.1.1 and Supplementary Section SA), can be leveraged to infer the outcome of manipulations based on observational data only. Assuming those assumptions are met (see also Discussion for examples), the SCM inferred from data can be modified using a family of operations named interventions to model the manipulations of the system described by the SCM (Pearl, 2000; Peters et al., 2017). Intervening typically refers to modifying the structural equation associated to one node in the SCM, to study the modifications it brings about in the system. When interventions are performed, the only affected mechanistic relations (represented by arrows in an SCM) are the ones between the intervened nodes and their parent nodes. For instance, one can impose a fixed deterministic value to a node, or that this node’s variable is drawn from a given distribution, independently from other variables in the SCM (Janzing et al., 2013; Correa and Bareinboim, 2020; Peters et al., 2017, Chapter 3). Both such interventions lead to an intervened causal graph where the arrows between the node intervened upon and its parents are removed. For a better understanding of the manipuation modelled in these two cases, consider examples from experimental manipulations in electrophysiology. The first case involves fixing a value, as in voltage clamp techniques used to study channel conductance. By fixing the membrane potential of a single neuron, the causal arrow between the membrane potential and extracellular ion concentration is broken. An example of the second case is injecting current in plasticity studies to maintain certain firing rates in a patch-clamped neuron, ensuring that the pharmacological shutdown of certain ion channels does not cause a change in neuronal firing patterns.

Importantly, while an intervention modifies the graph associated to an SCM, the variables’ joint distribution can still be obtained by exploiting interventional knowledge, (unintervened) observational data and prior assumptions related to the unaffected conditionals.

For the SCM introduced in Eq. 1, intervening on V _k consists in replacing its structural assignment by a new one: V k : = f k ̃ P A k ̃ , N k ̃ . (3) where the new function f k ̃ , set of parents P A k ̃ , and/or the distribution of the exogenous variable N k ̃ may be differ from the original ones. The resulting distribution P ̃ V is called intervention distribution and denoted P V  do ( V k : = f k ̃ ( P A k ̃ , N k ̃ ) ) (see e.g., Peters et al. (2017, chapter 6)), where the superscript indicates that we refer to the distribution resulting from the modification of the SCM’s kth equation by Eq. 3, which is called a “do operation”. Meanwhile, the other structural equations and the distribution of their associated exogenous variables are kept unchanged.

Typically, one brain region influences another through axonal propagation of spiking activity in afferent neurons, and synaptic transmission to dendrites of the target region. In the simplest bi-variate case (i.e., if we focus on two brain regions with direct synaptic projections), this relationship can be represented in an SCM by Figure 1A(left) where two nodes representing the neural activities of the two regions are linked by a uni-directional arrow. For the sake of simplicity, we ignore for now the influence of the past of X ² on itself, such that the only causal link is from the past (at t − 1) of region 2 to the present (at t) of region 1. The joint probability can be causally factorized as P ( X t − 1 2 , X t 1 ) = P ( X t 1 | X t − 1 2 ) P ( X t − 1 2 ) , where P ( X t 1 | X t − 1 2 ) reflects the stochastic map or causal mechanism between the child and parent of the arrow: for example, Figure 1A(bottom left) shows a typical example how a nonlinear structural equation (i.e., the arrow) induces dependencies between a normally-distributed X t − 1 2 and X t 1 .

A natural manipulation to study the causal mechanism is to shut down the synaptic transmission during event occurrence and compare the outcomes, e.g., via physically cutting the pre-synaptic dendrite attached to the synapse, or pharmacological blockade of the relevant ion channel. Suppose that the experiment could be done, the data obtained in this hypothetical manipulated scenario can be modeled by another SCM without the arrow between X t − 1 2 and X t 1 , as seen in Figure 1A(right middle), such that they are independent from each other with P ̃ ( X t − 1 2 , X t 1 ) = P ̃ ( X t 1 ) P ̃ ( X t − 1 2 ) due to the Markov properties (Supplementary Section SA). In a model-free setting and in absence of data where this (technically challenging) manipulation is actually done, the choice of factorized P ̃ , such as the marginal mean and covariances, is non-trivial (see illustration in Figure 1A(bottom right)). Janzing et al. (2013) introduce a well-grounded way to make this choice to emulate this experimental ablation of connectivity, that generalize to arbitrary causal graphs.

This approach will be thoroughly discussed in Section 2.2.3 but here we explain it briefly in this simplified example to provide the readers with an intuitive understanding of the principles. Figure 1B illustrates the intervention performed in the SCM: cut the causal arrow from X t − 1 2 to X t 1 and feed X t 1 with an independent copy of X t − 1 2 (denoted X t − 1 2 ′ ), where independent copy means that X t − 1 2 ′ is a random variable statistically independent from any exogenous variables within the graph, and with the same marginal distribution as X t − 1 2 . The idea of this intervention is to achieve the independence of variables in the post-manipulation world by exploiting the observational conditional and marginal probabilities available to us. That is:

• choose P ̃ ( X t − 1 2 ) to be the observed marginal distribution of the cause, i.e., P ̃ ( X t − 1 2 ) = P ( X t − 1 2 ) , because it is the only one available to us in a model-free setting. Consider alternative choices: we could set the cause to a constant, but which constant to choose is unclear without additional information on the system. For example, even taking the observational average would not be realistic if X t − 1 2 is binary,

The strength of the causal arrow is then quantified by the Kullback-Leibler (KL) divergence D _KL between original (unintervened) and intervened joint densities D K L ( P ‖ P ̃ ) = D K L ( P ( X t − 1 2 , X t 1 ) ‖ P ( X t − 1 2 ) P ( X t 1 ) ) , which in this simple case thus boils down to the mutual information between the two nodes. Figure 1B(bottom) illustrates the contrast between the actual joint distribution and the intervention distribution for such intervention: by comparing these two distributions one could quantify how much the causal mechanism tilts the Gaussian shape, while they still largely overlap. In a context where nodes correspond to single neurons, this can be thought of as a proxy for the experiment of cutting the axon of afferent neurons, while injecting a current to maintain the baseline level of excitation in the target neuron, such that it is kept in naturalistic conditions.

However, we will argue that this choice of intervention is not an ideal way to measure the event-based causal influences between two brain regions. Going back to the manipulation experiment, despite the transient activities occurring at the pre-synaptic neurons as an input to the synapse, due to the cutoff of afferents or blockage of ion channels, the activity of the post-synaptic neuron (i.e., the effect variable) is expected to remain at baseline level without being influenced by the event occurring in the cause. In this context, the operation of feeding the effect node X t 1 with an independent copy of the cause node X t − 1 2 at the same time t − 1 still implicitly incorporates the influence of the event-related transient changes undergone by X ² at the time t − 1 on X t 1 , as the distribution of X t − 1 2 may strongly differ from what it is during baseline activity (without the occurrence of events). Therefore we propose instead to reconstruct the baseline state by replacing the independent copy of the marginal distribution of the event-related activity X t − 1 2 by the marginal distribution of a baseline state in X t ref 2 , at a reference time point t _ref where the event of interest has not yet occurred (Figure 1D(Top)), such that the new intervention distribution becomes P ̃ ( X t − 1 2 , X t 1 ) = P ( X t − 1 2 ) ∫ P ( X t 1 | X t ref 2 ′ ) P ( X t ref 2 ′ ) d X t ref 2 ′ = P ( X t − 1 2 ) P ̃ ( X t 1 ) . The difference between feeding independent copies of different marginal distributions and the resulting baseline joint probability are illustrated in Figures 1B, C. We thus argue that this is a better reference scenario for testing the influence of an event between two brain regions because it accounts for the event-related variations of the input distribution relative to baseline activity.

Granger causality (GC), as well as its information-theoretic extension, Transfer Entropy (TE) is rooted in Wiener’s principle of causality. For the bivariate case, Granger (1969) defines the statement of (Granger-)causality from X ² to X ¹ whenever knowledge of X p , t 2 , in addition to X p , t 1 , yields a strictly better prediction of X t 1 . This can be interpreted as a comparison between two prediction scenarios:

• Scenario 1: predict X t 1 using both X p , t 1 and X p , t 2 ,

The predictive model describing the first scenario is referred to as the full model (Geweke, 1984) and corresponds to Eq. 2a of the SVAR model, where the first variable X ¹ is dependent on both variables X ¹ and X ². An estimate of the innovation variance of X t 1 ( σ 1 , t 2 in Eq. 2a) is the mean squared residual error ( σ ̂ 1 , t 2 ) of the forecast of X t 2 under the assumption that both X p , t 1 and X p , t 2 contribute to X t 1 . Under Scenario 2 where X t 1 is predicted only by X p , t 1 , we have a reduced model X t 1 = a t ′ ⊤ X p ′ , t 1 + η t 1 ′ , η t 1 ′ ∼ N k ′ 1 , σ ′ 1 , t 2 . (4) where the model order p′, the coefficient a′, the innovations mean k ′ 1 and innovations variance σ ′ 1 2 are different from the corresponding terms in Eq. 2a and are classically re-estimated.

If X ² Granger-causes X ¹, then the full model should fit the data more accurately compared to the reduced model as measured by the estimated variance σ ′ ̂ 1 , t 2 , which should be larger than σ ̂ 1 , t 2 . Then the amount of Granger causality can be defined as the log ratio of the residual variance between the reduced model and the full model, which leads to estimating the magnitude of Granger causality as GC X t 2 → X t 1 = 1 2 log σ ′ ̂ 1 , t 2 σ ̂ 1 , t 2 , (5) where the factor 1/2 is chosen for consistency with TE (see Section 2.2.2). While the above linear SVAR model is the most widely used, Granger causality has been extended to non-linear models following the same predictive approach (e.g., Marinazzo et al., 2008; Marinazzo et al., 2011b; Diks and Wolski, 2016; Wismüller et al., 2021, for a recent review see Shojaie and Fox, 2022).

TE is an information-theoretic implementation of Wiener’s principle, where a comparison between the prediction performance of the above two scenarios is quantified with conditional entropy. TE quantifies to which amount X ² causes X ¹ in the Granger sense and is defined by the entropy difference TE X t 2 → X t 1 = H X t 1 | X p , t 1 − H X t 1 | X p , t 1 , X p , t 2 . (6) Interestingly, using the Kullback-Leibler (KL) divergence D _KL between two probability densities D K L ( p ‖ q ) = ∫ p ( x ) log p ( x ) q ( x ) d x , TE can be rewritten as an expected KL-divergence between the corresponding conditional probabilities, thereby contrasting the two above mentioned scenarios: TE X t 2 → X t 1 = E X p , t 1 , X p , t 2 D K L p X t 1 | X p , t 1 , X p , t 2 ‖ p X t 1 | X p , t 1 . (7) As noticed by Barnett et al. (2009), under stationary Gaussian SVAR assumptions the analytic expression of Gaussian entropy applied to Eq. 6 leads to GC ( X t 2 → X t 1 ) = TE ( X t 2 → X t 1 ) in the limit of unbiased variance estimation, such that TE appears as a strict generalization of GC, and can be estimated by GC in the context of Gaussian SVAR models. TE and GC statistics are two commonly used measures of causal strength for investigating interactions between brain regions (e.g., Vicente et al., 2011; Besserve et al., 2010, 2015; Ding et al., 2006; Wibral et al., 2013, 2014; Barrett et al., 2010; Wen et al., 2013; Shorten et al., 2021; Cekic et al., 2018, with several widely applied toolboxes such as Barnett and Seth, 2014; Montalto et al., 2014; Lizier, 2014; Wollstadt et al., 2019). They generalize easily to more than two signals by including also the past of additional signals in the prediction equations when assessing causality for a specific pair, as is done in conditional pairwise Granger causality (Barrett et al., 2010; Faes et al., 2011; Runge et al., 2012; Barnett and Seth, 2014). Based on the observational conditional distribution of the neural signals being analyzed, these two measures estimate a quantity that is easily interpretable from a forecasting perspective. However, they have some limitations with regard to their interpretability as interventions in the SCM framework and in the time varying setting that interests us in the present paper.

A key issue is that the reduced model ignores but does not remove the influence of past values of X ² ( X p , t 2 ) on X t 1 by marginalizing with respect to them. It can be shown that such change does not preserve the SCM structure, and leads to violations of the Markov properties due to the implicit dependency on the mechanisms relating X p , t 2 and X p , t 1 , which manifest themselves through the p ( X p , t 2 | X p , t 1 ) term of the marginalization equation (Ay and Polani, 2008; Janzing et al., 2013): p X t 1 | X p , t 1 = ∫ p X t 1 | X p , t 1 , X p , t 2 p X p , t 2 | X p , t 1 d X p , t 2 . (8)

As a consequence, the reduced model cannot be generally interpreted as an intervention on the original SCM that would result in a model where arrows associated to the causal influence of interest would be removed. In addition, in case of bi-directional coupling, the reduced model of Eq. 4 is misspecified (in a generic case) for any finite order. This can be seen easily by exploiting the d-separation criterion (Supplementary Section SA), as illustrated in Figure 2. Figure 2B shows the estimation in the full model, where conditioning on both X p , t 1 and X p , t 2 blocks all the paths from X t − 3 1 to X t 1 such that X t − 3 1 and X t 1 are conditionally independent. For such a uni-directionally-coupled system, a finite order for the reduced model also guarantees such conditional independence, as seen in Figure 2C where all paths are blocked by conditioning. However, in the same system with bi-directional coupling, for any k > p (i.e., k > 2), there is always a path from X t − k 1 to X t 1 going through nodes of X ² that is unblocked by ( X t − p 1 , … , X t − 1 1 ) . As Figure 2D shows, 2 paths from X t − 3 1 to X t 1 are not blocked by conditioning on X p , t 1 . Under faithfulness assumptions, this implies that there is conditional dependence between X t 1 and its remote past samples, no matter how many finite past states we are conditioning on. This further implies that to minimize the forecast error of X t 1 in the reduced model one should ideally exploit the past information of this time series up to p = +∞.

This issue has been both raised and addressed in the literature, in particular by resorting to Autoregressive Moving Average models and state space models for defining an appropriate reduced model (e.g., (Barnett and Seth, 2015; Solo, 2016)). However, this remains an important limitation when extending TE to time-varying versions, where the model is assumed to be stationary at best locally in time. For example, when defining a non-stationary SVAR model as Eq. 2a, we assume a different linear model in each 1-point time window. The non-locality of TE is particularly problematic for such a time-varying model assumption because of the implicit influence of past activities on this quantity.

To overcome the limitations of TE and GC, Ay and Polani (2008) have proposed a measure of information flow to quantify the influence of some variables on others in a system, which has been further studied and generalized in Janzing et al. (2013) as a measure of the Causal Strength (CS) of an arbitrary set of arrows in a graphical model. In the present paper, we define CS in the context of time-inhomogeneous vector autoregressive processes and their associated unrolled causal graph, and thus call it Dynamic Causal Strength (DCS).

DCS can be naturally defined using the SCM interventional formalism (Pearl (2000); Peters et al. (2017), see Section 2.1.2). Briefly, interventions are performed on nodes in order to remove the specific arrows from the causal graph whose influence we wish to quantify. In agreement with Ay and Polani (2008) and Janzing et al. (2013), in the context of inhomogeneous SVAR models (as illustrated in Figure 2A), an appropriate intervention to remove the causal influence from X p , t 2 to X t 1 can be designed as the following intervention (shown in Figure 2E): remove the arrow X p , t 2 → X t 1 by injecting instead X p , t 2 ′ , an independent copy of X p , t 2 with the same joint distribution, in the original mechanism P ( X t 1 | X p , t 1 , X p , t 2 ) . The intervention distribution p ^DCS models the post-interventional world after removing the causal arrow from X p , t 2 to X t 1 and results in the entailed conditional probability p D C S X t 1 | X p , t 1 = p d o X t 1 : = f X p , t 1 , X p , t 2 ′ , η t 1 X t 1 | X p , t 1 , X p , t 2 = ∫ p X t 1 | X p , t 1 , X p , t 2 p X p , t 2 d X p , t 2 , which does not depend on p ( X p , t 2 | X p , t 1 ) anymore, in comparison to Eq. 8. DCS then quantifies the KL divergence between the distributions of X t 1 | ( X p , t 1 , X p , t 2 ) obtained in both worlds, such that DCS X t 2 → X t 1 = E X p , t 1 , X p , t 2 D K L p X t 1 | X p , t 1 , X p , t 2 ∣ p DCS X t 1 | X p , t 1 . (9) A parametric formulation under linear Gaussian model assumptions is given in Supplementary Section SD.5. Generalization to more than two time series is also straightforward following (Janzing et al., 2013): p ^DCS and DCS are simply computed by also including conditioning on the past of all other time series, in addition to X p , t 1 .

Remark: In contrast with Janzing et al. (2013), but in line with Ay and Polani (2008), we do not use jointly independent copies of each component of X p , t 2 , that is, the copy preserves the dependency between the successive past time points of X ². Indeed, Janzing et al. (2013) require having copies with jointly independent components in order to assess the individual strength of each arrow in the causal graph, which would correspond to the influence of each time lag in our setting. In contrast, this is not a requirement for us as we are only interested in assessing the overall effect of the whole past of a given time series on the another. One benefit of our choice is that it is consistent with the definition of TE: one can easily check that in absence of dependency of X t 1 on its own past, both TE (based on Eq. 6) and DCS reach the same value: the mutual information of X p , t 2 and X t 1 . Given that successive samples may be strongly correlated in practice, our choice avoids unnecessary discrepancies between these two measures to focus on their key difference. In additional, our choice can be seen as in line with Janzing et al. (2013) when considering a state representation of the time series’ causal graph, where the node of variable k at time t would be the vector [ X t k , X t − 1 k , … , X t − p + 1 k ] ⊤ .

Besides, it has also been pointed out that the definition of TE in Eq. 7 has some other non-intuitive implications (Ay and Polani, 2008; Janzing et al., 2013). In particular, there are situations in which TE(X ² → X ¹) almost vanishes, although the influence is intuitively clear. How frequent are the practical situations in which we have these detrimental effects is unclear; however, theoretical analysis suggests that this can happen when the time series are strongly correlated.

To see this, we can derive from Eq. 7 the case where X ² is a deterministic function of X ¹ such that TE vanishes. Take the special case where X t 2 is proportional to X t 1 such that X t 2 = k X t 1 , representing a time-wise synchronization of the two signals, the conditional variance will be Σ X p 2 | X p 1 = Σ X p 2 − Σ X p 2 X p 1 Σ X p 1 − 1 Σ X p 1 X p 2 = Σ X p 2 − k Σ X p 2 ⋅ 1 k 2 Σ X p 2 − 1 ⋅ k Σ X p 2 = 0

Plugging into Supplementary Equation S8 in Supplementary Section SD.4 yields, TE X t 2 → X t 1 = log b t ⊤ Cov X p , t 2 | X p , t 1 b t + σ 1 , t 2 σ 1 , t 2 = log σ 1 , t 2 σ 1 , t 2 = log ⁡ 1 = 0

However, a strong correlation between two observed time series does not necessarily imply that causal interactions between them are weak, from an SCM perspective. We will investigate this case in Section 3.1 and compare with the results of DCS to show that DCS does not suffer from this non-intuitive vanishing problem.

While several intuitive properties make DCS a good candidate to quantify causal influences, we exhibit a counterintuitive property common to TE and DCS in the context of peri-event time series. Transient neural events are mainly investigated in two types of analyses: 1) stimulus-triggered (or response-triggered) data that are temporally aligned by task (or response) onset and 2) event-triggered data where occurrences of a type of brain-activity pattern are detected along the time course of the recordings (manually or algorithmically) and used to create peri-event trials.

In both cases, neural activities are likely to have a deterministic component appearing in the peri-event ensembles, due the similarity of the response to successive stimuli in case 1), or due to the similarity of the neural patterns detected in the recordings in case 2). Here we will show that, in a linear setting, TE and DCS are insensitive to such a deterministic component. Specifically, TE and DCS values are unaffected by interventions on the innovations’ mean at any time point.

First, we exhibit the role played by a deterministic perturbation in an example.

Example 1. Consider the bi-variate SVAR(1) model in the following form X t 1 : = a X t − 1 1 + b X t − 1 2 + η t 1 , (10a) X t 2 : = η t 2 , (10b) with a, b ≠ 0 and a stationary innovation for X ¹ , η t 1 ∼ N ( 0,1 ) , but a non-stationary innovation for X ² , η t 2 ∼ N ( α δ t , t 0 , σ 2 , t 2 ) , with δ t , t 0 = 1 , for   t = t 0 , 0 , otherwise .

When varying α, this models a intervention on the second time series. Then it can be easily shown that the expected time course of X ¹ is E X t 1 = α b a t − t 0 + 1 , t ≥ t 0 + 1 0 , otherwise. This witnesses the causal influence of X t 0 2 on values of X t 1 at subsequent times, which for large α results in large deviations from the baseline expectation of X t 1 for t prior to t ₀. Intuitively, one may expect that a quantification of the magnitude (strength) of the causal influence of X ² on X ¹ should be larger for larger α, as a transient of larger magnitude propagates from X ² to X ¹. From a neuroscientific perspective, this could model an experimental setting where one brain region is electrically stimulated with increasing strength to detect whether it is anatomically connected to another. Obviously, the magnitude of the stimulation is expected to be critical to elicit a response in the target region. However, TE and DCS actually turn out to be insensitive to such stimulation.

We will show this in the more general setting of the SVAR(p) model of Eq. 2a and Eq. 2b.

Proposition 1. For linear SVAR models defined by Eq. 2a and Eq. 2b, TE and DCS measures are invariant to deterministic perturbations, i.e., to changes in the mean of the innovation’s distributions ( k t 1 , k t 2 ) .

Proof. Without loss of generality, we will show invariance to an elementary intervention at a single time t ₀ that transforms η t 0 2 to η t 0 2 + α , which boils down to changing the mean parameters of the innovation k t 0 2 in Eq. 2a and Eq. 2b. By linearity and symmetry of the problem for channel 1 and 2, invariance to deterministic perturbations results from combining several elementary interventions.

To compute how the intervention distribution of the new variables denoted ( X ̃ 1 , X ̃ 2 ) changes with respect to the distribution of the original variables, we can examine the difference with respect to (X ¹, X ²) that has the same innovations, except for η t 0 2 for which we remove a constant α. (X ¹, X ²) is then distributed according to the original distribution (before intervention), and the difference ( U , V ) = ( X ̃ 1 − X 1 , X ̃ 2 − X 2 ) follows the equations U t = a ⊤ U p , t + b ⊤ V p , t , V t = c ⊤ U p , t + d ⊤ V p , t + δ t , t 0 , which is a set of linear deterministic difference equations with a unique solution making X and X ̃ coincide before the intervention ¹ (U _t , V _t ). As a consequence, by linearity, the interventional density p ̃ is a shifted version of the original: p ̃ X t 1 , X p , t 1 , X p , t 2 = p X t 1 − U t , X p , t 1 − U p , t , X p , t 2 − V p , t , which implies the same for conditional marginal densities, e.g., p ̃ X t 1 | X p , t 1 , X p , t 2 = p X t 1 − U t | X p , t 1 − U p , t , X p , t 2 − V p , t

and p ̃ X t 1 | X p , t 1 = p X t 1 − U t | X p , t 1 − U p , t .

As a consequence TE on the intervention distribution writes TE X ̃ t 2 → X ̃ t 1 = ∫ p ̃ X t 1 , X p , t 1 , X p , t 2 log p ̃ X t 1 | X p , t 1 , X p , t 2 p ̃ X t 1 | X p , t 1 d X t 1 d X p , t 1 d X p , t 2 = ∫ p X t 1 − U t , X p , t 1 − U p , t , X p , t 2 − V p , t × log p X t 1 − U t | X p , t 1 − U p , t , X p , t 2 − V p , t p X t 1 − U t | X p , t 1 − U p , t d X t 1 d X p , t 1 d X p , t 2 . And by change of variable we get the invariance property: TE X ̃ t 2 → X ̃ t 1 f = ∫ p X t 1 , X p , t 1 , X p , t 2 log p X t 1 | X p , t 1 , X p , t 2 p X t 1 | X p , t 1 d X t 1 d X p , t 1 d X p , t 2 = TE X t 2 → X t 1 , which can be generalized to arbitrary deterministic perturbations. The same reasoning can be applied to DCS leading to invariance as well (Supplementary Section SB) and this concludes the proof.

Arguably, this result is not what we would expect from an event-related measure of influence, because in the above example of Eq. 10a and Eq. 10b, setting a large α intuitively leads to a large influence of X ² on X ¹ provided b ≠ 0. Provided that TE and DCS can be made arbitrarily small by reducing the innovation’s variance σ 2 , t 2 (according to their analytical expression in Supplementary Section SD), TE and DCS may detect no influence despite this strong effect on the mean of X t 2 . Although this invariance result is rigorously derived for linear SVAR models, it uncovers an issue for non-linear models as well, the magnitude of the causal influence associated to deterministic perturbation then depending chiefly on the non-linear properties of the system under study, and not on the magnitude of the changes triggered by the perturbation. Moreover, linear SVAR(p) models being able to approximate nonlinear dynamics, this suggests that deterministic causal influences cannot be detected by TE or DCS for a broad class of models in practice.

As elaborated above, this is in contrast to what would be expected in the neuroscientific context, and directly relates to the observational, event-related setting that we investigate: the deterministic component is due to the alignment of the data with respect to an event of interest, and we do not have a different condition to contrast the occurrence of this event with what would have happened in its absence. This analysis calls for building a synthetic baseline condition that would allow deterministic changes to be detected.

Following the guidelines for event-based causality (presented in Section 2.1), we propose a novel measure, the relative Dynamic Causal Strength (rDCS), as a modification of DCS. This measure aims at taking into account the influence of event-based changes in the cause signals independent from the connectivity (the mechanism), and notably those driven by deterministic exogenous inputs. In the specific problem we are investigating, the cause is the past states of X ², denoted X p , t 2 , while the mechanism can be represented by the model in Eq. 2a and symbolized by the corresponding causal arrow in the causal graph. DCS only deletes the causal arrow in the post-intervention scenario but preserves the event-related change in the cause itself.

In the case where X ² is driven by a deterministic exogenous input in a transient window, the cause exhibits significant changes relative to baseline; thus, intuitively, the causal effect should also be enhanced even if the causal arrow remains the same (i.e., the coefficient b stays unchanged). Apart from intervening on the causal arrow, further intervention can be implemented on the cause node to construct a post-intervention scenario where the cause receives no time-varying innovations. Therefore, inspired by causal impact (Section 2.1.3) which characterizes the difference between the current state and a baseline state, we propose (additionally to DCS) to replace the marginal of X p , t 2 by the marginal of X p , t ref 2 , that we denote p ref ( X p , t ref 2 ) , for a reference time t _ref . The reference time t _ref is typically chosen to be a stationary period before the occurrence of the transient deterministic perturbations and statistics of X p , t ref 2 can be averaged by statistics of X p , t 2 within this period. This leads to the relative Dynamic Causal Strength (rDCS) rDCS X t 2 → X t 1 = E X p , t 1 , X p , t 2 D K L p X t 1 | X p , t 2 , X p , t 2 p d o X t 1 : = f X p , t 1 , X p , t ref 2 , η t 1 X t 1 | X p , t 1 , X p , t 2 (11)

with p d o X t 1 : = f X p , t 1 , X p , t ref 2 , η t 1 X t 1 | X p , t 1 , X p , t 2 = ∫ p X t 1 | X p , t 1 , X p , t 2 p ref X p , t 2 d X p , t 2 (12)

The implementation of rDCS given a SVAR model is derived in Supplementary Section SD.6. Generalization to more than two time series can be done in the same way as for DCS, by including extra conditioning on the past of other time series for all quantities.

Intuitively, the term relative originates from the comparison between the current past state X p , t 2 and the reference past state X p , t ref 2 . It is then natural to predict that in the uni-directional case, rDCS(X ² → X ¹) = DCS(X ² → X ¹) for any reference time t _ref if X ² is stationary because stationarity implies that the marginal distributions of X p , t ref 2 and X p , t 2 are identical. As a particular case, this result implies that a transient loss of causal link from X ² to X ¹ will lead to rDCS = 0, while for a stationary bivariate system, DCS = rDCS is constant.

The definition of rDCS implies sensitivity to deterministic perturbations. Indeed, taking the example in Section 2.3.2, the reference state X p , t ref 2 is unaffected by the deterministic perturbation. Consequently, the translational invariance does not hold for the intervention distribution because rDCS X ̃ t 2 → X ̃ t 1 = ∫ p ̃ X t 1 , X p , t 1 , X p , t 2 log p ̃ X t 1 | X p , t 1 , X p , t 2 ∫ p ̃ X t 1 | X p , t 1 , X p , t 2 p ̃ ref X p , t 2 d X p , t 2 d X t 1 d X p , t 1 d X p , t 2 = ∫ p X t 1 − U t , X p , t 1 − U p , t , X p , t 2 − V p , t log p X t 1 − U t | X p , t 1 − U p , t , X p , t 2 − V p , t ∫ p X t 1 − U t | X p , t 1 − U p , t , X p , t 2 − V p , t p ref X p , t 2 d X p , t 2 d X t 1 d X p , t 1 d X p , t 2 = ∫ p X t 1 , X p , t 1 , X p , t 2 log p X t 1 | X p , t 1 , X p , t 2 ∫ p X t 1 | X p , t 1 , X p , t 2 p ref X p , t 2 + V p , t d X p , t 2 d X t 1 d X p , t 1 d X p , t 2 ≠ ∫ p X t 1 , X p , t 1 , X p , t 2 log p X t 1 | X p , t 1 , X p , t 2 ∫ p X t 1 | X p , t 1 , X p , t 2 p ref X p , t 2 d X p , t 2 d X t 1 d X p , t 1 d X p , t 2 = rDCS X t 2 → X t 1 , because p ̃ ref is not translated by the deterministic perturbation in the way p ̃ ( X p , t 2 ) is (as the perturbation is happening after the reference time), such that the denominators do not allow equating the integrated terms by change of variables in the generic case. Therefore rDCS is capable of uncovering transient causal influences between stimulus-triggered events exhibiting a deterministic waveform.

We simulated a non-stationary uni-directionally-coupled autoregressive system defined in Eq. 2a and Supplementary Eq. 2b. The causal direction is X ² → X ¹. The system is designed as a bivariate SVAR(4) process with a time-invariant coefficient matrix: a ^⊤ = [−0.55, − 0.45, − 0.55, − 0.85], b ^⊤ = [1.4, − 0.3, 1.5, 1.7], c ^⊤ = [0, 0, 0, 0] and d ^⊤ = [0.9, − 0.25, 0, 0.25]. These coefficients were randomly generated and kept after checking the stability of the SVAR(4) system. Uni-directional interactions are ensured by setting the autoregressive coefficients associated to interactions in the opposite direction (i.e., c) to zero for all lags.

We enforce non-stationarity of η t 2 , the innovations of the ground truth cause process { X t 2 } . Both innovations η t 1 and η t 2 are drawn from a Gaussian distribution with unit variance (with no correlation in between, i.e., C o v [ η t 1 , η t 2 ] = 0 ); the difference is that E [ η t 1 ] = k t 1 = 0 while E [ η t 2 ] = k t 2 is non-zero and time-varying. We designed the time-varying profile of k t 2 as a Morlet-shaped waveform to mimic the oscillatory properties of neural event signals: k t 2 = H ⁡ exp ( − ( α x ) 2 / 2 ) cos ( 5 α x ) , where α = 2/25 is a constant controlling the event duration, and H = 4 is the amplitude of the highest peak in the center of the event. The total duration of the Morlet-shaped waveform is 101 ms. The innovation’s mean designed for X ² is shown in Figure 5B (top left panel).

We generated this bi-variate SVAR(4) process for 1300s consisting of 5,000 trials of perturbation events by transiently varying η t 2 , detecting event occurrence based on the cause X t 2 , as illustrated in Figure 5A. The central peaks of these Morlet events are used as the ground-truth reference points for which peri-event time t′ = 0, and used to extract a dataset of multi-trial events ensemble with a 200-ms peri-event window such that t′ ranges, from −99 ms to +100 ms (i.e., there is no alignment procedure that could lead to selection bias, see Section 2.5). The event waveforms of the cause variable X t 2 and the effect variable X t 1 are illustrated in Figure 5B (bottom left, middle left). The whole process is repeated 100 times to obtain variabilities plotted in the figure.

The designed deterministic innovation (i.e., identical across trials), can be seen as imposing a hidden state evolving across the peri-event interval. The event ensembles obtained by this ground truth model define a dataset where no event selection and alignment is needed. We can compare the SVAR model estimation and causality measures resulting from this dataset to the outcomes obtained by selecting and aligning events based on either variable X ² or X ¹, as discussed in Section 2.5.

To validate the recoverability theory in the presence of selection bias due to the event detection procedure, we test the single-time selection setting (see Section 2.5) where sub-threshold trials are removed from the ground truth peri-event dataset (as illustrated in Figure 3D(Top right)), thus preserving the ground-truth hidden states (Figure 5B(Top middle)). The peri-event trials having reference point values higher than a threshold d ₀ = 3SD for the chosen variable are selected, where the standard deviation is computed from the whole signal. The selected event ensembles are shown in Figure 5A(middle center) for thresholding based on X t 1 and Figure 5A(middle right) for thresholding based on X t 2 . Notably, this kind of selection is only feasible when the hidden state in known, which is not realistic practically for real data.

Next we demonstrate the appropriateness of the approach performed on real data (i.e., selection and smoothed alignment based on putative cause), we set d ₀ as a threshold and performed smoothed alignment over the original signal itself. We obtain an event ensemble by selecting local peaks for points over d ₀ as new reference points, which is shown in Figure 5B(middle right and bottom right). This can be seen as a smoothed version of the ground-truth dynamics, which is also confirmed by checking the aligned hidden states (Figure 5B(Top right)).

While inferring SVAR model parameters of the event ensembles according to Shao et al. (2022), the true model order 4) can be recovered for all five ensembles. Figure 5C demonstrates the recoverability of conditional probabilities for ensembles aligned by the putative cause. Coupling strengths from the putative cause to the putative effect are plotted in red. As described in the simulation procedure in Section 3.2.1, the coupling strength is constant over time, which is reflected in Figure 5C(left). Consistent with the theory in Section 2.5, biased selection of event trials on the samples at t′ = 0 leads to unbiased estimation of the coupling strength X t 2 → X t 1 aligned by the cause X ² (denoted also as “∣X ²” in Figure 5C(middle)). By comparison, the coupling strength in the other direction is slightly biased at negative peri-event times (t′) but still relatively close to its true value 0). This contrast holds for alignment with local peaks over threshold, as seen in Figure 5C(right).

Figure 5D (Top, middle, bottom) shows the corresponding results of how causality measures perform in the three alignment scenarios. For clearer visualization of TE and DCS with a zoomed vertical scale, see Supplementary Figure S8. During the periods where no transient events occur, all three measures are able to infer a time-invariant causal effect in the ground-truth direction (X ² → X ¹) compared to the opposite direction. Besides, in line with theoretical predictions, DCS is higher than TE and is equal to rDCS. During the perturbation events, in the ground truth direction TE and DCS remain constant and rDCS exhibit a rhythmic pattern. These results match the theoretical predictions: TE and DCS measures the connectivity strength, which does not change, while rDCS measures the combined causal effect related to the connectivity and the event-based changes at the cause, yielding larger variations transmitted to the effect node.

In the transient time scale, thresholding leads to selection bias in estimating causality measures. In the case where event ensembles are aligned by single-time selection of the cause X t 2 , TE and DCS of the ground truth direction is underestimated while rDCS is slightly overestimated around t′ = 0. A bias appears in the opposite direction while aligned by the effect, but the direction of causation is detected correctly. The case of local peak alignment shows similar results, except the peak amplitude of the smoothed rDCS is less amplified. Notably, in the smoothed case, a transient increase is observed in both TE and DCS, resembling the envelope of the perturbation. This is likely an effect of the smoothing procedure but is quickly interrupted by a negative bias due to thresholding, making the results unreliable in detecting transient changes. We also showed a negative example with putative effect alignment in Supplementary Figure S5, where the coupling strength in the causal direction is much weaker than the model in Figure 5. The coupling strength of the causal direction undergoes a sharp decrease at peri-event time t′ = 0, leading to a transient underestimation of TE, DCS and rDCS for both single-time and local peak scenarios. The close-to-zero value of rDCS is misleading for the inference of transient causal interactions, thus illustrating the unreliability of putative effect alignment.