Oscillating systems with cointegrated phase processes.

We present cointegration analysis as a method to infer the network structure of a linearly phase coupled oscillating system. By defining a class of oscillating systems with interacting phases, we derive a data generating process where we can specify the coupling structure of a network that resembles biological processes. In particular we study a network of Winfree oscillators, for which we present a statistical analysis of various simulated networks, where we conclude on the coupling structure: the direction of feedback in the phase processes and proportional coupling strength between individual components of the system. We show that we can correctly classify the network structure for such a system by cointegration analysis, for various types of coupling, including uni-/bi-directional and all-to-all coupling. Finally, we analyze a set of EEG recordings and discuss the current applicability of cointegration analysis in the field of neuroscience.

Introduction

Since the first scientific discovery of two pendulums synchronizing by Christiaan Huygens in the seventeenth century, this naturally occurring phenomenon has now been observed in diverse areas such as fireflies synchronizing their flashing behavior, a theatre audience applauding after a show and also in chemical and biological systems, such as the brain and the heart beats of a mother and her fetus, where coupled oscillators appear, see also Pikovsky et al. (2001). Due to it’s pervasive presence, understanding synchronization is of key interest for researchers to understand biological networks, such as the connectivity of the nervous system, circadian rhythms or the cardiovascular system. To a statistician this presents a fascinating challenge of modelling complex behavior in large scale systems and how to infer the data-generating mechanisms. To this day, synchronization is not fully understood, but has been the centre of research for decades as evident in Ermentrout (1985), Kuramoto (1984), Strogatz (1987, 2000), Taylor and Holmes (1998), Winfree (1967), even the phenomenon of synchronizing pendulums as observed by Huygens, still attracts attention today, see Martens et al. (2013), Oliveira and Melo (2015). Many innovative ideas have been presented since Winfree (1967) began a mathematical treatment of the subject. When Kuramoto (1984) first presented his model of coupled oscillators, this made a huge impact in the field and spawned a new generation of research on synchronization. Kuramotos model is still considered among one of the most significant advancements in the study of synchronization in oscillating systems as acknowledged by Strogatz (2000), and the study of coupled oscillators still attracts a fair interest from researchers Ashwin et al. (2016), Burton et al. (2012), Fernandez and Tsimring (2014), Ly (2014), Ly and Ermentrout (2011).

A long standing problem in neuroscience is to recover the network structure in a coupled system. This could for example be to infer the functional connectivity between units in a network of neurons from multiple extracellularly recorded spike trains, or how traces of EEG signals from different locations on the scalp affect each other, which we will treat in this paper. To the authors knowledge, this challenge is still lacking a sound statistical framework to model and test for interaction in a system, as well as impose statistical hypotheses on the network structure. For this task, cointegration analysis offers a refined statistical toolbox, where detailed information on the connections can be inferred, such as the direction and proportional strength of the coupling. The theory of cointegration was originally conceived by Granger (1981), and has since then also been the subject of intense research, most notably within the field of econometrics. In the monograph by Johansen (1996), the full likelihood theory for linear cointegration models with Gaussian i.i.d. errors is derived, and a framework for estimation and inference on parameters using the quotient test is presented. This well acknowledged framework is popularly termed the Johansen procedure. Even though cointegration analysis has developed from within the field of econometrics, it may potentially be used for different models outside economics, such as biological models in continuous time as we explore here. It has also been applied in climate analysis, see Schmith et al. (2012).

In this paper, we demonstrate how to apply cointegration analysis to a system of linearly phase coupled oscillating processes. To display the applicability of the method, we present a simulation experiment, where we present a statistical analysis of phase coupled systems with varying network structures, including uni-/ bi-directional and all-to-all couplings. We show that we can identify the proportional coupling strengths and directions given by the estimated cointegration matrix parameter. Our work is inspired by Dahlhaus and Neddermeyer (2012), which also introduces cointegration analysis as a statistical toolbox to neuroscientists and new challenges for researchers in cointegration theory. However, in contrast to Dahlhaus and Neddermeyer (2012), we incorporate the fact that we are dealing with continuous systems and also ensure that the cointegration property of the system is well posed as a linear structure. This approach assures that the conclusion on the interaction in the data is accurate in terms of cointegration.

The paper is composed as follows. In Sect. 2 we define a class of phase coupled oscillators, in Sect. 3 we highlight some cointegration theory for the analysis including an extension to discretely observed, continuous time models. In Sect. 4 we present a statistical analysis of linearly phase coupled oscillating systems and in Sect. 5 we analyze EEG recordings from an epileptic subject experiencing a seizure, previously analyzed by Shoeb (2009). We discuss the model and findings, conclude on the research and give an outlook of the future direction of the research in Sect. 6. Technical details are presented in the appendix.

Throughout we use the following notation and conventions: unless explicitly stated otherwise, time is assumed continuous, and the process is assumed observed with corresponding polar coordinates . Here denotes transposition. For a matrix M, with , we denote the orthogonal complement , a matrix such that (zero matrix). Also denote by the subspace spanned by the columns of a matrix A, and let denote the rank of the matrix, i.e., the dimension of .

Oscillating systems

Studying biological rhythms corresponds to studying systems of periodical processes. Intuitively we define a single oscillator as a continuous time bi-variate process , , such that revolve around some arbitrary center. Such a process can be derived from an equivalent process in polar coordinates , where is the phase process and is the amplitude process, such thatWe then define the process to be an oscillator if the phase process has a monotonic trend.

Defining a class of coupled oscillators

Definition (1) naturally extends to a system of coupled stochastic oscillators, where we observe p oscillators that interact, i.e., . Define a class of oscillators with phase () and amplitude () processes given by the multivariate stochastic differential equations (SDE) where are real valued vector functions, possibly depending on both or constant, are multivariate standard Wiener processes and such that is a positive semi-definite covariance matrix for . Assume the properties of (2) and (3) are such thatandwhere denotes the mean. Since are interpreted as the amplitudes of the individual oscillators, Eq. (4) is a natural assumption and Eq. (5) ensures that the individual oscillators actually revolve (anti-clockwise) around the center and that they are not ”stuck” in some part of the phase space, i.e., their angular velocities are positive. Note that we have defined the phase-trend as positive, corresponding to counter-clockwise rotation in accordance with the standard interpretation of the phase. However, for a negative trending process, one can either look at or simply interpret rotations as clockwise.

