A practical method for estimating coupling functions in complex dynamical systems
Abstract
A foremost challenge in modern network science is the inverse problem on reconstruction (inference) of coupling equations and network topology from the measurements of the network’s dynamics. Of particular interest are the methods that can operate on real (empirical) data without interfering with the system. One of such earlier attempts [Tokuda et al., PRL, 2007] was a method suited for general limitcycle oscillators, yielding both oscillators’ natural frequencies and coupling functions between them (phase equations) from empirically measured time series. The present paper makes two advances along the same lines. First, it reviews the above method in a way comprehensive to domainscientists other than physics. Second, it presents three new extensions: (i) algorithm for inference of the phase sensitivity function, (ii) coupling function to approximate interaction among phasecoherent chaotic oscillators, (iii) analysis of the experimental data from a forced Van der Pol electric circuit. This reaffirms the range of applicability of the method and makes it accessible to a wider scientific community.
rsta.royalsocietypublishing.orgResearch
Article submitted to journal
Subject Areas:
Nonlinear dynamics, coupled oscillators, data analysis
Keywords:
Coupling function, phase dynamics, parameter estimation
Author for correspondence:
Isao T, Tokuda
A practical method for estimating coupling functions in complex dynamical systems
Isao T. Tokuda, Zoran Levnajic, and Kazuyoshi Ishimura
Ritsumeikan University, Japan
Faculty of Information Studies in Novo Mesto, Slovenia
Complex networks are representations of complex systems, where nodes (vertices) represent system’s units and links (edges) represent the interactions among those units [1, 2, 3, 4]. The functioning of a real network is a cumulative effect of its structure (topology of connections among nodes/units) and dynamics (interactions/relationships among these nodes) [3, 4]. Hence, in models of real networks, nodes are often assumed to be (simple) systems with their inherent
dynamics, whereas links mediate the dynamical coupling between the connected pairs of nodes. Using this paradigm, network science has made valuable contributions to all scientific disciplines that involve systems composed of many units, including biology, neuroscience, sociology, economics, etc. [1, 2, 3, 4, 5, 6].
To really grasp the functioning of a real network, we need information on both its structure and its dynamics. The inverse problem of reconstructing (or inferring) this information from the empirical data is a foremost challenge in modern network science. Namely, understanding the inner connectivity patterns of real networks not only enables us to grasp their operations, but also helps in controlling and predicting their behavior [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19].
The problem of network reconstruction can be seen as composed of two parts. The first part is the reconstruction of network topology, where one tries to learn which pairs of nodes are connected and which are not. This can (in some cases) be done separately from the second part of the problem, which is the reconstruction of the coupling functions that dictate how the connected nodes interact. Of course, two parts of the problem are inherently related, but which one to tackle first depends on what data are available, what assumptions can be reasonably made about the system, and what exactly we wish to learn.
Numerous reconstruction methods have been proposed over the last decade, both in physics [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19] and in computer science literature [16, 20, 21, 22, 23, 24, 25, 26, 27, 28]. While some methods tackle only one of the two above mentioned parts of the problem [10], other methods seek to address both parts at the same time. In physics literature, special emphasis is put on method aimed at oscillatory systems that exhibit synchronization as the most researched paradigm of collective dynamics. This includes methods focused on either topology, coupling functions, or both [7, 10, 8, 18].
On a somewhat different front, research effort was devoted to the problem of estimation of phase variables and phase equations from the data [29, 30, 31, 32, 33, 7, 34, 11, 35, 36]. Namely, isolated units in many real systems exhibit oscillatory nature, in the sense that they can be well approximated as limitcycle oscillators (oscillator whose dynamics after transients reduces to periodic or quasiperiodic orbit). Researchers showed that, even if the oscillatory behavior is very stochastic, there are robust ways to extract a well defined phase variable for each network node, and hence reconstruct the phase equations that describe the system dynamics. This paradigm found applications in diverse domain sciences, notably biology and neuroscience, where many systems have this nature. Estimating phase equations, however, is nothing but reconstructing coupling functions from data. While such a reconstruction approach is valid only in the approximation of phase variables, these methods require very little additional assumptions about the system. This means they can be almost immediately applied to empirical data [37, 29, 38, 7, 34, 39].
For a system of phase equations, a standard way to construct the coupling function is to measure phase sensitivity function of an individual oscillator element and obtain the coupling function by averaging method that computes the amount of phase shift induced through interaction with another oscillator element [40]. However, a precisely measured phase sensitivity function is not always available, since it requires application of external perturbations to an individual oscillator in an isolated condition [41, 42, 43, 44, 45, 46, 47, 48]. On the other hand, as a noninvasive approach, the coupling function can be inferred directly from timetrace data measured from coupled oscillators [29, 30, 31, 32, 38, 33, 7, 34, 39, 36]. One of them is a method by Tokuda et al. [29]: this approach utilized a multiple shooting method to realize robust parameter estimation of the coupling functions. The multiple shooting provides a general framework for fitting ordinary differential equations to recorded timetrace data. It is applicable to any system, where the dynamics of individual nodes can be approximated as those of limitcycle oscillators, yielding both oscillators’ natural frequencies and coupling functions between them (phase equations). Most importantly, the method was actually shown to operate very well with the data from a real experiment, which highlights its potential for practical use for physics problems and otherwise [29, 38, 49].
The contribution of the present paper is twofold. First, we review this method in a way that is understandable and approachable to communities outside physics. With this, we hope to make our method more useful to field such as biology and neuroscience, for which it was originally intended. Second, we show and discuss three new extensions of this method, specifically: (i) we extend this method to inference of the phase sensitivity function, which is vital for phase equations, (ii) coupling function is estimated for coupled chaotic oscillators to demonstrate how well the phase model approximates chaotic phase synchronization, (iii) using the data from a new experiment involving a system of Van der Pol electric circuits, we show another independent example of 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 original method in a comprehensive way. In section 3, we briefly discuss method’s extensions that were published so far. In section 4, we present novel contributions (extensions) mentioned above. In the last section, we discuss our findings and lay out perspectives for future work.
In this Section, we describe the original method in a more comprehensive way and show how it works for the case of coupled FitzHughNagumo oscillators.
We consider a network composed of interacting oscillator elements. In biology, such systems include a network of circadian cells in the suprachiasmatic nucleus [50], brain network composed of many spiking neurons [41, 44], cardiac muscle cells in the heart [51], etc. In terms of nonlinear dynamics, such systems are described as a system of coupled limit cycle oscillators:
(2.1) 
Here, and () represent state variables and autonomous dynamics of the th oscillator element, respectively. While represents an interaction function between th and th oscillators, strength of their interaction is determines by the coupling constant . The matrix describes connectivity between the oscillator elements. Without coupling (i.e., ), individual systems (i.e., ) generate periodic oscillations, after transients. Such oscillations are called limit cycles, which have intrinsic periods of . The equation (2.1) desciribes, to a good approximation, a variety of systems of interest in biology and neuroscience. Then the theory of phase reduction [52, 53] states that, as far as the coupling strength is weak enough, the network dynamics can be reduced to the following phase equations:
(2.2) 
where represents phase of the th oscillator, gives natural frequency of the th oscillator (i.e., ). stands for phase sensitivity function (also called “infinitesimal phase response curve”), which determine the amount of phase shift induced by the interaction with other oscillators (we will here not go in the detail of how equation (2.2) is obtained; an interested reader can refer to [52, 53, 54]). By averaging approximation [53], which integrates one cycle of the phase sensitivity function with the interaction function , the coupling function is derived as . Transformation of the original equations (2.1) to (2.2) provides a significant reduction in system’s dimensionality, in the sense that the original state variables , which can be highdimensional, are represented only by the single phase variable . This substantially simplifies the system’s modeling and enables its identification in a straightforward fashion.
As a recording data, individual oscillator states are simultaneously measured as (: observation function, : sampling time, : Data points). This corresponds to an experimental situation, under which states of individual oscillators (e.g., gene expression levels of individual cells, membrane potentials of neurons, etc.) are recorded simultaneously. Our purpose is to infer the phase equations from these measurement data under the conditions that the underlying system equations (2.1) are unknown.
The phase dynamics can be reconstructed in the following steps.

