Characterizing unidirectional couplings between point processes and flows

Experimental data comprising both time-continuous flows and point processes are recorded in many scientific disciplines. The characterization of causal interactions from such signals is key to an advanced understanding of the underlying dynamics. We therefore introduce a unified approach to characterize unidirectional couplings between point processes, between flows, as well as between point processes and flows. For this purpose we show and exploit the generality of the asymmetric state similarity conditioning principle. We use Hindmarsh-Rose neuron models and Lorenz oscillators to illustrate the high sensitivity and specificity of our approach.

Introduction. -The reliable detection of directional couplings from experimental signals is crucial for the study of many systems in nature. These signals can for example be given by sequences of discrete event times of point processes or by time-continuous variables measured from flows. Often experiments yield simultaneous recordings comprising both point processes and flows. A prominent example from neuroscience is the recording of spiking activity of individual neurons and local field potentials. The characterization of directional interactions from these different representations of neuronal dynamics can progress the understanding of neuronal information processing. While various approaches exist to detect directional couplings between point processes and between flows (see [1,2] and references therein), much less attention has been paid to the hybrid case of directional couplings between point processes and flows [3]. In this letter we therefore introduce a unified approach that allows us to characterize unidirectional couplings between pairs of point processes, pairs of flows, as well as between point processes and flows. For this purpose we generalize the concept of nonlinear interdependence measures, which was introduced to detect unidirectional couplings from delay reconstructions obtained from pairs of flows.
Nonlinear interdependence. -Nonlinear interdependence measures are based on the principle of asymmetric state similarity conditioning (see [2,4]). Suppose (a) E-mail: ralphandrzejak@yahoo.de that we have two stationary deterministic dynamics X and Y with a unidirectional coupling that is not strong enough to induce a synchronized motion of the two dynamics. Importantly we assume that the sampling is synchronous, i.e. we have states x(t i ) and y(t i ) at times t i for i = 1, . . . , N. Furthermore, suppose that we have a symmetric dissimilarity measure between states attained by individual dynamics at different times: d X ij = d X (x(t i ), x(t j )) = d X (x(t j ), x(t i )) for i = 1, . . . , N and j = 1, . . . , N. Let D X ij stand for the dissimilarity matrix containing all N × N pairwise dissimilarities d X ij . The quantities d Y ij and D Y ij are defined analogously, but importantly the definition of the dissimilarity measure d Y ij does not have to be the same as d X ij . To test for asymmetric state similarity conditioning we quantify, across i and j, the degree to which a low value of a certain d Y ij implies a low value of the corresponding d X ij . For couplings X → Y this is expected to be higher than for uncoupled dynamics and higher than the degree to which small elements in D X ij are mapped to small elements in D Y ij . Evidently, to explore this criterion based on simultaneous recurrences, the aforementioned synchronous sampling of X and Y is indispensable.
Various nonlinear interdependence measures quantify this signature of directional couplings. Here we use the measure L, which was shown to be of higher sensitivity and specificity than a number of previous approaches [2]. To calculate L we carry out the following steps. For each fixed i and j we denote by g ij the rank which d X ij takes in a sorted ascending list of all dissimilarities included in 50012-p1 the set {d X il } l=1,...N,|i−l|>W . The number of elements in this set is denoted by M i (see footnote 1 ). The parameter W is used to exclude pairs of states that are similar to each other simply because they are close to each other in time (see [5]). We now turn to the Y dynamics and look up the time index j 0 of the spatially nearest neighbor of the state with time index i by determining d Y ij0 = min{d Y il } l=1,...N,|l−i|>W . Again temporally close states within W are excluded. Merging information from both dynamics we define the conditioned rank G i (X|Y ) = g ij0 . Above we indicated that for a coupling X → Y , a low value of d Y ij implies a low value of the corresponding d X ij . Accordingly, for couplings X → Y the quantity G i (X|Y ) will tend to be small, bounded by a minimum of 1. For independent dynamics, in contrast, j 0 can be regarded as a random sample from 1, . . . , M i . In consequence, the expected value of G i (X|Y ) for independent dynamics is G i = Mi+1 2 . Based on these considerations we define 2 : The quantity L(Y |X) is calculated by switching X and Y in the above definitions, and we use ∆L(X, Y ) = L(X|Y ) − L(Y |X) [2]. Values of L(X|Y ) and L(Y |X) are expected to be distributed around zero for independent realizations of independent dynamics. For couplings X → Synchronous distance matrices for point processes and flows. -In the following we introduce the dissimilarity matrices D X ij and D Y ij used as input for the nonlinear interdependence measure L. At this point we do not need to specify whether X and Y are flows or point processes. We only require that both dynamics are measured simultaneously for a duration of Q. Key is that the required synchronous sampling of both dynamics is reached by an isochronous segmentation of the time interval [0, Q] that is common to both dynamics. We therefore at first define a sequence of overlapping time intervals of length q and step size s with s q Q. The i-th interval spans [(i − 1)s, (i − 1)s + q] for i = 1, . . . , N, where N is obtained by rounding down Q−q s + 1. For each segment i we measure time relative to its beginning τ = t − (i − 1)s.
First we deal with a pair of point processes X P and Y P and assume that both processes are aperiodic and that 1 M i = N − 2W − 1 for W < i < N − W + 1. Below and above this range M i increases in steps of 1 and reaches M i = N − W − 1 at i = 1 and i = N . 2 Note that in the general definition of the measure L(X|Y ) and related approaches the index j 0 is determined not only for one nearest neighbor but for k nearest neighbors. See [2] and references therein.
individual event times in general do not coincide across the two processes. As an analog to delay coordinates, vectors of subsequent inter-event intervals can be used to reconstruct the dynamics underlying the individual point processes [6]. While these reconstructions can serve as basis for univariate nonlinear signal analysis of the individual dynamics [6,7], inter-event interval vectors cannot be used as input for nonlinear interdependence measures since they are not synchronous across X P and Y P . Therefore, instead of reconstructing the dynamics underlying a point process, we directly quantify the dissimilarity of segments of the event time sequence.
Various approaches can be applied to quantify the dissimilarity of the event time sequences, and we use the recently introduced so-called ISI spike train distance [8]. This particular approach evaluates the ratio of the instantaneous inter-event intervals and has the advantage of being parameter free and time scale independent. Let the point process X P comprise subsequent event times within the interval 0 t Q. For the span of relative time 0 τ q of segment i, let I x i (τ ) denote the length of the interval between the last previous and first subsequent event. These previous and subsequent events can be located in segment i or, in adaptation of the original approach, in previous or subsequent segments. Prior to the overall first event across all segments the start of the inter-event interval is defined as 0. Posterior to the overall last event across all segments the end of the inter-event interval is defined as Q. Let If at relative time τ the two segments with index i and j are similar to each other, in the sense that their instantaneous inter-event interval I x i (τ ) and I x j (τ ) are of similar length, the quantity R x ij (τ ) will tend to small values. Averaging across the relative time yields the dissimilarity between the two segments with indexes i and j: d values. The same steps of analysis are carried out for the point process Y P resulting in D Y P ij . Such a definition of an isochronous dissimilarity matrix for point processes has some analogy to [9] where recurrence plots were used to study univariate marked point processes. For marked point processes conventional spike train distances cannot be applied, and [9] proposed a corresponding adaptation.
Let us now consider a pair of flows X F and Y F . For the flow X F let x (i) (τ n ) denote the signal amplitudes recorded during segment i. The temporal resolution is determined by the sampling time of the flow as τ n = n∆t, with integer n. The dissimilarity between segment i and j is defined by the temporal average Again the same steps of analysis are carried out for the flow ij are based on the same segmentation of the signals, they can be used as input for the nonlinear interdependence measure L. We have to distinguish four different pairings. Apart from the cases when both dynamics are point processes (X P , Y P ) or flows (X F , Y F ), we have to distinguish the two mixed cases (X P , Y F ) and (X F , Y P ). This is necessary since in general results depend on whether the driver or the response dynamics is a point process or a flow. We identify specific cases through the argument of the measures, e.g. L(X P |Y F ). When we refer to findings that hold independent of the type of dynamics we use the general notation L(X|Y ), etc. In all examples we call the driver X and the response Y . Accordingly, we expect ∆L(X, Y ) > 0 for nonzero couplings.
Application to Hindmarsh-Rose neurons. -As first example we use Hindmarsh-Rose neuron models (see [10] and references therein). Here the dynamics X Ḣ drives the dynamics Y Ḣ y 3 (t) = 0.0021(−y 3 (t) + 4(y 1 (t) + 1.6)), with coupling strength ε via the coupling functioṅ with We integrated these equations applying a fourth-order Runge-Kutta algorithm with step size of 0.1 time units and sampling step t n = n∆t with n integer and ∆t = 0.2 time units. As unit interval we define T = 1000∆t. Apart from the uncoupled case (ε = 0) we used 60 coupling values within 10 −4 < ε < 0.24 equidistantly distributed on a logarithmic scale. For each ε we generated 100 independent realizations of the dynamics starting from random initial conditions. The time interval covering the first 500T was discarded to let transients of the dynamics die out. We use · to denote mean values across the 100 independent realizations. Being defined by a set of differential equations Hindmarsh-Rose neurons represent The thresholds for the definition of the events (see [10]) are depicted as horizontal lines. a deterministic flow. Therefore, we can use the first components of the driving X H dynamics and driven Y H dynamics as flow variables and denote them by X F H : x 1 (t n ) and Y F H : y 1 (t n ), respectively. On the other hand, as a model for the integrate-and-fire behavior of neurons this dynamics is a prototype of a point process where the events are given by the spikes. Accordingly, we define the event times of a point process to be denoted by X P H from the times of upward threshold crossings of x 1 (t n ). Analogously, Y P H is defined from threshold crossings of y 1 (t n ). In both cases a value of Θ = 0.6 is used for the threshold [10].
In the following we show detailed results for J x = 3.30 and J y = 3.28 (figs. 1, 2). For these parameters the driver 50012-p3 X H exhibits aperiodic spiking with 4.8 spikes per T . For the response Y H the spikes are fired in aperiodic bursts. With increasing coupling the characteristics of this bursting changes. In particular, the mean number of spikes per T increases from 6.3 at ε = 0 to 9.0 at ε = 0.24. For the calculation of the dissimilarity matrices we fix the segment length to q = T . The step size is set to s = 0.2T , resulting in an overlap of 80% of neighboring segments. To discard elements of the dissimilarity matrices that correspond to overlapping segments we set W = q s − 1 = 4. We start with results for the pair of point processes with a duration of Q = 400T . While L(X P H |Y P H ) increases with the coupling strength, L(Y P H |X P H ) remains close to zero ( fig. 2(A)). In consequence, ∆L(X P H , Y P H ) becomes significantly positive 3 for ε 0.002. Before interpreting these results as a correct detection of the coupling we have to rule out that they are caused by the aforementioned asymmetries between the characteristics of the firing patterns of X P H and Y P H or by the increasing mean firing rate of Y P H upon increasing ε. For this purpose we use replica surrogates. For each realization of X P H and Y P H we generated a replica X P * H and Y P * H by integrating the dynamics again but starting from different random initial conditions. Evidently, the characteristics of the firing patterns of X P * H and Y P * H and their dependence on the coupling for Y P * H are the same as for the original dynamics X P H and Y P H . On the other hand, Y P * H is independent from X P H . Hence, we can test the specificity of our approach by determining ∆L(X P H , Y P * H ). We find that ∆L(X P H , Y P * H ) is never significantly different from zero, i.e., no false positive detections of directional couplings are found for the replica surrogates ( fig. 2(E)).
So far we considered results for the case that both dynamics are point processes (X P H , Y P H ). Qualitatively similar results are obtained for the mixed cases (X P H , Y F H ) and (X F H , Y P H ), and also when both dynamics are flows (X F H , Y F H ) ( fig. 2(B)-(D)). The weakest increase is found in the latter case. Nonetheless, significantly positive ∆L(X F H , Y F H ) values are obtained for most ε 0.008. To further compare our results in dependence on the duration Q, we define ψ as the fraction of nonzero ε values for which ∆L(X|Y ) is significantly positive. Evidently, this ad hoc definition should only be used as a relative measure of sensitivity since the ψ values depend on the range and sampling of ε values as well as on the number of realizations and the α used for the Wilcoxon test. Across Q the highest sensitivity is obtained when both dynamics are point processes ( fig. 3). Overall lower ψ values are found for the two mixed cases and the lowest sensitivity is obtained when both dynamics are flows. Accordingly, for the Hindmarsh-Rose example the 3 To probe whether the mean ∆L(X, Y ) is significantly different from zero, we apply a Wilcoxon signed rank test. Our analysis involves several levels of multiple testing. Apart from the multiple ε values and the different types of dynamics we use different Q below. Therefore, we apply an arbitrary but very small significance level of α = 5 × 10 −6 for the Wilcoxon test.  coupling can be detected best when focussing on the event times, and the actual waveform of the Hindmarsh-Rose variables is less relevant. To reaffirm the specificity of our approach we repeated the analysis underlying fig. 3(F) for replica surrogates. Across all Q values and types of dynamics we found not a single nonzero ψ value. This very good specificity is inherited from the applied nonlinear interdependence measure L, which in contrast to a number of previous approaches has an expected value of zero for independent dynamics (cf. [2,11]). By varying J x and J y we tested other settings of unidirectionally coupled Hindmarsh-Rose neurons. These included settings where both the driver X and response Y were in the aperiodic bursting regime throughout the coupling range. In further settings, X was either in the aperiodic spiking or in the aperiodic bursting regime while for zero and low coupling Y was in the aperiodic spiking regime. Here the coupling caused that the response Y switched to aperiodic bursting. In all investigated cases we correctly detected the direction of the coupling from ∆L(X, Y ) > 0, regardless of the type of dynamics being a point process or flow. Results of the replica surrogates again confirmed the specificity of our approach.
Application to Lorenz dynamics. -To illustrate the generality of our approach, we present results for unidirectionally coupled Lorenz dynamics. Here the dynamics X L , drives the dynamics Y L , y 1 (t) = 10(y 2 (t) − y 1 (t)) + ε(x 1 (t) − y 1 (t)), For each coupling strength ε we integrate 100 independent realizations with a fourth-order Runge-Kutta algorithm. We used a step size of 0.005 time units and a sampling step of t n = n∆t with n integer and ∆t = 0.01 time units. The unit interval is again defined as T = 1000∆t. In contrast to the spiking of the Hindmarsh-Rose model the Lorenz dynamics exhibits no evident point process behavior. Nonetheless events can be defined, and the resulting inter-event interval vectors yield a reconstruction of the individual Lorenz dynamics (e.g. [6,7]). We apply two different event time definitions.
For the first definition we follow Sauer [6] and drive an integrate-and-fire process by the x 1 (t n ) variable of X L . Starting at S(t 1 ) = 0 we integrate S(t n ) = S(t n−1 ) + (x 1 (t n ) + 25)∆t. Whenever S(t n ) crosses the threshold of Θ = 12, this time t n0 is used as event time, S(t n0 ) is reset to zero, and the integration is restarted. The sequence of all threshold crossings defines the point process X P S . The underlying Lorenz variable is used as flow X F S : x 1 (t n ). The same steps of analysis are carried out for Y L resulting in Y P S and Y F S . We set Q = 350T , a duration that comprises on average approximately 2600 events for both X P S and Y P S . Here the coupling has only little influence on the number of events in Y P S . For the second definition we follow Ding and Wang [8] and trigger events for X P D by upward threshold crossings (Θ = 27) of the variable x 3 (t n ) of X L . The underlying Lorenz variable x 3 (t n ) is used as flow X F D . Again the same steps of analysis are carried out for Y L yielding Y P D and Y F D . Here we set q = 2T and Q = 350T , resulting in a total of approximately 2700 spikes for X P S . Depending on ε, between 2400 and 2700 spikes are obtained for Y P S during Q. For both definitions we set the segment length to q = 2T , the step size to s = 0.4T (80% overlap), and accordingly set W = 4.
Starting at coupling strengths ε around 1, high ∆L(X, Y ) values are obtained from x 1 (t n ) and events generated by the integrate-and-fire process driven by this component (fig. 4). This increase is smaller for x 3 (t n ) and events generated by threshold crossings of this component ( fig. 5). Nonetheless, the minimal ε leading to significantly positive ∆L(X, Y ) does not depend attain high values, and ∆L(X, Y ) decreases or even attains negative values. This is in agreement with the fact that for synchronized dynamics the coupling direction cannot be determined reliably (see [2] and references therein). Again no false positive coupling detection was found for replica surrogates.
Summary and discussion. -We propose a unified approach that allows us to detect with high sensitivity and specificity unidirectional couplings between pairs of point processes, pairs of flows, as well as between point processes and flows 4 . We used replica surrogates to probe the specificity of our approach. These replica surrogates test the null hypothesis that X and Y are independent, while no assumption is made about the properties of the individual dynamics. For experimental data it is in general not possible to generate replicas of the dynamics. However, in case multiple realizations of the experimental signals measured under the same experimental conditions are available, surrogates can be obtained by shuffling across these realizations. If that is not the case, surrogates can be obtained by applying constrained randomizations to the inter-event sequences and the time continuous variables (see, e.g., [12,13] and references therein). Alternatively, surrogates can be constructed by introducing time shifts between the signals from X and Y . Whatever surrogate approach is used, it is always essential to specify exactly what null hypothesis is tested by the surrogates. This null hypothesis is determined by the constraints that are applied for generating the surrogates.
None of the parameters of the dynamics, the event extraction methods, or the measure L were tuned to optimize our results. Key to the unification is a representation of the dynamics that is common for point processes and flows. In this letter this representation is given by the dissimilarity matrices D X ij and D Y ij which we evaluate using the nonlinear interdependence measure L. Analogously instantaneous phases can be extracted from both point processes and flows. Such instantaneous phases can be analyzed using, e.g., the directional evolution map approach [14]. This hybrid approach was applied in [3] to study cardiorespiratory interactions. A further representation that can be derived for both point processes and flows are Fourier spectra. In consequence, spectrum based Granger causality ( [1] and references therein) could be used to detect directional interactions between point processes and flows. To the best of our knowledge, however, this latter approach has not been tested yet. A priori, our approach has the advantage over these alternatives that it can be applied to irregular signals measured from nonlinear deterministic dynamics. Often a meaningful phase cannot be defined for such irregular signals, and spectrum based approaches might not be sensitive to couplings of nonlinear deterministic dynamics. * * * RGA acknowledges grant FIS-2010-18204 of the Spanish Ministry of Education and Science. TK acknowledges support by the Italian Ministry of Foreign Affairs regarding the activity of the Joint Italian-Israeli Laboratory on Neuroscience. We thank D. Chicharro and C. Rummel for useful discussions.