To emphasize the implication of inducing interaction in a system, for the data generating process (DGP) in the xy-plane, we derive a DGP from (2)–(3), see “Appendix 1”. Assuming that and we find thatHence, with the definitions (2)–(5) we have introduced a general class of coupled oscillators, where the specifications of f and g define the properties of the system, such as interaction. This broad definition of oscillating systems covers among others the Kuramoto model, see example (Sect. 2.5) below and other standard oscillators such as the FitzHugh–Nagumo and the Duffing oscillator. In this paper we will analyze phase coupled oscillators, and therefore we assume that , such that there is no feedback from the phase process into the amplitude and the k’th amplitude is not dependent on the rest. Hence, interaction in the system is solely through , such that the phase processes are attracted by some interdependent relation.

Linear coupling

The arbitrary function f enables us to choose any transformation of the variables to obtain a coupled system, including unidirectional coupling between phases or periodic forcing of the system if we extend f to depend on t as well, intermittent synchronization dependent on a threshold in process differences, etc.

Studying the general case where is nonlinear in and is a complex exercise. In this paper we restrict ourselves to models where f is composed of a linear mapping of and a function of , with components,for a real matrix and constant vector . With this restriction, the interaction between oscillators is linear in the phase, and the k’th oscillator is only dependent on the intrinsic amplitude through . We will refer to such a system as linearly phase coupled.

Although we impose the linear restriction on the interaction between phases, we can still model a broad set of coupling structures as we show with examples below. Since the interaction is given by , we note that the coupling strength in the system is given as the absolute values of the entries of and that row k of define how oscillator k depends on the rest. Note also that defines the attracting state for the phase relations, see example (Sect. 2.3) below. Normally is restricted to a constant, but in Sect. 4 we will relax this and investigate systems where is only approximately linear and has a sufficiently low variance. This implies a misspecified model, but as we will show, we can still identify the coupling structure, although inference on itself is less meaningful.

Example: Linearly phase coupled system with a degenerate process

Let f be defined as in (7) and assume that is a constant (positive) process such that . Then f is of the formwhere are constant vectors. For reduced rank matrices (2) is a continuous time cointegrated process (see Sect. 3) and f admits a linearly phase coupled system with intrinsic rotating frequencies . Note that if then there is no interaction in the system, and the individual oscillators will rotate according to their own , and we refer to the system as independent.

The linear specification implies that at most one attracting point can exist. As an illustration of this, assume a system composed of two coupled oscillators, withwhere . Since define an attracting state of the phase difference , then with the system is attracted towards being in-phase, whereas would imply that the system is attracted towards being in anti-phase. Considering that neither or depend on time, the system cannot switch to a different attracting regime.

To illustrate possible coupling structures, consider again the system of two oscillators and assume that . Then with and the coupling between is uni-directional where the arrow denote the direction of interaction. Likewise, if and then . However, if both then and the coupling is bi-directional. In general, if appears in the expression for oscillator , then . If the opposite is true, then and if both directions exist, then . For oscillator k is (oneway) independent from the rest, but it can still possibly influence others.

For systems where is a degenerate process, then and . With then (6) simplifies towhere . Note that if then (8) is simply a constant trend and hence (9) is a rotating process. One can show that the eigenvalues of the deterministic drift matrix in (9) in this case are complex conjugates, , where , implying that the solutions to (9) oscillate for . The oscillations are damped by the negative real part, but sustained by the noise term.

When is a constant vector process the properties of the system are fully identified by (2). Furthermore, if the noise level of the phases is sufficiently small, we can use the Hilbert transform1 to derive the phase process from observations of either or . This is a commonly used technique in signal processing and has been applied to oscillating systems as well, see Dahlhaus and Neddermeyer (2012), Pikovsky et al. (2001). For systems where is very noisy, this method is less applicable.

Example: Winfree oscillator

Let for a vector and for such thatWith these definitions (6) becomesThis example is taken from Winfree (2001) and extended with noise and phase interaction, and therefore we will refer to (10) as the (noisy) Winfree oscillator. Note that the formulation of implies that the amplitude fluctuates around . Due to this, we can for sufficiently small noise insist that for and therefore analyze the Winfree oscillator using the cointegration toolbox, assuming a constant in . In Sect. 4 we analyze the range of noise, , where the cointegration analysis still performs well.

Example: Kuramoto model

Choose such thatthen (2) is the Kuramoto model extended with a stochastic noise term, for phase coupled oscillators, where denotes the coupling strength between the k’th and j’th oscillators. In the classic version, , such that for a certain threshold , then with the oscillators exhibit synchronization. For an arbitrary process we cannot simplify (6), but with a degenerate we obtain the same expression as in (9) with as in (11).

For the Kuramoto model f is a nonlinear function, hence it is not directly applicable to a standard cointegration analysis where f is assumed linear. To emphasize this fact, consider the special case , where the Kuramoto model is particularly simple and (11) can be written explicitly as,where and . If at and the values of are large enough, then , such that and we can write a crude linear approximation of the sine function: , such thatThis is a coarse, but linear, approximation of the Kuramoto model and we can perform a cointegration analysis assuming this approximation is satisfactory. However, one must be cautious with this approximation. Consider , when . In this case , and henceand we see that not only do we add a term with , but the interaction also reverses sign. Recall that which implies a stationary relation in the system in (12), see Sect. 3.2. In (13) this condition is reversed, in the sense that will imply stationarity. If , (13) leads to an explosive system, which is not covered in this paper. Therefore, an essential requirement for an approximation of the Kuramoto model is a regime switching ability of (2). For a model with this property, we propose that cointegration analysis on a piecewise linear approximation of the Kuramoto model does make sense and can lead to correct conclusions regarding the network structure. In this paper we will not deal with non-linear cointegration of oscillating systems, but leave this direction open for future research. For a statistical analysis of nonlinear cointegrated systems of the form , i.e. time varying, or regime switching coefficients, see Bec and Rahbek (2004) and Kristensen and Rahbek (2013).

Note that with a general coupling constant , then the simple linear approximation to the Kuramoto model around is

Cointegration

Cointegration theory was originally developed for discrete time processes, however the ubiquitous use of continuous time models has inspired development of continuous time cointegration theory, see Kessler and Rahbek (2004, 2001). In order to present cointegration analysis as a framework for phase-processes, we therefore review some background on discrete time processes before entering into continuous time cointegrated models. The first part of this section is based on Johansen (1996) and Ltkepohl (2005).

Integrated process

Assume that is a discrete time vector autoregressive process,where , is a Gaussian white noise and is a deterministic term. The characteristic polynomial for (15) is the determinant of for , where is the p-dimensional identity matrix. If the roots of the characteristic polynomial are all outside the unit circle, then the initial values of can be given a distribution such that is stationary, see Johansen (1996).