From the measured data , phases are extracted as , where represents the time, at which th signal takes its th peak and [54].

Fit the phases to the phase equations:
(2.3) which capture the underlying phase dynamics (2.2). represent approximate values for the natural frequencies. The coupling function , which is in general nonlinear and periodic with respect to , is approximated by a Fourier series up to order of as .
The unknown parameters are estimated by the multipleshooting method. The multiple shooting, developed in physics and engineering, provides a general framework for fitting ordinary differential equations to recorded time series data [55]. The methodology is applicable to a situation, under which the system equations and the recorded data are in a good quantitative agreement. As far as the condition of phase reduction, i.e., weakly coupled limit cycle oscillators, holds, the observed phase should obey the phase equations (2.3) in a quantitative manner. To estimate the parameters , time evolution of the phase equations (2.3) starting from an initial condition is denoted by . Then, at each sampling time , the phase equation must satisfy the boundary conditions: . With respect to the unknown parameters , the nonlinear equations are solved by the generalized Newton method. To compute the gradients , which are needed for the Newton method, variational equations are solved simultaneously, where represents th equation of the phase dynamics (2.3) as . The evolution function as well as the variational equations are integrated numerically.

To avoid overfitting of the coupling function, crossvalidation is utilized to determine the optimum number of Fourier series, [56]. We divide the measurement data into two parts. For the first half data, the parameter values are estimated. Then, the estimated parameters are applied to the latter half data and measure the error . The order number providing the minimum error is considered as the optimum.
To illustrate how the method described above works, we apply the multiple shooting to a prototypical example of weakly coupled limit cycle oscillators. In the original study [29], coupled Rössler oscillators were analyzed. As another yet challenging example, which has more complex shape of coupling function, we consider the following network of FitzHughNagumo (FHN) oscillators [57, 58]:
(2.4)  
(2.5) 
where . The system of FitzHughNagumo oscillators can be seen as a simple model for interacting neurons. Under the parameter setting of , , , , individual FHN oscillators (without coupling ) gives rise to limit cycles of slowfast type. Inhomogeneity parameters, which control natural periods of the individual oscillators, were set as (), where yields natural oscillation period of .
We started with the case of . As the coupling matrix, alltoall connections were assumed (). The level of inhomogeneity was set to . The multivariate data were recorded at a coupling strength of , which is in a nonsynchronized regime. The sampling interval was set to be . Then, the phases were extracted and downsampled to a sampling interval of . Total of data points have been collected for the parameter estimation. As an initial condition, the unknown parameter values were all set to be zero, i.e., , . The data were divided into and , which were used for the parameter estimation and the crossvalidation error , respectively. By varying the Fourier order from to , the optimal value was found to be by the crossvalidation test.
The coupling function estimated by the present method is in a good agreement with the one computed by the adjoint method [59] (Fig. 1a). The errorbars were computed from inverse of the Hessian matrix of the squared error function, under the assumption that the phase data contain uncorrelated observational noise [60]. The estimated natural frequencies are distributed on 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 the original coupled FHN oscillators can be recovered, where the onset of synchronization was predicted at , which is close to the real onset of (Fig. 1c).
Next, we show how the estimation depends upon the problem setting. The primary factor that influences the estimation results is the coupling strength used to generate the time series. Fig. 1 (d) shows dependence of estimation error on the coupling strength. The estimation error is evaluated as deviation of the estimated coupling function from the one estimated by the adjoint method, i.e.,
(2.6) 
where the denominator represents normalization factor and . As the coupling strength is located close to the onset of synchronization, the estimation error increases significantly. Under the synchronized state, phase differences between the oscillators do not change in time , providing no information about the phase interaction. Increase in estimation error due to synchronized data is therefore reasonable.
Even in a synchronized regime, the coupling function can be recovered from transient data, during which phase differences evolve (transient data often reveals far more information about the underlying system, since it is recorded before the system ’settled’ into its dynamical equilibrium). To show this, the multivariate data were recorded after discarding only a short duration of transient process that starts from a random initial condition. Before the system is confined in a synchronized state, transient process of data (time interval of ) was collected. sets of such data were used for the parameter estimation. Fig. 1 (e) shows dependence of the estimation error on the transient duration. Note that the coupling strength is set to , that is in a 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 the system is in synchrony with a moderate coupling, application of perturbation that kicks the system 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 , varied from to . The level of inhomogeneity was set to . The multivariate data were recorded at a coupling strength of , which corresponds to nonsynchronized regime. Other settings were the same as those in the case of . Surprisingly, the estimation error remains in a low level. Even for , the coupling function has been precisely estimated, while the estimated natural frequencies are consistent with those obtained from the noncoupled simulations. This suggests that the system size does not pose a major limit on the estimation of phase dynamics as far as the data contain nonsynchronized phase information.
In this section, we discuss extensions of the original method so far. Although we have dealt with the case that all oscillator elements are connected to all the others in the previous section, heterogeneous connections are more common in nature and engineering. Another challenge of our technique is to infer the connectivity of the oscillator network from the measured data. Numerous approaches have been proposed up to date using information transfer [61], mutual predictability [62], recurrence properties [63], permutationbased asymmetric association measure [64], index for partial phase synchronization [65, 66, 67], and graph theory [68]. Response properties of the network dynamics to external stimuli have been also exploited [69, 8]. For weakly coupled limit cycle oscillators, to which phase reduction is applicable, the phase modeling approach is again quite effective for detecting the network topology [49, 11, 70, 71, 72]. As an application of the coupling function estimated in the previous section, multipleshooting method was extended to infer connectivity of an oscillator network [49]. Here, we present such method. For simplicity, we suppose that the connection matrix is composed of zeroorunity elements (i.e., or ). The estimation procedure is the same as before except that the connection matrix is estimated as the unknown parameters in the multiple shooting method. The coupling function was fixed to the one estimated in the previous estimation. Natural frequencies obtained from noncoupled original equations were also utilized in the phase model.
For a network of four () oscillators, two defects, i.e., , were introduced to the connection matrix. The level of inhomogeneity was set to , whereas the coupling strength was , i.e., in a nonsynchronized regime. As in the previous section, total of data points (sampling time: ) have been collected. The connection matrix was estimated as follows.
We see that the two defects were precisely identified as small values, whereas other matrix elements pointed almost unity.
For comprehensive analysis, the connection matrices with randomly generated defects were estimated for variable network size from to . For our analysis, the estimation error was evaluated as , where the estimated connectivity was digitized as for and otherwise. For each setting of the network size, instances of connection matrices were randomly generated and the average and the standard deviation of the estimation errors were plotted in Fig. 2. Panels (a) and (b) correspond to 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 are almost zero except 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 estimated by using a simple sine as the coupling function, i.e., . For small networks, the difference was not large between the precisely estimated (higherorder) coupling function (red solid line) and the simple sine function (blue dotted line). However, as the network size increases, the estimation error increases much more rapidly in the sine function than in the higherorder coupling function. This indicates that, for reliable detection of the connectivities, precisely estimated coupling function is of significant importance.
Finally, in this section, we present three new extensions of our original method. One subsection is devoted to each.
The phase sensitivity function plays a vital role in the studies of coupled oscillators [52, 53, 54]. Numerous approaches have been proposed to estimate the phase sensitivity function from experimental data [41, 42, 43, 44, 45, 46, 47, 48]. As an extension of our technique, the phase sensitivity function can be recovered from the coupling function [73]. As described earlier in the averaging approximation [53], the coupling function is given by a convolution of the phase sensitivity function and the input waveform as . It is straightforward to recover the phase sensitivity function by the spectral deconvolution [74]. Namely, in a frequency domain, the phase sensitivity function is given as , where , , and represent Fourier transforms of , , and , respectively. Inverse Fourier transform of yields the phase sensitivity . Fig. 3(a) shows phase sensitivity function (solid red line) obtained by the deconvolution of the coupling function estimated from coupled FHN oscillators () in section 2. Compared with the one computed by the adjoint method [59] (dashed blue line), the estimated function is somewhat deviated from the true curve. We consider that, due to 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 the Winfree formula [52] as follows. For simplicity, we consider a single phase oscillator receiving an external force :
(4.1) 
where represents natural frequency of the oscillator. The phase sensitivity function is described in terms of a Fourier series as . The unknown coefficients can be estimated by the multipleshooting method in a similar manner as the estimation of coupling function. We apply this technique to a single FHN oscillator that receives random impulses (stimulus duration: , stimulus strength: ) as external forcing . Parameter values of the FHN oscillator and the sampling time interval were set to be the same as those in the previous sections. For simplicity, natural frequency and the external signal were assumed to be known. Order of the Fourier series was set to . For impulse strength of and , the estimated phase sensitivity functions are plotted by solid red line in Fig. 3(b) and (c), respectively. As a comparison, estimation results of the leastsquare method [41, 44], applied to the same data sets as one of the standard estimation techniques, are plotted by dashed blue lines. In both panels (b) and (c), estimation results of the multiple shooting method are consistent with those of the adjoint method [59]. The leastsquare method, on the other hand, recovered the 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 the impulse strength is increased, the estimation error increases much more rapidly in the leastsquare method (dashed blue line) than the multiple shooting method (solid red line) (see Fig. 3d). The leastsquare method [41, 44] assumes that phase of the oscillator evolves linearly in time according to the natural frequency. Consequently, strong perturbation, which induces nonsmall phase shift, may increase the modeling error. On the other hand, in the present formula (4.1), the phase equation takes into account the phase shift induced by the external perturbations, reducing the estimation error.
As another application, the estimated coupling function can be utilized for modeling chaotic phase synchronization [75]. It has been known that phases of chaotic oscillators can be synchronized with each other, while their amplitudes remain irregular and uncorrelated. Especially for phasecoherent chaos, in which rotation center can be welldefined, the phase dynamics can be approximated as , where represents frequency modulation, which depends upon oscillation amplitude [75]. For chaotic amplitude , the term can be regarded as an effective noise. In many phasecoherent systems such as the Rössler equations [76], amplitudedependent frequency modulation is very small, so the noise term is negligible. Phase dynamics of such chaotic attractor becomes very similar to those of limit cycle oscillators.
To extract phaseinteraction between chaotic oscillators, we consider two coupled Rössler equations [76]:
Each Rössler oscillator gives rise to chaotic dynamics without coupling . The inhomogeneity parameters were set as , which yield average oscillation periods of and , respectively. The bivariate data were simulated under coupling strength of , which corresponds to nonsynchronized regime. The sampling interval was set to be for the extraction of the phases . Then, to apply the multipleshooting method, the data have been down sampled to and the total of data points were collected. The data were divided into and points, which were used for the parameter estimation and the crossvalidation test, respectively. By varying the Fourier order from to , the optimal value was found to be . The corresponding coupling function is shown by the solid red line in Fig. 4(c). The estimated function is in good agreement with the one obtained by the convolution of averaged phase sensitivity function (Fig. 4a) and the averaged input waveform (Fig. 4b). Using the estimated phase equations, synchronization diagram of the original two coupled Rössler equations can be recovered, where the onset of synchronization was predicted at , which is close to the real onset of (Fig. 4d). This suggests that our simple method of estimating the coupling function provides a good approximate of describing the phase dynamics of phasecoherent chaotic oscillators.
Finally, we apply our method to experimental data generated from Van der Pol electric circuit [77]. As shown in the diagram of Fig. 5a, the system is based on a LC circuit, composed of an inductor () and a capacitor (). To form a negativeresistance converter, three positive resistors (, , ) were connected to a voltagecontrolled voltage source (i.e., operational amplifier and its associated power supplies , ) [78]. External forcing was injected from a function generator (Keysight 33500B) to the Van der Pol circuit through a capacitor (). Physical parameters of the electric components used in the present experiment are summarized in Table 1. To obtain the phase sensitivity function, impulses (stimulus duration: s, stimulus strength: V) were randomly injected as the external force . The circuit output as well as the input impulses were simultaneously measured with a sampling frequency of kHz. First, the phase sensitivity function was estimated by fitting the measured data to equation (4.1) with the multipleshooting method. Natural frequency Hz (i.e., ), measured before the stimulus experiment, was used in the phase dynamics. Order of the Fourier series was set to . As shown in Fig. 5b, the estimated phase sensitivity fits to the experimental observation of phase response data well.
Next, a sinusoidal forcing (forcing frequency: Hz, forcing amplitude: V) was applied to the Van der Pol circuit. The circuit output as well as the forcing waveforms were simultaneously measured with a sampling frequency of kHz. By the multipleshooting method, which fits the measured data to the phase equations (2.3), the coupling function (Fourier order: ) was estimated. In Fig. 5c, the estimated coupling function is compared with the one obtained by the averaging of the phase sensitivity function , estimated from the impulse stimuli, and the input sine waveform . Despite a slight difference in the initial phase, the coupling functions agree quite well with each other. In Fig. 5d, the estimated phase equations recovered synchronization diagram of the experimental system, where the onset of synchronization was predicted at V, which is very close to the real onset of V.
[mH]  
[uF]  
[K]  
[k]  
[k]  
5 [V]  
5 [V]  
OPAMP  LF412CN 
[nF] 
The multipleshooting method has been focused on as a noninvasive approach to estimate coupling functions from observed multivariate time series [29]. Among various methods developed so far [37, 30, 31, 32, 38, 33, 7, 34, 39, 36], which are based on the Bayesian estimation and other variants, the multipleshooting provides a straightforward approach to fit the phase equations to phase data measured from an oscillator network. Despite its simplicity, the method was shown to be capable of precisely estimating the coupling function of the coupled FHN oscillators including higherharmonic terms. The estimation was found effective for a large network of up to oscillators. Utilization of the transient part of data successfully enlarged applicability of the estimation technique even in a synchronized regime of coupled oscillators. The estimated coupling function was further applied to inference of network topology and chaotic phase synchrony. Precise estimation of the coupling functions was shown to improve the reconstruction of network topology. As another intriguing issue, estimation of the phase sensitivity function was also discussed. Although the phase sensitivity function obtained by deconvolution of the estimated coupling function was slightly deviated from the true function, refinement has been made by extending the multiple shooting method directly to the phase data of a driven limit cycle oscillator. Finally, efficiency of the present approach was demonstrated with the experimental data measured from the Van der Pol electric circuit with a sinusoidal forcing. To our knowledge, phase response property of the Van der Pol circuit has not been measured experimentally, adding a new insight into the experimental system.
Beyond experimental system in physics, chemistry, and engineering, we foresee that our method will be applicable to system of rhythmic, interacting elements such as cellular gene expressions in the suprachiasmatic nucleus (SCN) [50], electrical activities of cardiac pacemakers [51], inferior olive neurons in the cerebellum [79] and can give insights useful for domainscientists in biology and neuroscience.
While considering our method of potentially practical use for various systems, its usefulness is not without limitations. The main among them is the assumption that the studied system can be approximated as a network of weakly coupled limit cycles [53]. This, however, is not true for all systems encountered in nature. For instance, in gene regulatory networks, phases of the clock component genes are tightly connected to each other [80]. It has been known that cortical neurons fire with a strong synchrony during epileptic seizure [81]. Such strongly coupled systems should be carefully distinguished and avoided as a target of modeling the phase dynamics. Even for tightly connected systems, there still remains a possibility that our method could, to some approximation, be implemented to model their transient dynamical process, where phase relationship between the oscillator components may change dynamically after strong external stimuli [82].
Another limitation is the length of the available time series: namely, experimental measurements, for a variety of realistic reasons, could produce the data (time series) of only a very short length. For instance, time resolved data on gene regulation are not likely to yield time series with much more than 10 cycles. In this case, our method might be of limited use. Also, realistic data are almost always noisy. The noise strength, depending on the experimental scenario, could be quite severe. Especially, the phase extraction process in our modeling is rather sensitive to noise. Temporal fluctuation and noise in natural frequencies of the oscillator elements may also cause estimation error in the coupling functions. In this respect, noise tolerance should be carefully examined, before the application to data contaminated with observational/dynamical noise.
Also, networks in real world are large and only partials of the dynamics elements are observable. Although our method was shown to be robust against system size as far as the oscillator elements are uniformly connected and they desynchronized, the effect of unobserved oscillator states should be examined carefully. Heterogeneity and hierarchy in the coupling functions may require further extension of the present approach.
Finally, we conclude the paper with a brief discussion of how our method’s performance compares to the performance of other methods that reconstruct coupling functions in oscillatory systems. Unfortunately, such comparison is not simple to make, since various available methods depart from very different hypotheses and knowledge about the system. Stronger hypotheses lead to better inferences, but the information on whether the hypotheses are met is not always available. This renders hard any independent comparison of reconstruction methods. One could argue that methods aimed at only network topology are more useful and precise, but such methods neglect the entire dynamical nature of many real networks. On the other hand, certain methods give excellent results, but are limited to dynamical systems with specific properties. In fact, our method belong to this category, since it assumes the limit cycle nature of the individual units. Furthermore, methods can be divided into invasive ones (that interfere with system’s ongoing dynamics) and noninvasive ones (that do not). Again, their real merits is hard to compare, since invasive methods, although often non practical, will almost always give better results. Therefore, we here conclude that our reconstruction concept, although limited by the assumption of limit cycles, is a promising – and above all practical – approach implementable in real experiments.
Authors’ Contributions. ITT designed the study. ITT performed the numerical simulations and the data analysis. KI carried out the circuit experiments. ITT and ZL wrote the manuscript. All authors read and approved the manuscript.
Competing Interests. The author(s) declare that they have no competing interests.
Funding. ITT acknowledges support by GrantinAid for Scientific Research (No. 17H06313, No. 16H04848, No. 16K00343, No. 18H02477) from Japan Society for the Promotion of Science (JSPS). ZL acknowledges support by EU via H2020MSCAITN2015 COSMOS 642563, and by Slovenian research agency via P10383 and J58236.
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 Connected World. Cambridge University Press.
 6 Boccaletti S, Bianconi G, Criado R, Del Genio CI, GómezGardenes J, Romance M, SendinaNadal I, Wang Z, Zanin M. 2014 The structure and dynamics of multilayer networks. Physics Reports 544, 1–122.
 7 Blaha KA, Pikovsky A, Rosenblum M, Clark MT, Rusin CG, Hudson JL. 2011 Reconstruction of twodimensional phase dynamics from experiments on coupled oscillators. Physical Review E 84, 046201.
 8 Levnajić Z, Pikovsky A. 2011 Network reconstruction from random phase resetting. Physical review letters 107, 034101.
 9 Wang WX, Yang R, Lai YC, Kovanis V, Harrison MAF. 2011 Timeseries–based prediction of complex oscillator networks via compressive sensing. EPL (Europhysics Letters) 94, 48006.
 10 Stankovski T, Duggento A, McClintock PV, Stefanovska A. 2012 Inference of timeevolving coupled dynamical systems in the presence of noise. Physical review letters 109, 024101.
 11 Kralemann B, Pikovsky A, Rosenblum M. 2014 Reconstructing effective phase connectivity of oscillator networks from observations. New Journal of Physics 16, 085013.
 12 Levnajić Z, Pikovsky A. 2014 Untangling complex dynamical systems via derivativevariable correlations. Scientific reports 4, 5030.
 13 Timme M, Casadiego J. 2014 Revealing networks from dynamics: an introduction. Journal of Physics A: Mathematical and Theoretical 47, 343001.
 14 Han X, Shen Z, Wang WX, Di Z. 2015 Robust reconstruction of complex networks from sparse data. Physical review letters 114, 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 and complex dynamical systems. Physics Reports 644, 1–76.
 17 Leguia MG, Andrzejak RG, Levnajić Z. 2017 Evolutionary optimization of network reconstruction from derivativevariable correlations. Journal of Physics A: Mathematical and Theoretical 50, 334001.
 18 Rosenblum M, et al. 2017 Reconstructing networks of pulsecoupled oscillators from spike trains. Physical Review E 96, 012209.
 19 Simidjievski N, Tanevski J, Ženko B, Levnajić Z, Todorovski L, Džeroski S. 2018 Decoupling approximation robustly reconstructs directed dynamical networks. New Journal of Physics 20, 113003.
 20 Zanin M, Papo D, Sousa PA, Menasalvas E, Nicchi A, Kubik E, Boccaletti S. 2016 Combining complex networks and data mining: why and how. Physics Reports 635, 1–44.
 21 Bridewell W, Langley P, Todorovski L, Džeroski S. 2008 Inductive process modeling. Mach. Learn. 71, 1–32.
 22 Džeroski S, Todorovski L. 2008 Equation discovery for systems biology: finding the structure and dynamics of biological networks from time course data. Current Opinion in Biotechnology 19, 360–368.
 23 Schmidt M, Lipson H. 2009 Distilling freeform natural laws from experimental data. science 324, 81–85.
 24 Brunton SL, Proctor JL, Kutz JN. 2016 Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the National Academy of Sciences 113, 3932–3937.
 25 Simidjievski N, Todorovski L, Džeroski S. 2016 Modeling dynamic systems with efficient ensembles of processbased models. PLoS One 11, e0153507.
 26 Tanevski J, Todorovski L, Džeroski S. 2016 Learning stochastic processbased models of dynamical systems from knowledge and data. BMC systems biology 10, 30.
 27 Tanevski J, Todorovski L, Džeroski S. 2016 Processbased design of dynamical biological systems. Scientific reports 6, 34107.
 28 Čerepnalkoski D, Taškova K, Todorovski L, Atanasova N, Džeroski S. 2012 The influence of parameter fitting methods on model structure selection in automated modeling of aquatic ecosystems. Ecological Modelling 245, 136–165.
 29 Tokuda IT, Jain S, Kiss IZ, Hudson JL. 2007 Inferring phase equations from multivariate time series. Physical review letters 99, 064101.
 30 Kralemann B, Cimponeriu L, Rosenblum M, Pikovsky A, Mrowka R. 2007 Uncovering interaction of coupled oscillators from data. Physical Review E 76, 055201.
 31 Kralemann B, Cimponeriu L, Rosenblum M, Pikovsky A, Mrowka R. 2008 Phase dynamics of coupled oscillators reconstructed from data. Physical Review E 77, 066205.
 32 Cadieu CF, Koepsell K. 2010 Phase coupling estimation from multivariate phase statistics. Neural computation 22, 3107–3126.
 33 Kralemann B, Pikovsky A, Rosenblum M. 2011 Reconstructing phase dynamics of oscillator networks. Chaos: An Interdisciplinary Journal of Nonlinear Science 21, 025104.
 34 Zhu Y, Hsieh YH, Dhingra RR, Dick TE, Jacono FJ, Galán RF. 2013 Quantifying interactions between real oscillators with information theory and phase models: application to cardiorespiratory coupling. Physical Review E 87, 022709.
 35 Pikovsky A. 2018 Reconstruction of a random phase dynamics network from observations. Physics Letters A 382, 147–152.
 36 Suzuki K, Aoyagi T, Kitano K. 2018 Bayesian estimation of phase dynamics based on partially sampled spikes generated by realistic model neurons. Frontiers in computational neuroscience 11, 116.
 37 Miyazaki J, Kinoshita S. 2006 Determination of a coupling function in multicoupled oscillators. Physical review letters 96, 194101.
 38 Tokuda IT, Wagemakers A, Sanjuán MA. 2010 Predicting the synchronization of a network of electronic repressilators. International Journal of Bifurcation and Chaos 20, 1751–1760.
 39 Stankovski T, Pereira T, McClintock PV, Stefanovska A. 2017 Coupling functions: universal insights into dynamical interaction mechanisms. Reviews of Modern Physics 89, 045001.
 40 Kiss IZ, Zhai Y, Hudson JL. 2005 Predicting mutual entrainment of oscillators with experimentbased phase models. Physical review letters 94, 248301.
 41 Galán RF, Ermentrout GB, Urban NN. 2005 Efficient estimation of phaseresetting curves in real neurons and its significance for neuralnetwork modeling. Physical review letters 94, 158101.
 42 Ota K, Omori T, Aonishi T. 2009 Map estimation algorithm for phase response curves based on analysis of the observation process. Journal of computational neuroscience 26, 185.
 43 Nakae K, Iba Y, Tsubo Y, Fukai T, Aoyagi T. 2010 Bayesian estimation of phase response curves. Neural Networks 23, 752–763.
 44 Hong S, Robberechts Q, Schutter ED. 2012 Efficient estimation of phaseresponse curves via compressive sensing. Journal of neurophysiology 108, 2069–2081.
 45 Goldberg JA, Atherton JF, Surmeier DJ. 2013 Spectral reconstruction of phase response curves reveals the synchronization properties of mouse globus pallidus neurons. Journal of neurophysiology 110, 2497–2506.
 46 Saifee TA, Edwards MJ, Kassavetis P, Gilbertson T. 2015 Estimation of the phase response curve from parkinsonian tremor. Journal of neurophysiology 115, 310–323.
 47 Imai T, Ota K, Aoyagi T. 2017 Robust measurements of phase response curves realized via multicycle weighted spiketriggered averages. Journal of the Physical Society of Japan 86, 024009.
 48 Cestnik R, Rosenblum M. 2018 Inferring the phase response curve from observation of a continuously perturbed oscillator. Scientific Reports 8, 13606.
 49 Tokuda IT, Wickramasinghe M, Kiss IZ. 2013 Detecting connectivity of small, dense oscillator networks from dynamical measurements based on a phase modeling approach. Physics Letters A 377, 1862–1867.
 50 Yamaguchi S, Isejima H, Matsuo T, Okura R, Yagita K, Kobayashi M, Okamura Hs. 2003 Synchronization of Cellular Clocks in the Suprachiasmatic Nucleus. Science 302, 1408–1412.
 51 Verheijck EE, Wilders R, Joyner RW, Golod DA, Kumar R, Jongsma HJ, Bouman LN, van Ginneken AC. 1998 Pacemaker synchronization of electrically coupled rabbit sinoatrial node cells. The Journal of general physiology 111, 95–112.
 52 Winfree AT. 2001 The geometry of biological time. Springer, New York.
 53 Kuramoto Y. 1984 Chemical oscillations, waves, and turbulence. Springer, Berlin.
 54 Pikovsky A, Rosenblum M, Kurths J. 2003 Synchronization: a universal concept in nonlinear sciences, volume 12. Cambridge university press.
 55 Baake E, Baake M, Bock H, Briggs K. 1992 Fitting ordinary differential equations to chaotic data. Physical Review A 45, 5524.
 56 Stone M. 1977 An asymptotic equivalence of choice of model by crossvalidation and akaike’s criterion. Journal of the Royal Statistical Society. Series B (Methodological) pp. 44–47.
 57 FitzHugh R. 1961 Impulses and physiological states in theoretical models of nerve membrane. Biophysical journal 1, 445–466.
 58 Nagumo J, Arimoto S, Yoshizawa S. 1962 An active pulse transmission line simulating nerve axon. Proceedings of the IRE 50, 2061–2070.
 59 Ermentrout B. 1996 Type i membranes, phase resetting curves, and synchrony. Neural computation 8, 979–1001.
 60 Horbelt W, Timmer J, Bünner M, Meucci R, Ciofini M. 2001 Identifying physical properties of a co 2 laser by dynamical modeling of measured time series. Physical Review E 64, 016222.
 61 Schreiber T. 2000 Measuring information transfer. Physical review letters 85, 461.
 62 Rosenblum MG, Pikovsky AS. 2001 Detecting direction of coupling in interacting oscillators. Physical Review E 64, 045202.
 63 Romano MC, Thiel M, Kurths J, Grebogi C. 2007 Estimation of the direction of the coupling by conditional probabilities of recurrence. Physical Review E 76, 036211.
 64 Hempel S, Koseska A, Kurths J, Nikoloski Z. 2011 Inner composition alignment for inferring directed networks from short time series. Physical review letters 107, 054101.
 65 Schelter B, Winterhalder M, Dahlhaus R, Kurths J, Timmer J. 2006 Partial phase synchronization for multivariate synchronizing systems. Physical review letters 96, 208103.
 66 Nawrath J, Romano MC, Thiel M, Kiss IZ, Wickramasinghe M, Timmer J, Kurths J, Schelter B. 2010 Distinguishing direct from indirect interactions in oscillatory networks with multiple time scales. Physical review letters 104, 038701.
 67 Wickramasinghe M, Kiss IZ. 2011 Phase synchronization of three locally coupled chaotic electrochemical oscillators: Enhanced phase diffusion and identification of indirect coupling. Physical Review E 83, 016210.
 68 Runge J, Heitzig J, Petoukhov V, Kurths J. 2012 Escaping the curse of dimensionality in estimating multivariate transfer entropy. Physical review letters 108, 258701.
 69 Timme M. 2007 Revealing network connectivity from response dynamics. Physical review letters 98, 224101.
 70 Stankovski T, Ticcinelli V, McClintock PV, Stefanovska A. 2015 Coupling functions in networks of oscillators. New Journal of Physics 17, 035002.
 71 Stankovski T, Petkoski S, Raeder J, Smith AF, McClintock PV, Stefanovska A. 2016 Alterations in the coupling functions between cortical and cardiorespiratory oscillations due to anaesthesia with propofol and sevoflurane. Phil. Trans. R. Soc. A 374, 20150186.
 72 Stankovski T, Ticcinelli V, McClintock PV, Stefanovska A. 2017 Neural crossfrequency coupling functions. Frontiers in systems neuroscience 11, 33.
 73 Kralemann B, Pikovsky A, Rosenblum M. 2013 Detecting triplet locking by triplet synchronization indices. Physical Review E 87, 052904.
 74 Wiener N. 1949 Extrapolation, interpolation, and smoothing of stationary time series: with engineering applications. MIT press Cambridge, MA.
 75 Rosenblum MG, Pikovsky AS, Kurths J. 1996 Phase synchronization of chaotic oscillators. Physical review letters 76, 1804.
 76 Rössler OE. 1976 An equation for continuous chaos. Physics Letters A 57, 397–398.
 77 Van der Pol B, Van Der Mark J. 1927 Frequency demultiplication. Nature 120, 363.
 78 Kennedy MP. 1992 Robust op amp realization of chua’s circuit. Frequenz 46, 66–80.
 79 Tokuda IT, Hoang H, Kawato M. 2017 New insights into olivocerebellar circuits for learning from a small training sample. Current opinion in neurobiology 46, 58–67.
 80 Zhang R, Lahens NF, Ballance HI, Hughes ME, Hogenesch JB. 2014 A circadian gene expression atlas in mammals: implications for biology and medicine. Proceedings of the National Academy of Sciences 111, 16219–16224.
 81 Velazquez JLP, Carlen PL. 2000 Gap junctions, synchrony and seizures. Trends in neurosciences 23, 68–74.
 82 Schmal C, Ono D, Myung J, Pett JP, Honma S, Honma KI, Herzel H, Tokuda IT. 2019 Weak coupling between intracellular feedback loops explains dissociation of clock gene dynamics. bioRxiv p. 542852.