A practical method for estimating coupling functions in complex dynamical systems
aa r X i v : . [ n li n . C D ] O c t rsta.royalsocietypublishing.org Research
Article submitted to journal
Subject Areas:
Nonlinear dynamics, coupledoscillators, data analysis
Keywords:
Coupling function, phase dynamics,parameter estimation
Author for correspondence:
Isao T, Tokudae-mail: [email protected]
A practical method forestimating coupling functionsin complex dynamicalsystems
Isao T. Tokuda , Zoran Levnajic , andKazuyoshi Ishimura Department of Mechanical Engineering, RitsumeikanUniversity, Kusatsu, Japan Complex Systems and Data Science Lab, Faculty ofInformation Studies in Novo Mesto, Novo Mesto,Slovenia
A foremost challenge in modern network scienceis the inverse problem of reconstruction (inference)of coupling equations and network topology fromthe measurements of the network dynamics. Ofparticular interest are the methods that can operateon real (empirical) data without interfering with thesystem. One of such earlier attempts [Tokuda et al. ,PRL, 2007] was a method suited for general limit-cycle oscillators, yielding both oscillators’ naturalfrequencies and coupling functions between them(phase equations) from empirically measured timeseries. The present paper reviews the above method ina way comprehensive to domain-scientists other thanphysics. It also presents applications of the method to(i) detection of the network connectivity, (ii) inferenceof the phase sensitivity function, (iii) approximationof the interaction among phase-coherent chaoticoscillators, (iv) experimental data from a forced Vander Pol electric circuit. This reaffirms the rangeof applicability of the method for reconstructingcoupling functions and makes it accessible to a muchwider scientific community.
1. Introduction
Complex networks are representations of complexsystems, where nodes (vertices) represent system’sunits and links (edges) represent the interactionsamong those units [1–4]. The functioning of a realnetwork is a cumulative effect of its structure (topologyof connections among nodes/units) and dynamics(interactions/relationships among these nodes) [3,4].Hence, in models of real networks, nodes are often c (cid:13) The Authors. Published by the Royal Society under the terms of theCreative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author andsource are credited. r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A .................................................................. assumed to be (simple) systems with their inherent dynamics, whereas links mediate thedynamical coupling between the connected pairs of nodes. Using this paradigm, network sciencehas made valuable contributions to all scientific disciplines that involve systems composed ofmany units, including biology, neuroscience, sociology, economics, etc. [1–6].To really grasp the functioning of a real network, we need information on both its structureand its dynamics. The inverse problem of reconstructing (or inferring) this information from theempirical data is a foremost challenge in modern network science. Namely, understanding theinner connectivity patterns of real networks not only enables us to grasp their operations, butalso helps in controlling and predicting their behavior [7–19].The problem of network reconstruction can be seen as composed of two parts. The first partis the reconstruction of network topology, where one tries to learn which pairs of nodes areconnected and which are not. This can (in some cases) be done separately from the second part ofthe problem, which is the reconstruction of the coupling functions that dictate how the connectednodes interact. Of course, two parts of the problem are inherently related, but which one to tacklefirst depends on what data are available, what assumptions can be reasonably made about thesystem, and what exactly we wish to learn.Numerous reconstruction methods have been proposed over the last decade, both in physics[7–21] and in computer science literature [16,22–30]. While some methods tackle only one of thetwo above mentioned parts of the problem [10], other methods seek to address both parts atthe same time. In physics literature, special emphasis is put on the methods aimed at oscillatorysystems as the most researched paradigm of collective dynamics. This includes methods focusedon either topology, coupling functions, or both [7,8,10,18].On a somewhat different front, research effort was devoted to the problem of estimation ofphase variables and phase equations from the data [7,11,31–38]. Namely, isolated units in manyreal systems exhibit oscillatory nature, in the sense that they can be well approximated as limit-cycle oscillators (oscillator whose dynamics after transients reduces to periodic or quasi–periodicorbit). Researchers showed that, even if the oscillatory behavior is very stochastic, there are robustways to extract a well defined phase variable for each network node, and hence reconstructthe phase equations that describe the system dynamics. This paradigm found applications indiverse domain sciences, notably biology and neuroscience, where many systems have thisnature. Estimating phase equations, however, is nothing but reconstructing coupling functionsfrom data. While such a reconstruction approach is valid only in the approximation of phasevariables, these methods require very little additional assumptions about the system. This meansthey can be almost immediately applied to empirical data [7,31,36,39–41].For a system of phase equations, a standard way to construct the coupling function is tomeasure phase sensitivity function of an individual oscillator element and obtain the couplingfunction by averaging method that computes the amount of phase shift induced throughinteraction with another oscillator element [42]. However, a precisely measured phase sensitivityfunction is not always accessible, since it requires application of external perturbations to anindividual oscillator, which cannot always be isolated from the rest of the system [43–52].On the other hand, as a non-invasive approach, the coupling function can be inferred directlyfrom time-trace data measured from coupled oscillators [7,31–36,38,40,41]. One of them is amethod by Tokuda et al. [31]: this approach utilized a multiple shooting method to realizerobust parameter estimation of the coupling functions. The multiple shooting provides a generalframework for fitting ordinary differential equations to recorded time-trace data. It is applicableto any system, where the dynamics of individual nodes can be approximated as those of limit-cycle oscillators, yielding both oscillators’ natural frequencies and coupling functions betweenthem (phase equations). Most importantly, the method was actually shown to operate very wellwith the data from a real experiment, which highlights its potential for practical use for physicsproblems and otherwise [31,40,53].The contribution of the present paper is two-fold. First, we review this method in a way thatis understandable and approachable to communities outside physics. With this, we hope to make r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A .................................................................. our method more useful to field such as biology and neuroscience, for which it was originallyintended. Second, we show and discuss applications of this method, specifically: (i) we utilize theestimated coupling function for detecting the connectivity of oscillator networks, (ii) the methodis extended to inference of the phase sensitivity function, which is vital for phase equations, (iii)the coupling function is estimated for coupled chaotic oscillators to demonstrate how well thephase model approximates chaotic phase synchronization, (iv) using an experimental data froma system of Van der Pol electric circuits, we show how the method can be applied to real data.The rest of the paper is organized as follows. In the next section, we review the originalmethod in a comprehensive way. In section 3, we discuss the problem of inferring the networkconnectivities. In section 4, we present further applications mentioned above. In the last section,we discuss our findings and lay out perspectives for future work.
2. The Original Method
In this section, we describe the original method in a more comprehensive way than the originalliterature [31] and show how it works for the case of coupled FitzHugh-Nagumo oscillators. (a) Multiple Shooting Method
Our approach is based on the multiple shooting method, which has been developed in physicsand engineering to provide a general framework for fitting ordinary differential equationsto recorded time series [54]. The methodology is applicable to a situation, where the systemequations are known a priori . When the equations and the recorded data are in a good quantitativeagreement, unknown parameters of the system can be precisely estimated as follows.We consider a nonlinear system ˙ x = F ( p , x ) , (2.1)where x , p and F represent state variables, parameters, and autonomous dynamics of the system,respectively. The system may generate nonlinear dynamics such as limit cycles, torus, strangeattractors, and transient dynamics to these attractors. The equation (2.1) may describe a varietyof systems of interest in science and engineering such as electric circuits and lasers. Empiricaldata consists of oscillators’ states measured as { x ( n∆t )) : n = 1 , · · · , M } ( ∆t : sampling time, M :data points). This corresponds to an experimental situation, in which the system state ( e.g. ,current and voltage of electric circuits, laser, etc. ) is fully recorded. Then the parameter values p that underly the measurement data can be estimated by fitting the original equations (2.1) tothe recording data. First, time evolution of the original equations (2.1) starting from an initialcondition x (0) is denoted by φ t ( x (0) , p ) . Then, at each sampling time t = i∆t , the equations mustsatisfy the boundary conditions: x (( n + 1) ∆t ) = φ ∆t ( x ( n∆t ) , p ) . With respect to the unknownparameters p , the nonlinear equations are solved by the generalized Newton method [55]. Tocompute the gradients ∂φ i /∂p , which are needed for the Newton method, variational equations ddt ( ∂φ i ∂p j ) = ∂f i ∂p j + P Nk =1 ∂f i ∂φ k ∂φ k ∂p j are solved simultaneously, where f i represents i th equation ofthe original dynamics (2.1) as ˙ x i = f i ( x , p ) . The evolution function φ t as well as the variationalequations are integrated numerically, using whichever integration scheme ( e.g. , 4th–order Runge-Kutta ). It has been shown that, when the equations and the experimental data are in goodquantitative agreement, the unknown parameters can be precisely estimated for real–worldsystems including electric circuits and lasers. All steps in the above computational procedurecan be realized relatively easily with standard programming knowledge. (b) Problem and Method
Equipped with the knowledge of multiple shooting method, we now explain in detail how it canbe utilized for inferring the coupling functions. We begin by considering a network composed ofinteracting oscillator elements. In biology, such systems include a network of circadian cells in the r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A .................................................................. suprachiasmatic nucleus [56], brain network composed of many spiking neurons [43,46], cardiacmuscle cells in the heart [57], etc. In terms of nonlinear dynamics, such systems are described as asystem of N coupled limit cycle oscillators: ˙ x i = F i ( x i ) + CN N X j =1 ,j = i T i,j G ( x i , x j ) . (2.2)Here, x i and F i ( i = 1 , , · · · , N ) represent state variables and autonomous dynamics of the i thoscillator element, respectively. While G represents an interaction function between i th and j thoscillators, strength of their interaction is determined by the coupling constant C . The matrix { T i,j } describes connectivity between the oscillator elements. For simplicity, we suppose thatthe connection matrix is composed of zero-or-unity elements ( i.e. , T i,j = 0 or ). We assume that,without coupling ( i.e. , C = 0 ), individual systems ( i.e. , ˙ x i = F i ( x i ) ) generate periodic oscillations,after transients. Such closed trajectories in phase space are called limit cycles , which have intrinsicperiods of τ i . The equation (2.2) describes, to a good approximation, a variety of systems ofinterest in biology and neuroscience. Then the theory of phase reduction [58,59] states that, asfar as the coupling strength C is weak enough, the network dynamics can be reduced to thefollowing phase equations: ˙ θ i = ω i + CN N X j =1 ,j = i T i,j Z ( θ i ) G ( θ i , θ j ) = ω i + CN N X j =1 ,j = i T i,j H ( θ j − θ i ) , (2.3)where θ i represents phase of the i th oscillator and ω i gives natural frequency of the i th oscillator( i.e. , ω i = 2 π/τ i ). Z stands for phase sensitivity function (also called “infinitesimal phase responsecurve”), which determines the amount of phase shift induced by the interaction G with otheroscillators (we will here not go in the detail of how equation (2.3) is obtained; an interestedreader can refer to [58–60]). By averaging approximation [59], which integrates one cycle ofthe phase sensitivity function Z with the interaction function G , the coupling function isderived as H ( θ i − θ k ) = π R π Z ( θ i + θ ′ ) G ( θ i + θ ′ , θ k + θ ′ ) dθ ′ . Transformation of the originalequations (2.2) to (2.3) provides a significant reduction in system’s dimensionality, in the sensethat the original state variables x i , which can be high-dimensional, are represented only bythe single phase variable θ i . This substantially simplifies the system’s modeling and enables itsidentification in a straightforward fashion.The individual oscillator states are simultaneously measured as { ξ i ( n∆t ) = η ( x i ( n∆t )) : n =1 , · · · , M } Ni =1 ( η : observation function, ∆t : sampling time, M : data points). This correspondsto an experimental situation, under which states of individual oscillators ( e.g. , gene expressionlevels of individual cells, membrane potentials of neurons, etc. ) are recorded simultaneously. Ourpurpose is to infer the phase equations from these measurement data under the conditions thatthe underlying system equations (2.2) are unknown.The phase dynamics can be reconstructed via the following steps.(i) From the measured data { x i ( t ) } , phases are extracted as θ i ( t ) = 2 πk + 2 π ( t − t k ) / ( t k +1 − t k ) , where t k represents the time, at which i th signal takes its k th peak and t k ≤ t < t k +1 [60]. Note that this method is limited to the case of simple waveform, wherea single peak appears during one oscillation cycle.(ii) Fit the phase equations: ˙ θ i = ˜ ω i + CN N X j =1 ˜ T i,j ˜ H ( θ j − θ i ) (2.4)to the phase data { θ i ( t ) } . Here, { ˜ ω i } represent approximate values for the naturalfrequencies. The coupling function ˜ H , which is in general nonlinear and periodic withrespect to π , is approximated by a Fourier series of pre-selected order D as ˜ H ( ∆θ ) = P Dj =1 a j sin j∆θ + b j (cos j∆θ − . For simplicity, we consider a specific type of coupling, r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A .................................................................. under which the interaction disappears as the phase difference becomes zero, i.e., ˜ H (0) =0 . This type of coupling arises quite often in diffusively coupled oscillator networks[59,61]. [Although general coupling can be also considered, more than one data setsassociated with different coupling strength are required to avoid parameter redundancy.As a simplified demonstration, this study deals with this specific coupling.]The unknown parameters p = { ˜ ω i , a j , b j } are now estimated by the above describedmultiple-shooting method (the connection matrix ˜ T i,j and the coupling strength C = 0 are supposed to be known here).(iii) To avoid over-fitting of the coupling function, cross-validation is utilized to determinethe optimal number of Fourier components D [62]. We divide the measurement datainto two parts. For the first half data, the parameter values p are estimated. Then,the estimated parameters are applied to the latter half data and measure the error e cv = P n | θ (( n + 1) ∆t ) − φ ∆t ( θ ( n∆t ) , p ) | , where φ ∆t ( θ ( n∆t ) , p ) represents ∆t -timefurther state of the phase dynamics (2.4) starting from an initial condition θ ( n∆t ) . Theorder number D providing the minimum error is considered as the optimum. (c) Application to Coupled FitzHugh-Nagumo Oscillators To illustrate how the method described above works, we apply the multiple shooting to aprototypical example of weakly coupled limit cycle oscillators. In the original study [31], coupledRössler oscillators were analyzed. As another yet challenging example, which has more complexshape of coupling function due to the nature of relaxation oscillations, we consider the followingnetwork of FitzHugh-Nagumo (FHN) oscillators [63,64]: ˙ v i = α i ( v i − v i / − w i + I ) + CN N X j =1 T i,j ( x j − x i ) , (2.5) ˙ w i = α i ǫ ( v i + a − by i ) , (2.6)where i = 1 , · · · , N . The system of FitzHugh-Nagumo oscillators can be seen as a simplemodel for interacting neurons. Under the parameter setting of a = 0 . , b = 0 . , ǫ = 0 . , I = 0 . ,individual FHN oscillators (without coupling C = 0 ) gives rise to limit cycles of slow-fast type.Inhomogeneity parameters, which control natural periods of the individual oscillators, were setas α i = 1 + ( i − ∆α ( i = 1 , · · · , N ), where α i = 1 yields natural oscillation period of . .We started with the case of N = 16 . We consider all-to-all coupling matrix ( T i,j = 1 ). The levelof inhomogeneity was set to ∆α = 0 . . The multivariate data { x i ( t ) } i =1 were recorded at acoupling strength of C = 0 . , which is in a non-synchronized regime. The sampling interval wasset to be ∆t = 0 . . Then, the phases { θ i ( t ) } were extracted and down-sampled to a samplinginterval of ∆t = 1000 · . . Total of data points have been collected for the parameterestimation. As an initial condition, the unknown parameter values were all set to be zero, i.e. , ˜ ω i = 0 , a j = b j = 0 . The data were divided into and , which were used for theparameter estimation and the cross-validation error e cv , respectively. By varying the numberof Fourier components from D = 1 to D = 10 , the optimal value was found to be D = 7 by thecross-validation test.The coupling function ˜ H ( ∆θ ) estimated by the present method is in a good agreement with theone computed by the adjoint method [65] (Fig. 1a). The error-bars were computed from inverseof the Hessian matrix of the squared error function, under the assumption that the phase datacontain uncorrelated observational noise [66]. The estimated natural frequencies are distributedon a diagonal line with the true frequencies obtained from simulations of the individual (isolated)FHN oscillators (Fig. 1b). Using the estimated phase equations, synchronization diagram of theoriginal coupled FHN oscillators can be recovered, where the onset of synchronization waspredicted at C = 0 . , which is very close to the real onset of C = 0 . (Fig. 1c).Next, we show how the estimation depends upon the problem setting. The primary factorthat influences the estimation results is the coupling strength C used to generate the time series. r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A .................................................................. Fig. 1 (d) shows dependence of estimation error on the coupling strength. The estimation error e cf is evaluated as deviation of the estimated coupling function ˜ H s ( ∆θ ) from the one ˜ H p ( ∆θ ) estimated by the adjoint method, i.e. , e cf = q R π { ˜ H s ( ∆θ ) − ˜ H p ( ∆θ ) } d∆θ q R π { ˜ H p ( ∆θ ) − h ˜ H p i} d∆θ , (2.7)where the denominator represents normalization factor and h ˜ H p i = (1 / (2 π )) R π ˜ H p ( ∆θ ) d∆θ .As the coupling strength is located close to the onset of synchronization, the estimation errorincreases significantly. Under the synchronized state, phase differences between the oscillators donot change in time ∆θ = const , providing no information about the phase interaction. Increase inestimation error due to synchronized data is therefore reasonable.Even in a synchronized regime, the coupling function can be recovered from transientdata, during which phase differences evolve (transient data often reveals far more informationabout the underlying system, since it is recorded before the system ’settled’ into its dynamicalequilibrium). To show this, the multivariate data were recorded after discarding only a shortduration of transient process that starts from a random initial condition. Transient data (timeinterval of ) were collected before the system reached the final synchronized state. sets ofsuch data were used for the parameter estimation. Fig. 1 (e) shows dependence of the estimationerror on the transient duration. Note that the coupling strength is set to C = 0 . , that is ina synchronized regime. Although the error increases as the transient duration is increased,relatively good estimation was realized for a short transient time. This suggests that, even thesystem is in synchrony with a moderate coupling, application of perturbation that kicks thesystem out of synchrony is an efficient way of inferring the underlying phase dynamics.Fig. 1 (f) shows dependence of the estimation error on the network size N , varied from to . The level of inhomogeneity was set to ∆α = 0 . /N . The multivariate data { y i ( t ) } Ni =1 wererecorded at a coupling strength of C = 0 . , which corresponds to non-synchronized regime.Other settings were the same as those in the case of N = 16 . Surprisingly, the estimation errorremains in a low level. Even for N = 512 , the coupling function ˜ H ( ∆θ ) has been preciselyestimated, while the estimated natural frequencies { ω i } i =1 are consistent with those obtainedfrom the non-coupled simulations. This suggests that the system size does not pose a majorlimit on the estimation of phase dynamics as far as the data contain non-synchronized phaseinformation.Although the coupling function has been reliably estimated for networks with all-to-allconnections ( T i,j = 1 ), the estimation error may increase when oscillators are heterogeneouslyconnected to each other. We deal with such situations in the next section.
3. Application to Network Inference
Although we have dealt with the case that all oscillator elements are connected to all the others inthe previous section, heterogeneous connections are more common in nature and engineering. Asanother challenge of our technique [53], this section discusses a problem of inferring connectivityof the oscillator network from the measured time series. Numerous approaches have beenproposed up to date using information transfer [67], mutual predictability [68], recurrenceproperties [69], permutation-based asymmetric association measure [70], index for partial phasesynchronization [71–73], and graph theory [74]. Response properties of the network dynamicsto external stimuli have been also exploited [8,75]. For weakly coupled limit cycle oscillators, towhich phase reduction is applicable, the phase modeling approach is again quite effective fordetecting the network topology [11,53,76–78].In our approach [53], the multiple-shooting method is again applied to fit the phaseequations (2.4) to the phase data { θ i ( t ) } . The fitting procedure is the same as before except thatthe connection matrix is estimated as the unknown parameters p = { ˜ T i,j } . For simplicity, the r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A .................................................................. coupling function ˜ H and the natural frequencies { ω i } Ni =1 of the oscillator elements were assumedto be known (the general case that both coupling function and connection matrix are unknownhas been dealt with in the previous study [53]). As coefficients { a j , b j } for the coupling function,the ones estimated in the previous section were utilized. Natural frequencies { ω i } Ni =1 were alsoobtained from the simulations of non-coupled original equations.As the target system, the network of FHN oscillators (2.5),(2.6) were studied. For a network offour ( N = 4 ) oscillators, two defects were introduced to the connection matrix as T , = T , = 0 (here, defect means that one oscillator is not connected to another). The level of inhomogeneitywas set to ∆α = 0 . , whereas the coupling strength was C = 0 . , i.e. , in a non-synchronizedregime. As in the previous section, total of data points (sampling time: ) have been collected.By the multiple-shooting method, the connection matrix was estimated as follows. ˜ T , ˜ T , ˜ T , ˜ T , ˜ T , ˜ T , ˜ T , ˜ T , ˜ T , ˜ T , ˜ T , ˜ T , = . ± .
01 1 . ± .
01 1 . ± . . ± .
02 1 . ± .
02 0 . ± . − . ± .
01 1 . ± .
01 1 . ± . . ± .
01 0 . ± .
01 1 . ± . . We see that the two defects ( ˜ T , , ˜ T , ) were precisely identified as small values, whereas othermatrix elements pointed almost unity.For comprehensive analysis, the connection matrices with randomly generated defects wereestimated for variable network size from N = 2 to N = 16 . For our analysis, the estimation errorwas evaluated as e cm = N ( N − P Ni =1 P Nj =1 ,j = i | ˘ T i,j − T i,j | , where the estimated connectivitywas digitized as ˘ T i,j = 0 for ˜ T i,j < . and ˘ T i,j = 1 otherwise. For each setting of the networksize, instances of connection matrices { T i,j } were randomly generated and the average and thestandard deviation of the estimation errors were plotted in Fig. 2. Panels (a) and (b) correspondto the cases that defect ratios ( i.e. , percentage of zero elements in the connection matrix) are % and %, respectively. For variable network sizes, the estimation errors e cm are almost zeroexcept N = 11 , in the case of low defect ratio. Although the errors increase for high defect ratio,their overall level is less than 0.2.To examine the effect of coupling function, the connection matrices were also estimatedby using a simple sine as the coupling function, i.e. , ˜ H ( ∆θ ) = a sin ∆θ . For small networks,the difference was not large between the precisely estimated (higher-order) coupling function(red solid line) and the simple sine function (blue dotted line). However, as the network sizeincreases, the estimation error increases much more rapidly in the sine function than in the higher-order coupling function. This indicates that, for reliable detection of the connectivities, preciselyestimated coupling function is of significant importance.In Fig. 2 (c), dependence of the network inference on data length M is indicated. For networksizes of N = 6 and N = 8 , we have varied the data length and studied how it affected theestimation results of the network connectivity. The defect ratio was set to %. The networkinference was reliable for data length longer than 200 points ( i.e. , about 20 cycles). For shorterdata length, the estimation error gradually increased. It is therefore crucial to utilize enough datalength for precisely detecting the network connectivity.Fig. 2 (d) shows dependence of the network inference on Gaussian noise N (0 , (2 πσ ) ) addedto the phase data. The defect ratio and the data length were set to % and M = 400 , respectively.The estimation error increased gradually as the noise level was increased, where σ = 0 . % and σ = 2 % of phase noises caused severe damages on the network inference for system sizes of N = 8 and N = 6 , respectively. This suggests that our estimation technique is rather sensitive tothe phase noise and, for reliable estimation of the connection matrix, phase information shouldbe accurately extracted from the observed time series.
4. Further Applications
In this section, we discuss further applications of the multiple-shooting technique. r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A .................................................................. (a) Inferring Phase Sensitivity Function First, we apply the multiple-shooting method to the estimation of phase sensitivity function Z ( θ ) .The phase sensitivity function Z ( θ ) plays a vital role in the studies of coupled oscillators, since itdescribes one of the most fundamental properties of the oscillator element [58–60]. Numerousapproaches have been proposed to estimate the phase sensitivity function from experimentaldata [43–52]. As an extension of our technique, the phase sensitivity function can be recoveredfrom the coupling function [79]. As described earlier in the averaging approximation [59], thecoupling function is given by a convolution of the phase sensitivity function Z ( θ ) and theinput waveform G ( θ ) as H ( θ ) = π R π Z ( ψ ) G ( θ + ψ ) dψ = ( Z ∗ G )( θ ) . It is straightforward torecover the phase sensitivity function by the spectral deconvolution [80]. Namely, in a frequencydomain, the phase sensitivity function is given as ˆ Z ( ω ) = ˆ H ( ω ) / ˆ G ( ω ) , where ˆ Z ( ω ) , ˆ H ( ω ) , and ˆ G ( ω ) represent Fourier transforms of Z , H , and G , respectively. Inverse Fourier transform of ˆ Z ( ω ) yields the phase sensitivity Z ( θ ) . Fig. 3(a) shows phase sensitivity function (solid red line)obtained by the deconvolution of the coupling function estimated from coupled FHN oscillators( N = 16 ) in section 2. Compared with the one computed by the adjoint method [65] (dashed blueline), the estimated function is somewhat deviated from the true curve. We consider that, dueto the averaging effect, where the effect of input signal is averaged over one oscillation cycle,information on the spontaneous phase response has been lost.To improve the situation, the phase sensitivity can be estimated more directly by using theWinfree formula [58] as follows. For simplicity, we consider a single phase oscillator receiving l -th external force G l ( t ) ( l = 1 , , . . ., L ): ˙ θ l = ω + ˜ Z ( θ l ) G l ( t ) , (4.1)where θ l and ω represent phase and natural frequency of the oscillator. Without loss of generality,the initial phase can be set to zero ( i.e. , θ l (0) = 0 ). The external force G l ( t ) is typically composed ofa short pulse, which lasts within one oscillator cycle of T = 2 π/ω . The phase sensitivity function ˜ Z is described in terms of a Fourier series as ˜ Z ( θ ) = c + P Dj =1 c j sin jθ + d j cos jθ . The unknowncoefficients p = { c j , d j } can be estimated by the multiple-shooting method in a similar manner asthe estimation of coupling function. Given the external force G l ( t ) , the phase oscillator model (4.1)can be integrated as φ T ( θ l (0) , G l , p ) . The parameters p can be optimized in such a way thatthe phase model (4.1) satisfies the boundary conditions: θ l ( T ) = φ T ( θ l (0) , G l , p ) , where θ l ( T ) represents the oscillator phase observed at t = T .Below, we compare the performance of multiple-shooting method with that of least squaresas the standard method of estimating the phase sensitivity function [43,46]. Here, the phasemodel (4.1) is integrated as Z T dθ l = Z T ωdt + Z T ˜ Z ( θ l ) G l ( t ) dt,θ l ( T ) − θ l (0) − π = Z T ˜ Z ( ωt ) G l ( t ) dt, where the oscillator phase is approximated as θ l ( t ) ≈ ωt under the assumption that the externalforce G l ( t ) is weak in equation (4.1). By expanding the external force into Fourier series as G l ( t )= g l, + P Dj =1 g l,j sin jωt + h l,j cos jωt , we obtain M p = D , r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A .................................................................. where M = g , / g , h , g , h , · · · g ,D h ,D g , / g , h , g , h , · · · g ,D h ,D ... ... ... ... ... ... ... g L, / g L, h L, g L, h L, · · · g L,D h L,D , p = c c d ... c D d D , D = ∆θ ∆θ ... ∆θ L .∆θ l = θ l ( T ) − θ l (0) represents phase shift induced by the l -th external force G l ( t ) . The unknowncoefficients p can be obtained as p = DM − .We apply the two methods to a single FHN oscillator that receives random impulses(stimulus duration: τ = 20 , stimulus strength: V = 0 . , . , .., . ) as external forcing G ( t ) .Parameter values of the FHN oscillator and the sampling time interval were set to be the sameas those in the previous sections. For simplicity, natural frequency ω and the external signal G ( t ) were assumed to be known. Number of the Fourier components was set to D = 10 . Theintegration time was set to T = 150 . For impulse strength of E = 0 . and E = 0 . , the estimatedphase sensitivity functions are drawn in Fig. 3(b) and (c), respectively. In both panels (b) and (c),estimation results of the multiple shooting method (solid red lines) are consistent with those ofthe adjoint method [65]. The least-square method (dashed blue line), on the other hand, recoveredthe phase sensitivity function faithfully for a small impulse strength in (b). The estimate is,however, deviated from the other two curves for a large impulse strength in (c). In fact, as theimpulse strength is increased, the estimation error increases much more rapidly in the least-square method (dashed blue line) than the multiple shooting method (solid red line) (see Fig. 3d).The least-square method [43,46] assumes that phase of the oscillator evolves linearly in timeaccording to the natural frequency. This approximation is effective as far as the external force isweak. If stronger perturbations are applied, inducing non-small phase shifts, this approximationincreases the estimation error. The multiple-shooting method, on the other hand, takes intoaccount the phase shift induced by the external perturbations by faithfully integrating the phaseequation (4.1). The estimation error has been therefore reduced by the multiple-shooting method. (b) Chaotic Phase Synchronization Next we show how the estimated coupling function can be utilized for modeling chaotic phasesynchronization [81]. It has been known that phases of chaotic oscillators can be synchronizedwith each other, while their amplitudes remain irregular and uncorrelated. Especially for phase-coherent chaos, in which rotation center can be well-defined, the phase dynamics can beapproximated as ˙ θ = ω + Γ ( A ) , where Γ ( A ) represents frequency modulation, which dependsupon oscillation amplitude A [81]. For chaotic amplitude A , the term Γ ( A ) can be regarded as aneffective noise. In many phase-coherent systems such as the Rössler equations [82], amplitude-dependent frequency modulation is very small, so the noise term Γ ( A ) is negligible. Phasedynamics of such chaotic attractor becomes very similar to those of limit cycle oscillators.To extract phase-interaction between chaotic oscillators, we consider two coupled Rösslerequations [82]: ˙ x , = − α , y , − z , , ˙ y , = α , x , + 0 . y , + C ( y , − y , ) , ˙ z , = 0 . z , ( x , − . Each Rössler oscillator gives rise to chaotic dynamics without coupling C = 0 . The inhomogeneityparameters were set as α , = 1 ∓ . , which yield average oscillation periods of . and . ,respectively. The bivariate data { y i ( t ) } i =1 were simulated under coupling strength of C = 0 . ,which corresponds to non-synchronized regime. The sampling interval was set to be ∆t = 0 . for the extraction of the phases { θ i ( t ) } . Then, to apply the multiple-shooting method, the data r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A .................................................................. have been down sampled to ∆t = 1000 · . and the total of data points were collected. Thedata were divided into and points, which were used for the parameter estimation andthe cross-validation test, respectively. By varying the number of Fourier components from D = 1 to D = 5 , the optimal value was found to be D = 4 . The corresponding coupling function ˜ H ( ∆θ ) is shown by the solid red line in Fig. 4(c). The estimated function is in good agreement with theone obtained by the convolution of averaged phase sensitivity function (Fig. 4a) and the averagedinput waveform (Fig. 4b). Using the estimated phase equations, synchronization diagram of theoriginal two coupled Rössler equations can be recovered, where the onset of synchronization waspredicted at C = 0 . , which is close to the real onset of C = 0 . (Fig. 4d). This suggests that oursimple method of estimating the coupling function provides a good approximate of describingthe phase dynamics of phase-coherent chaotic oscillators. (c) Application to Circuit Experiment Finally, we apply our method to experimental data generated from Van der Pol electric circuit [83]to demonstrate the performance of our method in a realistic experimental setting. As shown inthe diagram of Fig. 5a, the system is based on a LC circuit, composed of an inductor ( L ) anda capacitor ( C ). To form a negative-resistance converter, three positive resistors ( R , R , R )were connected to a voltage-controlled voltage source ( i.e. , operational amplifier and its associatedpower supplies V DD , V SS ) [84]. External forcing G ( t ) was injected from a function generator(Keysight 33500B) to the Van der Pol circuit through a capacitor ( C ). Physical parameters of theelectric components used in the present experiment are summarized in Table 1. To obtain thephase sensitivity function, impulses (stimulus duration: τ = 380 µ s, stimulus strength: V = 3 V) were randomly injected as the external force G ( t ) . The circuit output as well as the inputimpulses were simultaneously measured with a sampling frequency of . kHz. First, the phasesensitivity function ˜ Z was estimated by fitting the phase model (4.1) to the measured data with themultiple-shooting method. Natural frequency f n = 110 . Hz ( i.e. , ω = 2 πf n ), measured before thestimulus experiment, was used in the phase dynamics. Number of the Fourier components wasset to D = 4 . As shown in Fig. 5b, the estimated phase sensitivity ˜ Z ( θ ) fits to the experimentalobservation of phase response data well.Next, a sinusoidal forcing G ( t ) = V sin( Ωt ) (forcing frequency: Hz, forcing amplitude: V =0 . V) was applied to the Van der Pol circuit. The circuit output as well as the forcing waveformswere simultaneously measured with a sampling frequency of . kHz. By the multiple-shootingmethod, which fits the phase equations (2.4) to the measured data, the coupling function ˜ H (number of Fourier components: D = 1 ) was estimated. In Fig. 5c, the estimated coupling functionis compared with the one obtained by the averaging of the phase sensitivity function ˜ Z , estimatedfrom the impulse stimuli, and the input sine waveform G ( t ) . Despite a slight difference in theinitial phase, the coupling functions agree quite well with each other. In Fig. 5d, the estimatedphase equations recovered synchronization diagram of the experimental system, where the onsetof synchronization was predicted at V = 0 . V, which is very close to the real onset of V = 7 V.
5. Discussions and Conclusions
The multiple-shooting method has been focused on as a non-invasive approach to estimatecoupling functions from multivariate time series measured from a real or synthetic complexdynamical system [31]. Among various methods developed so far [7,32–36,38–41], whichare based on the Bayesian estimation and other variants, the multiple-shooting provides astraightforward approach to fit the phase equations to phase data measured from an oscillatornetwork. Despite its simplicity, the method was shown to be capable of precisely estimatingthe coupling function of the coupled FHN oscillators including higher-harmonic terms. Theestimation was found effective for a large network of up to oscillators. Utilization of thetransient part of data successfully enlarged applicability of the estimation technique even ina synchronized regime of coupled oscillators. The estimated coupling function was further r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A .................................................................. Table 1.
Parameters of Van der Pol circuit. L [mH] C . [uF] R . [K Ω ] R . [k Ω ] R [k Ω ] V DD V SS -5 [V]OPAMP LF412CN C [nF]applied to inference of network topology and chaotic phase synchrony. Precise estimation of thecoupling functions was shown to improve the reconstruction of network topology. As anotherintriguing issue, estimation of the phase sensitivity function was also discussed. Although thephase sensitivity function obtained by deconvolution of the estimated coupling function wasslightly deviated from the true function, refinement has been made by extending the multipleshooting method directly to the phase data of a driven limit cycle oscillator. Finally, efficiency ofthe present approach was demonstrated with the experimental data measured from the Van derPol electric circuit with a sinusoidal forcing.Beyond experimental system in physics, chemistry, and engineering, we foresee that ourmethod will be applicable to system of rhythmic, interacting elements such as cellular geneexpressions in the suprachiasmatic nucleus (SCN) [56], electrical activities of cardiac pacemakers[57], inferior olive neurons in the cerebellum [85] and can give insights useful for domain-scientists in biology and neuroscience.While considering our method of potentially practical use for various systems, its usefulnessis not without limitations. The main among them is the assumption that the studied system canbe approximated as a network of weakly coupled limit cycles [59]. This, however, is not truefor all systems encountered in nature. For instance, in gene regulatory networks, phases of theclock component genes are tightly connected to each other [86]. It has been known that corticalneurons fire with a strong synchrony during epileptic seizure [87]. Such strongly coupled systemsshould be carefully distinguished and avoided as a target of modeling the phase dynamics. Inthe case that the system property is not well understood, it is nontrivial to judge only from therecorded data whether the coupling is weak enough to apply the phase modeling to the oscillatornetwork. It is an important open problem to provide a criterion to assess whether the phase modelis suitable for analyzing the observed time series without prior knowledge on the underlyingdynamical equations.Another limitation is the length of the available time series: namely, experimentalmeasurements, for a variety of realistic reasons, could produce the data (time series) of only a veryshort length. For instance, time resolved data on gene regulation are not likely to yield time serieswith much more than 10 cycles. In this case, our method might be of limited use. Also, realisticdata are almost always noisy. The noise strength, depending on the experimental scenario, couldbe quite severe. Especially, the phase extraction process in our modeling is rather sensitiveto noise. Temporal fluctuation and noise in natural frequencies of the oscillator elements mayalso cause estimation error in the coupling functions. In this respect, noise tolerance should becarefully examined, before the application to data contaminated with observational/dynamicalnoise.Also, networks in real world are large and only partials of the dynamics elements areobservable. Although our method was shown to be robust against system size as far asthe oscillator elements are uniformly connected and they are desynchronized, the effect of r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A .................................................................. unobserved oscillator states should be examined carefully. Heterogeneity and hierarchy in thecoupling functions may require further extension of the present approach.Finally, we conclude the paper with a brief discussion of how our method’s performancecompares to the performance of other methods that reconstruct coupling functions in oscillatorysystems. Unfortunately, such comparison is not simple to make, since various available methodsdepart from very different hypotheses and knowledge about the system. Stronger hypotheseslead to better inferences, but the information on whether the hypotheses are met is not alwaysavailable. This renders hard any independent comparison of reconstruction methods. One couldargue that methods aimed at only network topology are more useful and precise, but suchmethods neglect the entire dynamical nature of many real networks. On the other hand, certainmethods give excellent results, but are limited to dynamical systems with specific properties. Infact, our method belong to this category, since it assumes the limit cycle nature of the individualunits. Furthermore, methods can be divided into invasive ones (that interfere with system’songoing dynamics) and non-invasive ones (that do not). Again, their real merits is hard tocompare, since invasive methods, although often non practical, will almost always give betterresults. Therefore, we here conclude that our reconstruction concept, although limited by theassumption of limit cycles, is a promising – and above all practical – approach implementable inreal experiments.Data Accessibility. Experimental data generated from electric circuit are available in Dryad dataset(https://doi.org/10.5061/dryad.z34tmpg80).
Authors’ Contributions.
ITT designed the study. ITT performed the numerical simulations and the dataanalysis. KI carried out the circuit experiments. ITT and ZL wrote the manuscript. All authors read andapproved the manuscript.
Competing Interests.
The author(s) declare that they have no competing interests.
Funding.
ITT acknowledges support by Grant-in-Aid for Scientific Research (No. 17H06313, No. 16H04848,No. 16K00343, No. 18H02477) from Japan Society for the Promotion of Science (JSPS). ZL acknowledgessupport by EU via H2020-MSCA-ITN-2015 COSMOS 642563, and by Slovenian research agency via P1-0383and J5-8236.
References
1. Newman M. 2007
Networks, An Introduction .Oxford University Press.2. Estrada E. 2011
The Structure of Complex Networks: Theory and Applications .Oxford University Press.3. Porter M, Gleeson J. 2016
Dynamical Systems on Networks .Springer.4. Barabasi AL. 2016
Network Science .Cambridge University Press.5. Easley D, Kleinberg J. 2012
Networks, Crowds, and Markets: Reasoning About a Highly ConnectedWorld .Cambridge University Press.6. Boccaletti S, Bianconi G, Criado R, Del Genio CI, Gómez-Gardenes J, Romance M, Sendina-Nadal I, Wang Z, Zanin M. 2014 The structure and dynamics of multilayer networks.
Physics Reports , 1–122.7. Blaha KA, Pikovsky A, Rosenblum M, Clark MT, Rusin CG, Hudson JL. 2011 Reconstructionof two-dimensional phase dynamics from experiments on coupled oscillators.
Physical Review E , 046201.8. Levnaji´c Z, Pikovsky A. 2011 Network reconstruction from random phase resetting. Physical review letters , 034101.9. Wang WX, Yang R, Lai YC, Kovanis V, Harrison MAF. 2011 Time-series–based prediction ofcomplex oscillator networks via compressive sensing.
EPL (Europhysics Letters) , 48006.10. Stankovski T, Duggento A, McClintock PV, Stefanovska A. 2012 Inference of time-evolvingcoupled dynamical systems in the presence of noise. r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A .................................................................. Physical review letters , 024101.11. Kralemann B, Pikovsky A, Rosenblum M. 2014 Reconstructing effective phase connectivity ofoscillator networks from observations.
New Journal of Physics , 085013.12. Levnaji´c Z, Pikovsky A. 2014 Untangling complex dynamical systems via derivative-variablecorrelations. Scientific reports , 5030.13. Timme M, Casadiego J. 2014 Revealing networks from dynamics: an introduction. Journal of Physics A: Mathematical and Theoretical , 343001.14. Han X, Shen Z, Wang WX, Di Z. 2015 Robust reconstruction of complex networks from sparsedata. Physical review letters , 028701.15. Nitzan M, Casadiego J, Timme M. 2017.Revealing physical network interactions from statistics of collective dynamics sci.16. Wang WX, Lai YC, Grebogi C. 2016 Data based identification and prediction of nonlinear andcomplex dynamical systems.
Physics Reports , 1–76.17. Leguia MG, Andrzejak RG, Levnaji´c Z. 2017 Evolutionary optimization of networkreconstruction from derivative-variable correlations.
Journal of Physics A: Mathematical and Theoretical , 334001.18. Rosenblum M, et al. Physical Review E , 012209.19. Simidjievski N, Tanevski J, Ženko B, Levnaji´c Z, Todorovski L, Džeroski S. 2018 Decouplingapproximation robustly reconstructs directed dynamical networks. New Journal of Physics , 113003.20. Papana A, Kyrtsou C, Kugiumtzis D, Diks C. 2013 Simulation study of direct causalitymeasures in multivariate time series. Entropy , 2635–2661.21. Koutlis C, Kugiumtzis D. 2016 Discrimination of coupling structures using causality networksfrom multivariate time series. Chaos: An Interdisciplinary Journal of Nonlinear Science , 093120.22. Zanin M, Papo D, Sousa PA, Menasalvas E, Nicchi A, Kubik E, Boccaletti S. 2016 Combiningcomplex networks and data mining: why and how. Physics Reports , 1–44.23. Bridewell W, Langley P, Todorovski L, Džeroski S. 2008 Inductive process modeling.
Mach. Learn. , 1–32.24. Džeroski S, Todorovski L. 2008 Equation discovery for systems biology: finding the structureand dynamics of biological networks from time course data. Current Opinion in Biotechnology , 360–368.25. Schmidt M, Lipson H. 2009 Distilling free-form natural laws from experimental data. science , 81–85.26. Brunton SL, Proctor JL, Kutz JN. 2016 Discovering governing equations from data by sparseidentification of nonlinear dynamical systems. Proceedings of the National Academy of Sciences , 3932–3937.27. Simidjievski N, Todorovski L, Džeroski S. 2016 Modeling dynamic systems with efficientensembles of process-based models.
PLoS One , e0153507.28. Tanevski J, Todorovski L, Džeroski S. 2016 Learning stochastic process-based models ofdynamical systems from knowledge and data. BMC systems biology , 30.29. Tanevski J, Todorovski L, Džeroski S. 2016 Process-based design of dynamical biologicalsystems. Scientific reports , 34107.30. ˇCerepnalkoski D, Taškova K, Todorovski L, Atanasova N, Džeroski S. 2012 The influence ofparameter fitting methods on model structure selection in automated modeling of aquaticecosystems. Ecological Modelling , 136–165. r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A ..................................................................
31. Tokuda IT, Jain S, Kiss IZ, Hudson JL. 2007 Inferring phase equations from multivariate timeseries.
Physical review letters , 064101.32. Kralemann B, Cimponeriu L, Rosenblum M, Pikovsky A, Mrowka R. 2007 Uncoveringinteraction of coupled oscillators from data. Physical Review E , 055201.33. Kralemann B, Cimponeriu L, Rosenblum M, Pikovsky A, Mrowka R. 2008 Phase dynamics ofcoupled oscillators reconstructed from data. Physical Review E , 066205.34. Cadieu CF, Koepsell K. 2010 Phase coupling estimation from multivariate phase statistics. Neural computation , 3107–3126.35. Kralemann B, Pikovsky A, Rosenblum M. 2011 Reconstructing phase dynamics of oscillatornetworks. Chaos: An Interdisciplinary Journal of Nonlinear Science , 025104.36. Zhu Y, Hsieh YH, Dhingra RR, Dick TE, Jacono FJ, Galán RF. 2013 Quantifying interactionsbetween real oscillators with information theory and phase models: application tocardiorespiratory coupling. Physical Review E , 022709.37. Pikovsky A. 2018 Reconstruction of a random phase dynamics network from observations. Physics Letters A , 147–152.38. Suzuki K, Aoyagi T, Kitano K. 2018 Bayesian estimation of phase dynamics based on partiallysampled spikes generated by realistic model neurons.
Frontiers in computational neuroscience , 116.39. Miyazaki J, Kinoshita S. 2006 Determination of a coupling function in multicoupledoscillators. Physical review letters , 194101.40. Tokuda IT, Wagemakers A, Sanjuán MA. 2010 Predicting the synchronization of a network ofelectronic repressilators. International Journal of Bifurcation and Chaos , 1751–1760.41. Stankovski T, Pereira T, McClintock PV, Stefanovska A. 2017 Coupling functions: universalinsights into dynamical interaction mechanisms. Reviews of Modern Physics , 045001.42. Kiss IZ, Zhai Y, Hudson JL. 2005 Predicting mutual entrainment of oscillators withexperiment-based phase models. Physical review letters , 248301.43. Galán RF, Ermentrout GB, Urban NN. 2005 Efficient estimation of phase-resetting curves inreal neurons and its significance for neural-network modeling. Physical review letters , 158101.44. Ota K, Omori T, Aonishi T. 2009 Map estimation algorithm for phase response curves basedon analysis of the observation process. Journal of computational neuroscience , 185.45. Nakae K, Iba Y, Tsubo Y, Fukai T, Aoyagi T. 2010 Bayesian estimation of phase response curves. Neural Networks , 752–763.46. Hong S, Robberechts Q, Schutter ED. 2012 Efficient estimation of phase-response curves viacompressive sensing. Journal of neurophysiology , 2069–2081.47. Goldberg JA, Atherton JF, Surmeier DJ. 2013 Spectral reconstruction of phase response curvesreveals the synchronization properties of mouse globus pallidus neurons.
Journal of neurophysiology , 2497–2506.48. Saifee TA, Edwards MJ, Kassavetis P, Gilbertson T. 2015 Estimation of the phase responsecurve from parkinsonian tremor.
Journal of neurophysiology , 310–323.49. Imai T, Ota K, Aoyagi T. 2017 Robust measurements of phase response curves realized viamulticycle weighted spike-triggered averages.
Journal of the Physical Society of Japan , 024009.50. Cestnik R, Rosenblum M. 2018 Inferring the phase response curve from observation of acontinuously perturbed oscillator. Scientific Reports , 13606. r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A ..................................................................
51. Okada M, Kaburagi T, Tokuda IT. 2019 Acoustic measurements of the infinitesimal phaseresponse curve from a sounding organ pipe.
Physics Letters A , 1733–1741.52. Nakae K. 2019 Statistical estimation of phase response curves using data transformation.
Journal of the Physical Society of Japan , 084003.53. Tokuda IT, Wickramasinghe M, Kiss IZ. 2013 Detecting connectivity of small, dense oscillatornetworks from dynamical measurements based on a phase modeling approach. Physics Letters A , 1862–1867.54. Baake E, Baake M, Bock H, Briggs K. 1992 Fitting ordinary differential equations to chaoticdata.
Physical Review A , 5524.55. Press WH, Teukolsky SA, Vetterling WT, Flannery BP. 2007 Numerical recipes 3rd edition: Theart of scientific computing .Cambridge university press.56. Yamaguchi S, Isejima H, Matsuo T, Okura R, Yagita K, Kobayashi M, Okamura Hs. 2003Synchronization of Cellular Clocks in the Suprachiasmatic Nucleus.
Science , 1408–1412.57. Verheijck EE, Wilders R, Joyner RW, Golod DA, Kumar R, Jongsma HJ, Bouman LN, vanGinneken AC. 1998 Pacemaker synchronization of electrically coupled rabbit sinoatrial nodecells.
The Journal of general physiology , 95–112.58. Winfree AT. 2001
The geometry of biological time .Springer, New York.59. Kuramoto Y. 1984
Chemical oscillations, waves, and turbulence .Springer, Berlin.60. Pikovsky A, Rosenblum M, Kurths J. 2003
Synchronization: a universal concept in nonlinearsciences , volume 12.Cambridge university press.61. Acebrón JA, Bonilla LL, Vicente CJP, Ritort F, Spigler R. 2005 The kuramoto model: A simpleparadigm for synchronization phenomena.
Reviews of modern physics , 137.62. Stone M. 1977 An asymptotic equivalence of choice of model by cross-validation and akaike’scriterion. Journal of the Royal Statistical Society. Series B (Methodological) pp. 44–47.63. FitzHugh R. 1961 Impulses and physiological states in theoretical models of nerve membrane.
Biophysical journal , 445–466.64. Nagumo J, Arimoto S, Yoshizawa S. 1962 An active pulse transmission line simulating nerveaxon. Proceedings of the IRE , 2061–2070.65. Ermentrout B. 1996 Type i membranes, phase resetting curves, and synchrony. Neural computation , 979–1001.66. Horbelt W, Timmer J, Bünner M, Meucci R, Ciofini M. 2001 Identifying physical properties ofa co 2 laser by dynamical modeling of measured time series. Physical Review E , 016222.67. Schreiber T. 2000 Measuring information transfer. Physical review letters , 461.68. Rosenblum MG, Pikovsky AS. 2001 Detecting direction of coupling in interacting oscillators. Physical Review E , 045202.69. Romano MC, Thiel M, Kurths J, Grebogi C. 2007 Estimation of the direction of the couplingby conditional probabilities of recurrence. Physical Review E , 036211.70. Hempel S, Koseska A, Kurths J, Nikoloski Z. 2011 Inner composition alignment for inferringdirected networks from short time series. Physical review letters , 054101.71. Schelter B, Winterhalder M, Dahlhaus R, Kurths J, Timmer J. 2006 Partial phasesynchronization for multivariate synchronizing systems.
Physical review letters , 208103.72. Nawrath J, Romano MC, Thiel M, Kiss IZ, Wickramasinghe M, Timmer J, Kurths J, Schelter r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A .................................................................. B. 2010 Distinguishing direct from indirect interactions in oscillatory networks with multipletime scales.
Physical review letters , 038701.73. Wickramasinghe M, Kiss IZ. 2011 Phase synchronization of three locally coupled chaoticelectrochemical oscillators: Enhanced phase diffusion and identification of indirect coupling.
Physical Review E , 016210.74. Runge J, Heitzig J, Petoukhov V, Kurths J. 2012 Escaping the curse of dimensionality inestimating multivariate transfer entropy. Physical review letters , 258701.75. Timme M. 2007 Revealing network connectivity from response dynamics.
Physical review letters , 224101.76. Stankovski T, Ticcinelli V, McClintock PV, Stefanovska A. 2015 Coupling functions in networksof oscillators. New Journal of Physics , 035002.77. Stankovski T, Petkoski S, Raeder J, Smith AF, McClintock PV, Stefanovska A. 2016 Alterationsin the coupling functions between cortical and cardio-respiratory oscillations due toanaesthesia with propofol and sevoflurane. Phil. Trans. R. Soc. A , 20150186.78. Stankovski T, Ticcinelli V, McClintock PV, Stefanovska A. 2017 Neural cross-frequencycoupling functions.
Frontiers in systems neuroscience , 33.79. Kralemann B, Pikovsky A, Rosenblum M. 2013 Detecting triplet locking by tripletsynchronization indices. Physical Review E , 052904.80. Wiener N. 1949 Extrapolation, interpolation, and smoothing of stationary time series: withengineering applications .MIT press Cambridge, MA.81. Rosenblum MG, Pikovsky AS, Kurths J. 1996 Phase synchronization of chaotic oscillators.
Physical review letters , 1804.82. Rössler OE. 1976 An equation for continuous chaos. Physics Letters A , 397–398.83. Van der Pol B, Van Der Mark J. 1927 Frequency demultiplication. Nature , 363.84. Kennedy MP. 1992 Robust op amp realization of chua’s circuit.
Frequenz , 66–80.85. Tokuda IT, Hoang H, Kawato M. 2017 New insights into olivo-cerebellar circuits for learningfrom a small training sample. Current opinion in neurobiology , 58–67.86. Zhang R, Lahens NF, Ballance HI, Hughes ME, Hogenesch JB. 2014 A circadian geneexpression atlas in mammals: implications for biology and medicine. Proceedings of the National Academy of Sciences , 16219–16224.87. Velazquez JLP, Carlen PL. 2000 Gap junctions, synchrony and seizures.
Trends in neurosciences , 68–74. r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A .................................................................. (a) -0.5-0.25 0 0.25 0.5 0.75 0 π I n t e r ac ti on F un c ti on Phase Difference [rad]Adjoint methodMS method (b) E s ti m a t e d p e r i od s o f o s c ill a t o r s Natural periods of oscillators (c) S t a nd a r d d e v i a ti on o f fr e qu e n c i e s Coupling Strength Real simulationEstimated model (d) E s ti m a ti on e rr o r Coupling strength (e) -1 0 1 2 3 4 5 200 300 400 500 600 700 800 900 E s ti m a ti on e rr o r Transient time (f) E s ti m a ti on e rr o r Number of oscillators
Figure 1.
Results for a network of N = 16 FHN oscillators. (a) Coupling functions ˜ H ( ∆θ ) estimated by the presentmethod (solid red line) and the adjoint method (dashed blue line). (b) Estimated natural frequencies (ordinate) { ω i } i =1 ofFHN oscillators vs. those obtained from non-coupled simulation (abscissa). (c) Synchronization diagrams of the estimatedmodel (solid red line) and the original coupled oscillators (dashed blue line). (d) Dependence of estimation error on thecoupling strength C used to generate multivariate data. The estimation error e is defined as the deviation of the estimatedcoupling function from the one computed by the adjoint method. (e) Dependence of the estimation error on the transienttime, after which the multivariate data were sampled. The coupling strength was set to C = 0 . . (f) Dependence of theestimation error on the number of oscillators N . r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A .................................................................. (a) -0.1 0 0.1 0.2 0.3 0.4 0.5 2 4 6 8 10 12 14 16 E s ti m a ti on e rr o r Number of oscillatorsSine functionHigher-order coupling function (b) -0.1 0 0.1 0.2 0.3 0.4 0.5 2 4 6 8 10 12 14 16 E s ti m a ti on e rr o r Number of oscillatorsSine functionHigher-order coupling function (c) E s ti m a ti on e rr o r Data points N=8N=6 (d) E s ti m a ti on e rr o r Noise level N=8N=6
Figure 2.
Estimation errors of the network connectivity. (a) Percentage of non-connected pairs of oscillators is %. Thecoupling function is composed of higher-order ( D = 5 ) Fourier components in solid red line, while it is based on a simplesine function in dashed blue line. (b) Percentage of non-connected pairs of oscillators is %. (c) Dependence of theestimation error on data length. Percentage of non-connected pairs of oscillators is %, while number of the oscillatorsis set to N = 6 (solid red line) and N = 8 (dashed blue line). (d) Dependence of the estimation error on noise level σ ,where Gaussian noise N (0 , (2 πσ ) ) is added to the phase data. Percentage of the non-connected pairs of oscillators is %, while number of the oscillators is set to N = 6 (solid red line) and N = 8 (dashed blue line). r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A .................................................................. (a) -1.5-1-0.5 0 0.5 1 1.5 0 π 2π P h a s e S h i f t [r a d ] Phase Difference [rad]Deconvolution of coupling functionAdjoint method (b) -1.5-1-0.5 0 0.5 1 1.5 0 π 2π P h a s e s h i f t Phase DifferenceLeast-square methodMS methodAdjoint method (c) -1.5-1-0.5 0 0.5 1 1.5 0 π 2π P h a s e s h i f t [r a d ] Phase [rad]Least-square methodMS methodAdjoint method (d) E s ti m a ti on e rr o r Strength of ImpulseLeast-square methodMS method
Figure 3. (a) Phase sensitivity function Z (solid red line) obtained by deconvolution of the coupling function estimated inFig. 1a. Compared is the estimate by the adjoint method (dotted black line). (b,c) Phase sensitivity functions Z obtainedby MS method (solid red) and the least-square method (dashed blue line). Strength of the impulse is E = 0 . in (b)and E = 0 . in (c). (d) Dependence of the estimation errors e of MS method (solid red line) and least-square method(dashed blue line) on strength E of the impulses injected to the FHN oscillator. r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A .................................................................. (a) -0.4 0 0.4 0.8 0 π 2π P h a s e s h i f t [r a d ] Phase [rad] (b) -15-10-5 0 5 10 0 π 2π y - v a r i a b l e Phase [rad] (c) -1.2-0.6 0 0.6 0 π 2π I n t e r ac ti on f un c ti on Phase difference [rad] Averaging methodMS method (d) F r e qu e n c y D i ff e r e n ce Coupling Strength Real simulationEstimated model
Figure 4. (a): Phase responses of chaotic dynamics observed from Rössler equations. By applying an impulse at variablephases, the phase shifts were measured as the difference in timing between the following peak of y -variable and the oneexpected from the average oscillation period. Bold black line represents the averaged phase response. (b): Waveforms of y -component of the Rössler equations. Bold red line represents the averaged waveform. (c) Coupling functions ˜ H ( ∆θ ) estimated by the present method (solid red line) and one (dashed blue line) obtained as the convolution of averagedphase response curve and the averaged waveform. (d) Synchronization diagram of the estimated phase model (solid redline) and the original coupled Rössler equations (dashed blue line). r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A .................................................................. (a) R R R C C L V DD V SS + - Functiongenerator (b) -300-200-100 0 100 200 300 0 π 2π P h a s e s e n s iti v it y f un c ti on Phase [rad]Experimental dataMS method (c) -40-20 0 20 40 0 π 2π C oup li ng F un c ti on Phase Difference [rad]Experimental dataMS method (d) F r e qu e n c y d i ff e r e n ce [ H z ] Coupling Strength ExperimentEstimated model
Figure 5.
Experiment of Van der Pol oscillator circuit. (a): Schematic illustration of the Van der Pol circuit, that is composedof an inductor ( L ), a capacitor ( C ), three resistors ( R , R , R ), an operational amplifier, and its associated powersupplies ( V DD , V SS ). External forcing is injected from a function generator (Keysight 33500B) through a capacitor( C ). (b): Phase sensitivity function estimated by the multiple-shooting method (red line) and the perturbation experiment(crosses). (c) Coupling functions ˜ H ( ∆θ ) estimated by the present method (solid red line) and one (dashed blue line)obtained by the averaging of the experimentally obtained phase sensitivity function and the sinusoidal input waveform. (d)Synchronization diagram of the estimated phase model (solid red line) and the experimental circuit system (dashed blueline). r X i v : . [ n li n . C D ] O c t Line 11: 1 Ritsumeikan University, Japan ⇒ ⇒ ⇒ network dynamicsLine 146: e.g. current ⇒ e.g., currentLine 151: the nonlinear ⇒ these nonlinearLine 153: P Nk =1 ∂f i ∂θ k ∂φ k ∂p j ⇒ P Nk =1 ∂f i ∂φ k ∂φ k ∂p j Line 186: oscillator, ω i gives ⇒ oscillator and ω i givesLine 187: i.e. ω i ⇒ i.e., ω i Line 188: determine ⇒ determinesLine 192: R π Z ( θ i ) G ( θ i , θ k ) dθ k ⇒ R π Z ( θ i + θ ′ ) G ( θ i + θ ′ , θ k + θ ′ ) dθ ′ Line 204: via following ⇒ via the followingLine 221: with the specific ⇒ with this specificLine 246: Please set the equal signs of (2.5) and (2.6) to locate horizontally atthe same position. ⇒ significantly (Figure 1d).Line 418: ˜ T , , ˜ T , = 0 ⇒ ˜ T , , ˜ T , Line 442: connectivity. Figure 2d, on the other hand, shows ⇒ connectivity. (start a new line) Figure 2d showsLine 458: oscillators, because ⇒ oscillators, sinceLine 464: R π Z ( θ ) G ( θ + ψ ) dψ ⇒ R π Z ( ψ ) G ( θ + ψ ) dψ Line 525: remove ”and”Line 530: is weak. ⇒ is weak in equation (4.1).Line 637–649: Please make the size of the table smaller.
Line 665: V = 6 V ⇒ V = 0 . ⇒ they are desynchronized with each other,Line 777–779: We have an experimental data generated from electric circuit. Wewould like to make these data available via a public repository, assoon as the manuscript is accepted for publication. Until then, wewould like to provide the data upon individual request. ⇒⇒