If the characteristic polynomial of (15) contains one or more roots at , then there is no stationary solution of , and we say that the process is integrated. In particular, see Johansen (1996), will have reduced rank and can be written as with matrices of rank r. Moreover, the process is integrated of order one, I(1) with r cointegrating relations under regularity conditions presented in Sect. 3.2. Note that the order of integration is a stochastic property and hence including deterministic terms in a model does not change the order of integration.

In this paper we will only deal with I(1) processes, so when we refer to as integrated, we implicitly mean that is integrated of order 1.

Cointegrated process

Let and rewrite (15) with asAs already noted if implies then has a stationary representation (as an I(0) process). In particular, P has full rank p and all linear combinations of are stationary. If the -dimensional matrix P has reduced rank then with dimensional matrices of rank r. Moreover, the process is integrated of order one, I(1) with r cointegrating stationary relations provided with denoting the spectral radius. This we refer to as the I(1) conditions in the following.

Note that if the process is I(1) with no cointegration, while if (and ) then is I(0), or p stationary linear combinations exist. Under the reduced rank r, the system is written as,with b containing the r cointegration vectors and a the loadings or adjustment coefficients. Note that the entries of a and b are not uniquely identified, since we can use any non-singular transformation to obtain similar results. Rather we identify the subspaces , that is, the subspaces spanned by the columns of a, b, where we use the normalizationof b in order to identify parameters uniquely. Furthermore, let denote the matrix such that is orthogonal to , then a necessary condition for an I(1) process is that . For more on estimation and inference in cointegration models, see “Appendix 2”.

Continuous time cointegrated models

Kessler and Rahbek (2001, 2004) derive a cointegration theory for continuous time models, and conclude that for a discretely observed process, using conventional methods for discrete time generally apply to inference on continuous time parameters. Consider (2) with f as in (8) and for simplicity . This is a p-dimensional Ornstein–Uhlenbeck process. The exact solution isNote that for the solution (17) to be stationary, then must be a full rank matrix, and all eigenvalues must have a strictly negative real part. This implies that if is not of full rank, then is necessarily not stationary.

Assuming discrete observations of (17) at equidistant timepoints with timestep , the corresponding vector autoregressive process issuch that the difference process can be written aswhere andResults (18) and (19) hold in general for multivariate processes. Thus, to obtain an estimate for the continuous time matrix, , from the discrete time estimate , a logarithmic transformation involving is requiredFor a univariate process (20) is unique, however this is not the case for a multivariate process, due to the non-uniqueness of the multivariate logarithm. Because of this, we cannot uniquely identify , even though we have a unique estimate .

For a continuous time process , however, Kessler and Rahbek (2001, 2004) conclude that this is cointegrated if and only if the discretely observed process (18) is cointegrated. In this case P is of reduced rank, and can be decomposed with of full rank . However, it also holds thatfor a non-singular matrix and matrices , such that given weak conditions on the sampling timestep (see below), the following relations holdsee Kessler and Rahbek (2001, 2004). Hence, for continuous time cointegrated processes, we can infer on the number of cointegration relations () from discrete time observations, and also identify the subspaces spanned by the columns of and . Note however that due to the unidentified scaling , we can only identify the subspaces, but not the parameters themselves. They are only unique up to a scaling (), even though we have imposed the normalization (23) and thus uniquely identified a and b.

In the numerical part, we will refer to estimates of and , implicitly referring to the discrete time estimates. In terms of subspaces, there is no difference between the discrete and continuous time, but in order to interpret the continuous time matrix, one must translate the discrete estimate to a continuous estimate using (19).

It is important to note that when working with continuous time models, one must be careful with regard to the relation (19) between discrete and continuous time and the sampling timestep . Kessler and Rahbek (2004) refer to this issue as the embedding problem, and to ensure that the continuous time model is appropriate, one must check for in (18) that it is non-singular, i.e., , and that it has no negative eigenvalues. If this is the case and the underlying process is in fact cointegrated, the results above hold.

Likelihood ratio test for

Consider discrete observations from the continuous process (17) and denote by the hypothesis for . Then the set of hypotheses is nested,and correspond to the unrestricted model. The likelihood ratio test (LRT) that compare and is applied sequentially for and continued until against cannot be rejected, and thus determine the number of cointegrating relations for . The LRT statistic is given bywhere are the solutions to the eigenvalue problem (49), see “Appendix 2”. The asymptotic distribution of (22) is non-standard and therefore it must be simulated. Here, to also improve on small-sample performance we perform bootstrap simulations as presented by Cavaliere et al. (2012), in order to determine critical values. Specifically, given the data bootstrap sequences for are simulated and for each sequence the LRT statistic is re-computed. The empirical quantiles of are then used for testing. With r determined, is given by the r eigenvectors corresponding to and the parameter estimates follow by ordinary least squares estimation as outlined in “Appendix 2”.

Inference for and

Since we identify subspaces for and , then a normalization is necessary to identify the parameters uniquely. If is known, then follows by OLS. Hence, if we impose a normalization on , we can identify all parameters. A common normalization, see Johansen (1996), iswhere is a matrix and is any version of the r eigenvectors corresponding to the r largest eigenvalues. This ensures thatExtending the idea of normalization to restrictions for , we can impose such under the hypothesis . Assume that and that the parameters and are all unrestricted within their corresponding subspaces, except for normalization (23). Possible hypotheses for are linear restrictions as given bywhere , , , . The known matrices A and B represent the linear hypotheses and and are parameters to be estimated. It is also possible to combine the hypotheses for and and we denote this .

As an example, assume a system of 3 oscillators with . If we believe that is independent of and , we can specify the hypothesissuch thatThis restriction imply that and do not contribute to the dynamics of , and hence that the latter is independent.

If we want to investigate a possible 1:1 coupling between and , we can specifyand obtainNote however, that under the interaction between and also influence if . Hence, the system admits the relations , and , where the restriction implies that the last two relations are uni-directional.

If we believe that and are bi-directionally coupled, , but is independent and does not contribute to either nor , we can phrase this hypothesis as a combination of (24) and (25). This leads to the restricted matrixOther hypotheses, such as equal or proportional coupling strength or l : n coupling, can be specified by appropriately designing the matrices A and B. Thus, a broad variety of linear hypotheses on the parameter can be investigated, notably inference on the coupling directions and the effect of system disequilibrium on individual oscillators.

Evaluation of the hypotheses , and all lead to similar likelihood ratio tests. To calculate the test statistic, solve again the eigenvalue problem (49) for the unrestricted model, and dependent on the restrictions A and/or B obtain eigenvalues for the restricted model. The LRT statistic is then given bywhere is a generic substitute for any of . Each of these statistics has an asymptotic distribution with varying degrees of freedom (df),where m and s are the column dimensions of the matrices A and B, respectively. This shows that once is determined, statistical inference for and can be carried out, relatively straightforward. As for the rank determination, an alternative to the approximation for inference on and is to perform bootstrapping for the test (26), see Boswijk et al. (2016).

Numerical simulations

General setup

We perform a series of experiments with a system of linearly coupled Winfree oscillators such that and . Hence, for , we have a DGP withSince we examine simulations from the Winfree oscillator, our cointegration model will be misspecified since the amplitude is not deterministic and linear, but rather stochastic and fluctuating. However, since the amplitude of the Winfree oscillator has a relatively steady level (of course this also depends on the noise level), due to the squared multiplicative term in the amplitude process, we can approximate it as a constant. Hence we will do so in terms of analyzing the phase process as a cointegrating system. This also implies in terms of parameter estimation for the phase process, the estimate of the constant is a pseudo estimate of the parameter for the amplitude process, and hence we will compare the estimates to the true value of .

For each experiment we simulate 1.000.000 iterations of the oscillator (10) using the Euler–Maruyama scheme with timestep and then subsample for , thus obtaining (equidistant and discrete) observations of for . Subsampling every 5000th values diminishes the discretization error of the simulation scheme.

We use the same initial conditions,and parametersfor all the experiments so that the only varying parameter is the coupling structure.

Note that the parameter for is set equal to to obtain similar simulated outcomes for some experiments to investigate whether we can distinguish between interaction and independence between these two. We set the cointegration parameters for each experiment individually to impose different coupling structures, and will refer to the relevant model by it’s matrix, where k defines the model structure (see Fig. 1).

Fig. 1

Graphical representation of the four systems, represented by the matrix. The arrows define the direction of interaction, hence implies that is influencing (uni-directional coupling), and denotes bi-directional coupling, i.e. influence eachother

The discrete time model fitted to the data is specified aswhere the estimate is used to obtain through (20). The reported estimate for is scaled by the timestep . Note that is not time-dependent and hence this model will fit a constant parameter to a varying quantity and thus it is misspecified as mentioned above. Model (29) is estimated for all 4 systems of three oscillators and we report the parameters and for each system. The latter is compared to , which is the level parameter of the process.

In addition to a cointegration analysis we apply the mean phase coherence measure, see Mormann et al. (2000), bilaterally to the wrapped phases (i.e., for )as an initial measure of synchronization between the phases in the system. If this implies synchronization ( means that oscillators are perfectly phase locked). On the contrary, implies that the distribution of phase differences is approximately uniformly distributed on the unit circle. Note that the mean phase coherence measure is symmetrical, like correlation, and therefore it cannot reveal uni-directional coupling. In order to determine the significance of the R measures, we bootstrapped critical values for the hypothesis . Hence, these values are the same for all experiments and presented along with the measured R values. We compare the resulting value of R to the conclusion of the cointegration analysis.

We use the same seed for all experiments so that the outcomes are fully comparable in terms of stochastic noise dW. First we run a simulation with uncoupled oscillators as a benchmark, and then continue with coupled systems as presented in Fig. 1. Figure 2 display the x-coordinates for from a simulation of these four systems.

Fig. 2

50 observations (x-coordinates only) from numerical simulation of the Winfree oscillator: the column displays the independent model (31), the column displays the uni-directional coupled model (32), the column displays the bi-directional coupled model (33) and the column displays the fully coupled model (34)

The data analysis is carried out using the free software package R (R Core Team 2015). The source code for simulation and bootstrapping procedures are written in C++ to decrease the runtime, utilizing the interface package Rcpp for R and linear algebra package RcppArmadillo for C++. The source code is available in the package cods as supplementary material.

Independent oscillators

This experiment is used as a reference example. We setso and there is no interaction in the system.

Simulating the model and unwrapping the phases, we obtain the top-left plot of Fig. 3.

Fig. 3

Unwrapped phase processes from numerical simulation of the Winfree oscillator: displays the independent model (31), displays the uni-directional coupled model (32), displays the bi-directional coupled model (33) and displays the fully coupled model (34). The dotted lines represent the corresponding phases from the independent model in

Visual inspection of the plot could lead to the conclusion that and are coupled, however the mean phase coherence measure R for the phases indicates that this is not the case.The distribution of the mean phase coherence measure is unknown, but can be approximated by bootstrapping for , that is for no synchronization present. 1000 bootstrap samples yield the reported 5% critical values in parentheses above 0.17, thus the mean phase coherence measure suggest no synchronization present, which is the case.

Performing now a rank test for the rank of in the system, we obtain the first part of Table 1.

Table 1

Rank tests for models with the selected models indicated in bold

Model	H_r	Test values	p value
Π_0	r=0	16.47	0.663
	r≤1	3.94	0.753
	r≤2	0.05	0.812
Π_1	r=0	118.39	0.000
	r≤1	4.30	0.568
	r≤2	0.00	0.958
Π_2	r=0	104.48	0.000
	r≤1	3.84	0.707
	r≤2	0.03	0.843
Π_3	r=0	157.81	0.000
	r≤1	63.82	0.000
	r≤2	0.00	0.947

The test values are given by Eq. (22) and p values are determined by bootstrapping

The test does not reject the hypothesis , thus suggesting that there is no cointegration present in the system. This in turn implies that the oscillators are independent in terms of synchronization, in accordance with the DGP for , and with the mean phase coherence measure.

Uni-directional coupling

In this experiment we analyze a system with a uni-directional coupling. Letsuch that , and we have the stationary relation . Since , then and are acting independently, whereas is influenced by . Hence, the only coupling is .

The unwrapped phases for the simulation of model are seen in the top-right of Fig. 3. The dashed lines indicate the independent phases from the top-left of Fig. 3, and we see that phases are equal to their independent versions, whereas we now clearly find that is attracted towards due to the coupling structure in the system.

Examining the mean phase coherence in Eq. (30) for the system (note that is equal to the value in the previous section),we find indications of some synchronization between the phases and in the system compared to in the independent model. The value is significant on a 5% level as seen by the reported critical values, whereas for and the reported values are not. However, the mean phase coherence measure does not recognize the uni-directional coupling as is the case here. Thus, it cannot distinguish between , and .

Results from the rank test are in the second part of Table 1. Here we see that is clearly rejected, whereas cannot be rejected with a p value of 0.568. This indicates the presence of a single cointegration relation, in accordance with the construction of the model.

Fitting the model with , we obtain the unrestricted MLE regression estimates in Table 2. The cointegration relations are close to their true values (approximately within 1 standard error), and both and are statistically insignificant. Moreover, the estimates of suggests a 1:1 coupling between and .

Table 2

Fitted model

Parameter	True value	Unrestricted estimates			Restricted α and β
		Estimate	SE	p value	Estimate	SE	p value
α_1	-0.5	-0.527	0.049	<0.001	-0.514	0.048	<0.001
α_2	0	-0.050	0.049	0.307	0
α_3	0	0.059	0.048	0.223	0
β_1	1	1			1
β_2	-1	-0.981			-1
β_3	0	-0.016			0
κ_1	0.75	0.765	0.076	<0.001	0.638	0.081	<0.001
κ_2	1	1.035	0.075	<0.001	1.063	0.080	<0.001
κ_3	1	1.119	0.074	<0.001	1.086	0.080	<0.001

Therefore, we perform a likelihood test for reducing the unrestricted model, with restrictions for both and so that A fix and B restricts to a 1:1 coupling. This yields the test statistic 3.617 which is distributed with 4 degrees of freedom and hence implies a p value of 0.460. Thus, we recover the true uni-directional coupling structure of the simulated phases. The fitted model is presented in the right of Table 2.

The conclusion is that we have successfully identified the coupling structure of uni-directional coupled phases in a three dimensional system, with two independent phases, and one dependent. Since is completely independent of and and we can discard when interpreting the cointegration in the system. Then we can interpret the cointegration parameter as the coupling strength and as the coupling scheme, here 1:1. If we had analyzed different data, with a estimate close to , we could then identify a n:1 coupling between and . This can be seen from the fact that in this case would be a stationary relation, and thus would rotate n times slower than .

A bi-directional coupling with one independent oscillator

We now look at a system withHence, and again we have 1 stationary relation , but now with only independent, and a bidirectional coupling .

Simulating the -model we obtain the bottom-left of Fig. 3. We have included the dashed lines again, as references for the independent system. If we contrast the bottom-left of Fig. 3 with the top-right of Fig. 3, we now find that and are attracting each other, and hence they are both different from their independent versions. Since , their coupling strength is equal, and the coupled phases lies roughly in the middle between the independent ones. If we look at the mean phase coherence measure for the pairwise interactions,we find relatively strong evidence of a coupling between the phases and , the value is higher than in the uni-directional case and it is (again) significant given the bootstrapped critical values. However, again we cannot distinguish between types of coupling structures.

Performing a rank test for cointegration in the system with , we see in the third part of Table 1 that is clearly rejected, and we find that the rank of is estimated to 1 with a p value of 0.707. Hence, we recover the correct dimension of the column space of , and fitting a model with yields the parameters in the left of Table 3.

Table 3

Fitted model

Parameter	True value	Unrestricted estimates			Restricted α and β			Restricted α and β with A^∗
		Estimate	SE	p value	Estimate	SE	p value	Estimate	SE	p value
α_1	-0.5	-0.530	0.071	<0.001	-0.506	0.069	<0.001	-0.475	0.069	<0.001
α_2	0.5	0.450	0.069	<0.001	0.443	0.067	<0.001	0.475	0.067	<0.001
α_3	0	0.087	0.070	0.214	0			0
β_1	1	1			1			1
β_2	-1	-0.970			-1			-1
β_3	0	-0.022			0			0
κ_1	0.75	0.754	0.072	<0.001	0.646	0.076	<0.001	0.660	0.076	<0.001
κ_2	1	0.943	0.070	<0.001	1.040	0.074	<0.001	1.053	0.074	<0.001
κ_3	1	1.103	0.071	<0.001	1.086	0.075	<0.001	1.086	0.075	<0.001

The only insignificant parameter for the model is , which is in accordance with the construction of the model. Specifying the hypothesisand performing a likelihood ratio test for the reduction yields a test statistic of 3.340, which follows a with 3 degrees of freedom, and result in a p value of 0.342. The fitted model is given in the middle of Table 3. If we instead of specifyimplying that , we obtain a test statistic of 3.880, with 4 degrees of freedom, and a p value of 0.423. Thus, we can also restrict the model to one where the coupling strengths are equal in magnitude. The fitted model is presented in the right part of Table 3.

Summing up, in a system of bi-directional coupled oscillators plus one independent, we can identify the correct coupling between them, including identifying the proportionally equal coupling strength between the coupled phases. Again we identify , and hence we can interpret the cointegration parameters as before, hence is the coupling strength, and the interaction, again 1:1 coupling.

Fully coupled system

We specify a system with full interaction between all phases.The and matrix are chosen, such thatinspired by the simplistic linearization of the Kuramoto model, as presented in Eq. (14). Note that .

The simulated phases are shown in the bottom-right of Fig. 3. Comparing to the dashed (independent) versions, we now find that all phases are different from their independent versions. It appears as if dominate the system, since is attracted closer to their independent versions than otherwise, but it is also a two-against one () scheme, and we roughly observe that is attracted 2/3 towards , whereas are attracted 1 / 3 towards . So by the construction of the system, this behavior seems natural. We find that the mean phase coherence measureindicates bilateral synchronization for all phases, and all values are significant. The rank test also gives clear evidence of cointegration and we identify , as seen in the bottom part of Table 1, where both the hypotheses and are rejected. Fitting a model with yields the left half of Table 4.

Table 4

Fitted model

Parameter	True value	Unrestricted estimates			Restricted α and β
		Estimate	SE	p value	Estimate	SE	p value
α_11	-0.50	-0.584	0.075	<0.001	-0.569	0.075	<0.001
α_21	0.25	0.232	0.073	<0.001	0.241	0.072	<0.001
α_31	0.25	0.326	0.072	<0.001	0.328	0.072	<0.001
α_12	0.25	0.223	0.067	<0.001	0.224	0.067	<0.001
α_22	-0.50	-0.423	0.064	<0.001	-0.423	0.064	<0.001
α_32	0.25	0.201	0.064	<0.001	0.199	0.064	<0.001
β_11	1	1			1
β_21	0	0			0
β_31	-1	-0.997			-1
β_12	0	0			0
β_22	1	1			1
β_32	-1	-0.999			-1
κ_1	0.75	0.712	0.076	<0.001	0.607	0.083	<0.001
κ_2	1	1.054	0.074	<0.001	1.061	0.080	<0.001
κ_2	1	1.023	0.073	<0.001	1.130	0.080	<0.001

The estimated ’s are close to their respective values, whereas some of the parameters deviate (more than 1 standard error) from their true values. If we inspect the estimated and compare with the true it looks better. The row sums are close to zero as they should be, and the signs are correct. The proportional coupling strengths are off though, especially between , but it seems that is relatively well estimated considering the identification issues. Recall that we can determine the subspaces and for continuous time cointegration models, see Kessler and Rahbek (2004), but that we have problems regarding the scaling of (see Sect. 3.3).

Inspired by the fitted values, we restrict both matrices and we find that the test statistic is 1.73, distributed with 4 degrees of freedom, and thus a p value of 0.785. Hence, we can reduce the model to one with restrictions that generates the true structure of . The estimated model parameters are presented in Table 4, and the corresponding isConcluding on the fully coupled system, we find that we can correctly identify the dimension of the cointegration relations. We can also determine the coupling structure as given by the parameters and . However, interpretation in this experiment is more informative in terms of , since with , the interpretation of cointegration parameters is not as intuitive as in the case of . We obtain a estimate that is reminiscent of the true matrix, with the true directions of the coupling, and strengths somewhat close to the actual values. Thus, we can interpret the system as fully coupled, in a simplistic (linear) Kuramoto type model.

Strength of coupling and identification of interaction

In this section, we compare the mean phase coherence measure to the cointegration analysis with respect to interactions in the system. More specifically, we look at how strong the coupling constants in must be in order for the two methods to conclude correctly on interaction in the system. We reuse the parameter settings (34) from the fully coupled experiment, but use a scaled matrix , for , where controls the coupling strength. The higher is, the stronger the coupling, and hence the attraction between phases. Note that corresponds to the model and corresponds to . The p values are calculated using bootstrapping as presented by Cavaliere et al. (2012) to obtain an estimate of the asymptotic distribution of the trace test statistics.

The aim is to investigate the rank test for varying compared to identification of interaction in the system using the mean phase coherence measure. Since low values of implies weak interaction, the expectation is that both methods will produce doubtful results in a low value regime. From the previous experiment on the fully coupled oscillators, the mean phase coherence measure produced low values on identifying the interaction of the system, hence we expect that the rank test will outperform for low values of .

The experimental setup is 100 repetitions for each value of , and in each repetition perform 500 bootstrap samples to estimate the p value for the hypotheses . Figure 4 presents the median p values for the rank test and median mean phase coherence measures against . The top row of the figure shows the p values for respectively, and the bottom row shows the mean phase coherence (R) measures for pairs of and . The dotted lines indicate the value, under which we reject the hypothesis. For the mean phase coherence measure, the 95% significance level for the hypothesis has been determined numerically using bootstrapping and is indicated by the dotted lines. If the R-measure falls below this line, independence cannot be rejected.

Fig. 4

Identification of interaction for varying coupling strengths for the model . Top row are hypotheses using the trace test. The orange bounds indicates empirical quantiles ranging from 2.5 to and the blue lines represents the median. The bottom row are mean phase coherence measures again with empirical quantiles illustrated. Dashed lines in the top row show the 5% confidence levels. Dashed lines in the bottom row show the 95% quantile for at , found by bootstrapping (color figure online)

Seen in the top row of Fig. 4, at least half the simulations reject for , and at least half the simulations reject for . The test does not reject for around 88% of the simulated samples for any values of . Thus, for , we can conclude that there is interaction present in the system, and in most of the simulations we also recognize the true .

If we turn to the bottom row of Fig. 4, where the mean phase coherence measures are shown, we find that half the simulations does not reject the hypothesis for and 0.35, respectively, for , and , thus clearly indicating an inferior detection of interaction for small values of equivalent to weak couplings.

Concluding on this experiment, we find that the rank test detects interaction in the system already at relatively weak coupling strengths. In contrast to this, the coupling must be significantly stronger for a sound conclusion on interaction in the system when using mean phase coherence as a measure of interaction. Furthermore, when detecting interaction in the system, the rank test is also very capable of identifying the true rank of the system, despite a misspecified model. Higher sample sizes will of course improve the inference results.

Consistency of the rank estimation

To investigate the consistency of the cointegration algorithm, we performed an experiment with 1000 repetitions of simulations for Winfree oscillators, the uni-directional coupling, the bi-directional and the fully coupled systems, respectively, and evaluating the rank test, using the same setup as in Sect. 4.1. Table 5 present the percentages of conclusions regarding hypotheses , for each model. Comparing with critical values at a level, obtained by bootstrapping, see Cavaliere et al. (2012), we find that comparing the percentage of simulations where the test correctly identifies the cointegration rank of 1 for uni- and bi-directional coupling are 76.8 and , respectively, at a 5% significance level. For the fully coupled system the percentage is , and for an independent system the percentage is .

Table 5

Percentage of conclusions on , at a 5% significance level for a sample size of 2000

Model	r=0	r≤1	r≤2	r≤3
Independent (%)	96.2	2.2	1.3	0.3
Uni-directional (%)	1.7	76.8	19.0	2.5
Bi-directional (%)	2.4	69.8	24.7	3.1
Fully coupled (%)	0.0	1.3	85.5	13.2

Note that the conclusion means that is of full rank and therefore invertible, hence . Correct conclusions in bold

These results show that identification of interaction in a system of coupled oscillators is quite precise, and the rank is underestimated in 2.5% of the simulations for any model. In the case of independent or full interaction, the method is very good, whereas for systems with directed interaction, or interaction among some oscillators the frequency of overestimating the rank is 20–25%. This discrepancy seems intuitively correct, since for the latter systems the true model is a subset of the model of higher rank. As before higher sample sizes will of course improve the inference results.

In Table 6 we compare, in percentages, the conclusions on interaction in the systems, for each model. The values for the rank test presented here, are the summed values from Table 5 for . We find that both methods are very adept in identifying interaction in these systems. The results, however, should be held against the previous section, where the rank test outperformed the mean phase coherence measure for weak coupling strength. Also noting the fact, that the mean phase coherence measure cannot account for uni-directional coupling, our overall conclusion is that in terms of identifying interaction in the system, the methods seem to perform equally well for stronger coupling, whereas in explaining the system architecture, a cointegration analysis leaves us with more information on how the network is constructed.

Table 6

Percentage of conclusions on interaction indicated by the rank test and the mean phase coherence measures, at a 5% significance level for a sample size of 2000

Model	Rank test	R(ϕ_1t,ϕ_2t)	R(ϕ_1t,ϕ_3t)	R(ϕ_2t,ϕ_3t)
Independent (%)	3.8	4.7	4.4	5.7
Uni-directional (%)	98.3	99.8	5.6	4.4
Bi-directional (%)	97.6	100.0	7.2	7.0
Fully coupled (%)	100	100.0	100	100.0

Analysis of EEG data

Electroencephalography (EEG) signals are recordings from electrodes distributed on the scalp of subjects. The recorded brainwave patterns are, among others, used for diagnosing sleep disorders, coma or epilepsy. A study on 22 subjects experiencing epileptic seizures from the Children’s Hospital Boston is presented by Shoeb (2009) with the aim of detecting seizures based on multiple hours of recordings for each individual. Figure 5 displays an EEG recording of a single subject during a period that include a seizure identified by Shoeb (2009) between 2996 and 3036 s. The seizure is marked by two red dashed lines in Fig. 5. The labels for the signals refer to the individual electrodes on the scalp. We analyze the four signals FP1-F7, FP1-F3, FP2-F4 and FP2-F8, where FP refer to the frontal lobes and F refer to a row of electrodes placed behind these. Even numbered electrodes are on the right side and odd numbered electrodes are on the left side. Smaller (larger) numberings imply that the electrode is placed closer to (further from) the center of the scalp. Hence FP1-F7, FP1-F3 are measurements from the left side, with F3 placed closer to the center than F7, and likewise for right side signals FP2-F4 and FP2-F8. The electrodes for these four signals mirror each other on the left/right side of the scalp. We analyze the seizure period of 40 s and the 40 s leading up to the seizure, i.e. we analyze the two intervals [2956; 2996] and [2996; 3036] respectively, and refer to these as prior to seizure and during seizure. With a sample frequency of 256 measurements each second there are a total of 10,240 measurements for each of the four signals during the 40 s intervals. For more details on the data, see Shoeb (2009). The objective is to compare two fitted cointegration models with interaction as in Eq. (8) for each period:discretely observed for in each of the two intervals.

Fig. 5

EEG recording leading op to a seizure and afterwards for a 11 year old female subject. The interval [2996;3036] s, as indicated by the vertical red dashed lines, is defined by Shoeb (2009) as a seizure. We analyze the four blue signals, FP1-F7, FP1-F3, FP2-F4 and FP2-F8 (color figure online)

The phase processes of the four signals are estimated using the Hilbert transform (see Sect. 2.3). Figure 6 shows the four signals in the two periods and their corresponding estimated unwrapped phase processes. Hence the offsets are in for the individual phase processes in each period. If we had not split the measurements at 2996 s, the phases in the bottom right of Fig. 6 would be continuations of the phases in the bottom left. A visual inspection of Fig. 6 shows that when transitioning to the seizure period, the phases change to a slower pace (the slopes decrease). Also, prior to the seizure all four phases are closer with no clear distinction between right side and left side phases. During the seizure, the phases split in two groups: right and left side respectively.

Fig. 6

FP1-F7, FP1-F3, FP2-F4 EEG signals and estimated phase processes for a 11 year old female subject. Top left EEG signals prior to a seizure. Top right EEG signals during a seizure. Bottom left estimated phase processes prior to a seizure. Bottom right estimated phase processes during a seizure

This indicates that the model regime changes when transitioning into the seizure period. Table 7 shows the mean phase coherence measures bilaterally for the 4 phase processes and the average of these. Comparing the columns we find no clear indication of a change in the phase regime when transitioning into the seizure period based on this measure, the average change is only 7.5%. However, the measure does indicate interaction in the system among all phases.

Table 7

Mean phase coherence measures for EEG phases prior to and during the seizure

	Prior to seizure	During seizure
R_(FP1-F3; FP2-F4)	0.480	0.542
R_(FP1-F3; FP1-F7)	0.535	0.644
R_(FP1-F3; FP2-F8)	0.295	0.184
R_(FP2-F4; FP1-F7)	0.321	0.350
R_(FP2-F4; FP2-F8)	0.486	0.342
R_(FP1-F7; FP2-F8)	0.525	0.379
Average	0.440	0.407

Table 8 displays the results of a rank test procedure for the system of the four EEG phase processes.

Table 8

Rank tests for EEG phases in the bottom of Fig. 6

H_r	Prior to seizure		During seizure
	Test values	p value	Test values	p value
r=0	105.87	0.000	1132.64	0.000
r≤1	42.82	0.000	41.68	0.008
r≤2	9.98	0.053	7.19	0.618
r≤3	0.46	0.439	0.72	0.786

The rank is determined to in both periods, although the conclusion is far stronger during the seizure. The significance of the statistics are found using 5000 bootstrap samples prior to the seizure due the border limit case of around 5%, during the seizure the p value is determined from 2000 bootstrap samples

In accordance with the indications from the mean phase coherence measure, the conclusion is a clear presence of cointegration during both periods. Prior to the seizure the rank test of is close to the usual 5% significance level, hence the p value here is determined using 5000 bootstrap samples, in contrast to the 2000 bootstrap samples used in the other interval, as the conclusion here is quite clear with a p value 0.62. In both cases we choose the rank for the system.

The fitted models are presented in Table 9 with the model fit prior to the seizure on the left side and the fit during the seizure on the right side. If we first note the estimated ’s, these are larger during the seizure and significantly so for FP1-F3 and FP2-F4, implying that these phase processes exhibit significantly higher intrinsic linear trends during the seizure. On the other hand, directly interpreting the cointegration parameters is not clear. Recall that these parameters specify subspaces, in this case within . We therefore look at the estimated matrices in Table 10 to compare the models for each period.

Table 9

Fitted model for EEG phases F7-T7, T7-P7 and FP1-F7

Parameter	Prior to seizure			During seizure
	Estimate	SE	p value	Estimate	SE	p value
α_FP1-F3,1	-0.100	0.018	<0.001	-0.462	0.028	<0.001
α_FP1-F7,1	-0.002	0.019	0.930	-0.308	0.032	<0.001
α_FP2-F4,1	-0.035	0.017	0.044	-0.722	0.035	<0.001
α_FP2-F8,1	-0.115	0.030	<0.001	-0.648	0.042	<0.001
α_FP1-F3,2	-0.117	0.016	<0.001	0.041	0.033	0.212
α_FP1-F7,2	-0.024	0.016	0.147	0.071	0.037	0.057
α_FP2-F4,2	-0.026	0.015	0.084	0.173	0.041	<0.001
α_FP2-F8,2	-0.049	0.026	0.063	0.468	0.049	<0.001
β_FP1-F3,1	1			1
β_FP1-F7,1	0			0
β_FP2-F4,1	-3.424			-0.036
β_FP2-F8,1	2.610			-0.573
β_FP1-F3,2	0			0
β_FP1-F7,2	1			1
β_FP2-F4,2	2.486			-0.840
β_FP2-F8,2	-3.631			0.188
μ_FP1-F3	25.210	2.162	<0.001	39.647	1.307	<0.001
μ_FP1-F7	30.648	2.252	<0.001	36.499	1.473	<0.001
μ_FP2-F4	39.058	2.107	<0.001	58.268	1.608	<0.001
μ_FP2-F8	48.853	3.615	<0.001	54.765	1.947	<0.001

Table 10

Fitted matrices for the two periods

	Π^ prior to seizure				Π^ during seizure
	FP1-F3	FP1-F7	FP2-F4	FP2-F8	FP1-F3	FP1-F7	FP2-F4	FP2-F8
FP1-F3	4.388	1.572	-11.120	5.743	-5.305	-11.021	9.447	0.971
FP1-F7	1.519	0.892	-2.985	0.725	-4.335	-7.285	6.275	1.116
FP2-F4	0.540	-0.050	-1.971	1.589	-10.265	-17.047	14.686	2.681
FP2-F8	-0.733	-1.658	-1.613	4.108	-14.907	-14.729	12.909	5.776

On the left side is the estimated matrix prior to the seizure, on the right side is the estimated matrix during the seizure

Here we can determine an all-to-all coupling during both periods and the estimated cointegration matrices show a clear difference for the two intervals. Prior to the seizure the right side signals FP2-F4 and FP2-F8 are much less influenced by the feedback in the system, whereas during the seizure both experience a much larger feedback from the left side signals FP1-F3 and FP1-F7 respectively. Surprisingly, the FP2-F8 signal does not seem to impose a large influence in the system in either interval. It is also interesting to note the changing signs in the two matrices. The two left side signals exhibit a positive feedback on themselves prior to the seizure, whereas during the seizure they impose a negative feedback both on themselves and the right side signals. This could possibly be part of an explanation of the slight kink seen in the phases around 3015–3020 s halfway through the seizure.

Concluding on this analysis we find, not surprisingly, a fully coupled 4 dimensional system with a clear change in the trends prior to and during the seizure. We find that during the seizure the interaction in the system is much stronger, suggesting the more distinctive phases shown in this interval. Including this temporal effect into a single cointegration model covering the full period by utilizing regime switching cointegration models, would be an interesting pursuit for future work.

Discussion

In this paper we have investigated the use of cointegration analysis to determine coupling structures in linearly phase coupled systems. Using these techniques we can with a good precision identify the coupling structure as a subspace for this type of model. A standard measure to identify synchronization in the literature is the mean phase coherence measure. Contrary to this standard measure, we can detect uni-directional coupling, and we can construct and test hypotheses on the model in form of linear restrictions in the estimated subspace. Furthermore, comparing the mean phase coherence measure with the cointegration analysis in Sect. 4.6, we found that cointegration detects interaction in a system more robustly and for weaker coupling strength than does the mean phase coherence measure. Combined with the fact that cointegration does not just provide a level of synchronization, but rather the structure of the synchronization mechanism, this technique can be used to infer system structures in a much more detailed manner. Of course this higher level of information comes at a cost, since the mean phase coherence measure is easily implemented for any system, whereas the cointegration analysis is more involved and time consuming.

Due to the linear nature of the cointegration theory used, we are not able to cover more complex models, such as the Kuramoto model. Thus, an important extension for future work would be to allow for nonlinear coupling functions. However, the linear structure appears naturally when considering a linearization around some phase-locked state, such as for systems showing synchrony or asynchrony. Another interesting pursuit is to extend the model framework to include nonlinear deterministic trends, such that also models like the FitzHugh–Nagumo or the van der Pol oscillator would be covered. The model considered in this paper was constructed from the starting point of the phase process in the spirit of the Kuramoto model, and noise was added on this level. Another approach would be to start from a biological model or a reduction thereof and introduce the noise on the DGP. This would also lead to non-linearities both in drift and diffusion of the phase process. Finally, high dimensional systems are a major challenge in the area of coupled oscillators, hence it would only be natural to investigate cointegration properties of high dimensional systems. A system of more than two synchronizing oscillators that are nonlinearly phase coupled, facilitate chaotic behavior since phases can then bilaterally attract and repel each other. When the number of oscillators increase, one quickly ends up with intuitive shortcomings. The number of parameters rapidly increase with the dimension of the system, possibly leading to a desirable reduction to a sparse interaction structure. This is a key issue with the cointegration framework, which take into account all individual oscillators, as opposed to a mean-field approach that does not run into the same curse of dimensionality. The quality of the estimators will rapidly decrease with increasing dimension of the parameter space or numerical problems may arise. This problem might be alleviated by imposing a sparse interaction structure through a LASSO penalization.

Cointegration to identify coupling of oscillators has been attempted before in a neuroscience context by Dahlhaus and Neddermeyer (2012). There, the Kuramoto model is approximated for strongly phase coupled oscillators by setting , since the phase differences are assumed to be small. We have used the idea from Dahlhaus and Neddermeyer (2012) of analyzing the unwrapped multivariate phase process. Contrary to Dahlhaus and Neddermeyer (2012), however, we have not linearized the sine function to replicate Kuramoto, since this will cause a discrepancy when the phase difference of two oscillators is closer to than 0 (or ). To mitigate this problem, we have instead taken the approach of designing a DGP with the properties we are interested in, and which allows for any phase differences. Furthermore, this DGP enables us to specify a cointegration model that comply with data from this DGP. Although it may not fully comply with a biological model, it can point to where necessary flexibility is needed in order to develop more realistic cointegration models for biological processes. A first attempt to analyze EEG signals with cointegration analysis with linear coupling structures has been presented. The results are promising, and reveal a finer dependence structure characterizing states of seizure and non-seizure in epileptic patients, which in this example was not possible from the simple Mean Phase Coherence measure. To fully explore the potential of the cointegration analysis for EEG signals, it would be useful to extend the model and analysis tools to allow for non-linearities and simultaneous treatment of many traces, as well as time varying coupling strengths.

Summing up, by applying cointegration as a technique to the field of coupled oscillators in biology, we open up for a whole new area of applications for this statistical theory. On the other hand, using cointegration methods, biologists can gain new insights into network structures, being able to fit models and carry out statistical hypothesis testing. If the cointegration framework presented in this paper can be extended to include the standard models currently used in the field, cointegration would prove a powerful analysis tool for researchers.