Noninvasive inference methods for interaction and noise intensities of coupled oscillators using only spike time data.

Measurements of interaction intensity are generally achieved by observing responses to perturbations. In biological and chemical systems, external stimuli tend to deteriorate their inherent nature, and thus, it is necessary to develop noninvasive inference methods. In this paper, we propose theoretical methods to infer coupling strength and noise intensity simultaneously in two well-synchronized noisy oscillators through observations of spontaneously fluctuating events such as neural spikes. A phase oscillator model is applied to derive formulae relating each of the parameters to spike time statistics. Using these formulae, each parameter is inferred from a specific set of statistics. We verify these methods using the FitzHugh-Nagumo model as well as the phase model. Our methods do not require external perturbations and thus can be applied to various experimental systems.

Coupled oscillators such as cardiac myocytes (1), heart pacemakers (2, 3), circadian clocks (2–5), neurons (6–9), electrochemical oscillators (10), spin torque oscillators (11–14), crystal oscillators (15), and nanomechanical oscillators (16) are found in many disciplines ranging from biology to engineering. Although these systems are subject to various types of noise, including thermal, quantum, and molecular noise, they can exhibit synchronization because of coupling between the oscillators. Thus, coupling and noise are crucial factors in the determination of multioscillator dynamics (17–19).

Since a noninvasive estimation is desired in many cases, it is important to develop methods to infer the coupling strength and noise intensity solely from temporal information on the oscillation. Such an attempt was made in an experiment with cultured cardiac myocytes beating spontaneously (1). Therein, the transition from a desynchronized state to a synchronized state between two cells was observed within the incubation time. This suggests that coupling between the cells should increase. However, this naive expectation is not generally fulfilled, because synchronization is facilitated not only by increased coupling strength but also by decreased noise intensity (17).

Fig. 1 displays spike time data generated with the FitzHugh–Nagumo model for cardiac and neural electrical activity (20) (precisely introduced later). For parameter sets i and ii, the typical values ζ of the spike time lag, which represent the degree of synchronization, are approximately equal. From this, the coupling strengths in the two cases may seem similar. However, the values actually differ by a factor of 2. Thus, an individual statistic derived from oscillation data can be misleading when attempting to infer coupling strength. The case of attempting to infer noise intensity is similar. Hence, in order to infer these properties, different types of statistics must be combined appropriately.

Fig. 1.

(A) Examples of spike timing generated for coupled cells with the FitzHugh–Nagumo model in Eq. . The SDs, ζ, of the spike time lag between two oscillators are similar in cases i and ii. (B) Simultaneous inferences of effective noise intensity aD and effective coupling strength for the FitzHugh–Nagumo model with method II. These inferences were achieved using only spike time data. Actual values are plotted as crosses. Inferred values are plotted as squares (for the lowest and third-lowest coupling strengths) and circles (otherwise). The actual values are very well approximated in all cases, including D = 0.0002, (case i) and D = 0.0004, (case ii). The average period is .

In this paper, we propose two methods of inferring coupling strength and noise intensity from data solely based upon the spike timing of two well-synchronized noisy oscillators. Method I requires spike timing data from only one of the oscillators, but we may infer the coupling strength and the noise intensity. Method II requires spike time data on both oscillators but provides more precise inferences. We demonstrate our methods using phase oscillator and FitzHugh–Nagumo models. An example of our inferences from the FitzHugh–Nagumo model is shown in Fig. 1. There, the coefficient of variation in periods (1.9 to 4.4%) and the number of observed spikes (160,000) are comparable to those in the abovementioned experiment on cardiac cells (1); hence, the demonstration in the figure is realistic. Although many inference methods work effectively using data taken from unsynchronized oscillators (21–23), external perturbations (24–26), and whole time series (21–23, 27, 28), ours do not require them. Moreover, we do not need to assume function form. Therefore, our methods are ready for application to synchronized coupled oscillator systems in various fields.

Formulations.

We introduce some statistical quantities based on the spike time data (Fig. 1). We assume that an oscillator spikes when its oscillatory variable passes a specific value. Let us denote the kth spike time of the oscillator by . In the case of the phase oscillator model, is defined as the time at which a phase first passes through , where is called the checkpoint phase. The kth m-cycles period and its variance are defined aswhere denotes the statistical average over k and τ is the average period given by . Note that in this paper, denotes both the statistical average over k and the ensemble average, which are identical in the steady state. Note further that V is calculated from the spike time data of one oscillator. To quantify the relationship between two oscillators, the SD of time lag between the spikes of the oscillators is defined aswhere is the kth spike time of the ith oscillator. It should be noted that in coupled noisy oscillators, phase slips can generally occur; this occurrence of phase slips demands redefinition of the time lag. To adequately define the time lag after the phase slips, we modify the definition of the th spike time when phase slips occur after the kth spikes of oscillator i and j as follows: when oscillator i generates successive spikes before the th spike generated by oscillator j, i.e., holds true, we regard and as the corresponding spikes; i.e., we redefine .

To derive an inference theory, we consider a pair of coupled phase oscillators subject to noise. When limit cycle oscillators are weakly coupled to each other and subject to weak noise, the dynamics can be described by (17, 29)where θ is the phase of oscillator i and is the coupling strength. The independent and identically distributed noise satisfies and . The positive constant D represents the noise intensity. The phase sensitivity function is a -periodic function that quantifies the phase response to noise. The -periodic function J(x, y) describes the interaction between oscillators that leads to synchronization. We assume that , which is satisfied in systems with diffusive coupling between chemical oscillators or gap junction coupling between cells. We focus on systems that are well synchronized in phase.

Our inference methods are based on the formula of period variability. In a previous work (30), the following expression for the variance V1 was derived from the system in Eq. 4 by means of linear approximation:where C1 and C2 are independent of and given by and . The negative constant corresponds to the average effective attractive force between the oscillators over one oscillation period. That is, , where . The -periodic function represents the phase distance from in-phase synchronization, where is the phase difference defined on the ring . If x is the value of x(t) when θ1 first passes through , then represents the average of x over k. Note that is proportional to and dependent on κ (30).

Through a derivation similar to that of Eq. , we derive that V is given bywhere . See for the derivation. Since a represents an average phase response to noise, the product aD represents the effective noise intensity (17). Our purpose is now to infer aD and , which are important values because they determine the strength of the phase diffusion and the time scale of the synchronization, respectively (17).

Method I.

We use only V1, V2, and V3 for one of the oscillators. Combining Eq. for m = 1, 2, 3, we can determine the three unknowns aD, , and . In particular, we obtainand

Note that as shown below, Eq. states that a temporal correlation decays exponentially with spike times, and the decay constant is given by the effective coupling strength . We define the temporal correlation as

Recall that . When n is sufficiently large, i.e., and for any k, we obtainwhere . Thus, the numerator and the denominator in Eq. represent the correlations G1 and G2, respectively, i.e., .

Method II.

We additionally use ζ, which is the SD of the spike time lags. When D and κ are sufficiently small, we can assume that . Using this approximation, we can express V aswhere is neglected. In terms of V1, V2, and ζ, the two unknowns aD and are given byand

Our formulae in Eqs. 7, , , and are independent of the checkpoint phase, whereas V and ζ are not (30).

Numerical Results.

We demonstrated the validity of the inference methods with numerical experiments. First, we again employed the phase oscillator model in Eq. . We assumed , which represents gap junction coupling or diffusive coupling (17, 29). We set for and for , with . The region satisfying mimics the refractory stage that exists for many chemical and biological oscillators. We set and . Under these assumptions, , and τ = 1. For and we assumed white Gaussian noise.

We prepared 16 parameter sets, each with a different combination of coupling strength and noise intensity given by and , where . We integrated Eq. using the Euler scheme with a time step of . The initial conditions were . In this simulation, we fixed the checkpoint phase at and observed the spike timing for . Three realizations were simulated for each parameter set. By using the V of one oscillator, we obtained three pairs of inferred parameters. By using the V of the other oscillator, we obtained three additional pairs. Thus, we have six pairs of inferred values for each parameter set.

The results of the simultaneous inferences of noise intensity and coupling strength with methods I and II are shown in Fig. 2 , respectively. In Fig. 2, the inferred values approximately reproduce the actual values even though only one oscillator was observed. The error in the inference increases as the coupling or noise intensity is increased. In Fig. 2, the inferences by method II are obviously an improvement on the results of method I.

Fig. 2.

Simultaneous inferences of effective noise intensity aD and effective coupling strength obtained with (A) method I and (B) method II. Actual values are plotted as crosses. Inferred values are plotted as squares (for ) and circles (otherwise).

We emphasize that a naive use of the statistical values V and ζ will not yield successful inferences of noise intensity and coupling strength. The correlation between V1 and aD is shown in Fig. 3, and that between and is shown in Fig. 3. We found that their correlation coefficients were 0.96 and 0.70, respectively. In contrast, the correlation coefficient between the actual and inferred noise intensities (coupling strengths) for method II was 0.99 (0.99), as shown in Fig. 3 (Fig. 3). This fact indicates that our methods are superior to the naive use of V and ζ. In addition, the naive use provides only relative intensities, whereas our methods directly infer the absolute values of aD and .

Fig. 3.

Raw data on (A) period variance V1 versus effective noise intensity aD and (B) inverse time lag versus effective coupling strength . The lines in A and B are drawn by using the least-squares method. For comparison, the inferred values of (C) aD and (D) with method II are plotted versus the actual values of aD and , respectively. While our methods achieved precise inferences, V1 and did not.

Next, we demonstrated the inference method for a more realistic system. Specifically, we employed the coupled FitzHugh–Nagumo oscillators, described byfor and (1, 2). We set , , and . Each was the same as that in the inferences discussed above. This system shows limit cycle oscillation with a period when noise and coupling are absent.

The actual values of a and c for this system were obtained as follows: to calculate a, we numerically integrated the function Z representing the phase response to noise, and to calculate c, we observed the relaxation of the phase difference between two oscillators in a system with a fixed κ but without noise. The phase difference is expected to exponentially decrease by a factor of each period. We adopted the value of obtained from this relaxation as the effective coupling strength.

We prepared 16 parameter sets with and , where . We integrated Eq. using the Euler scheme with a time step of . In this simulation, the checkpoint threshold was fixed at , and the kth spike time was defined as the time at which v first passes through in the kth oscillation. We observed the spike timing for . The observed number of oscillations was about , corresponding to a day in the experiment on cardiac myocytes (1). Three realizations were simulated for each parameter set.

The results of simultaneous inferences for the FitzHugh–Nagumo model with methods I and II are shown in Fig. 4 and Fig. 1, respectively. Fig. 4 reveals that precise inferences were obtained with method I, except in the case of the greatest coupling strength, even though data on only one oscillator were employed. Fig. 1 reveals that better inferences were obtained with method II because the additional information raised the precision of each inference. The errors tend to be large when the coupling or noise intensity is large, as also observed in the case of the phase oscillator model.

Fig. 4.

Simultaneous inference of effective noise intensity aD and effective coupling strength for the FitzHugh–Nagumo model with method I. Actual and inferred values are plotted in the same manner as Fig. 1. The inferences were successful overall, although the errors became larger as the coupling strength was increased.

Discussion and Conclusion

We developed methods to infer coupling and noise intensities in two well-synchronized oscillators. Our methods successfully provide precise inference not only for the phase oscillator model but also for the FitzHugh–Nagumo model, which is a representative model of neural and cardiac oscillations. Here we discuss the robustness, limitations, and possible extensions of our methods.

One of the main factors that hamper the inference is the finiteness of sample data. In our numerical demonstrations, the number n of the observed spikes was or 106. This number would be feasible for the observation of cardiac myocytes as it approximately corresponds to the number of beats in a day. However, this number may be difficult to achieve in some systems; therefore, it is important to check the robustness of our methods against smaller sample sizes. In , we present the dependence of the inferred noise and coupling intensities on n in the FitzHugh–Nagumo model with two parameter sets. It indicates that although the absolute values of the inferred values are imprecise, relative relation between the two parameter sets can be inferred even for . There, we also present a convenient method to estimate the inference error from a single dataset by utilizing the bootstrap method (31).

The inference error is expected to be significant when the system is far from the in-phase state because our methods are based on the linear approximation around the in-phase state. This failure of the linear approximation typically arises when noise and/or oscillator inhomogeneity are strong. We numerically checked the robustness of method I against strong noise using the phase oscillator model with a fixed coupling strength; a detailed description has been provided in . The noise effect on synchronization is conveniently captured by the order parameter r for the synchronization (see for details). The actual and inferred intensities versus the observed order parameter r are shown in . As observed, the inference error is not significantly large for . With this noise strength, we obtained , at which the typical phase difference between oscillators is approximately ; the system is not significantly close to the in-phase state. Moreover, with such strong noise, phase slips occur occasionally; therefore, there is a small discrepancy between the time-average frequencies of the oscillators for a finite observation time. Despite the phase slip occurrence, method I is still adequately effective. We also investigated the effect of oscillator inhomogeneity in . We numerically confirmed that method I provides precise inference unless the difference between the oscillator frequencies is too large. These results indicate that method I is effective even when the synchronization is somewhat violated by the noise or inhomogeneity. It should be noted that method II uses the time lag between two oscillators; our definition of the time lag would not be appropriate for the case where the phase slips occur frequently. Therefore, we leave it as an open problem for investigation in future studies.

Strong coupling also hampers the inference due to the smallness of terms in Eqs. and . Because these terms rapidly decay as κ increases, the inference of coupling intensities become difficult for large κ. As is clearly observed in Fig. 2, the inference gets worse for higher κ values. By comparing Fig. 2 , we also notice that the error in Fig. 2 (thus method II) is smaller for high κ values. One reason for this is that method II avoids using .

For strong coupling and/or noise, the trajectory of each oscillator largely deviates from the limit cycle orbit of the unperturbed oscillator. Subsequently, the inference error is expected to increase owing to the following reasons. First, the phase oscillator models do not approximate the dynamics of limit cycle oscillators well (17). Our formulae are based on the phase oscillator model and, thus, cannot precisely predict the behavior of the limit cycle oscillators. Second, a spike time is defined in different manners between a phase oscillator and a limit cycle oscillator: a spike time in the phase oscillator is defined by a time passing a certain phase, whereas that of the limit cycle oscillator is defined as a time at which an oscillator passes a certain linear cross-section defined arbitrarily in the one-oscillator state space (32). According to the phase reduction theory (17), the linear cross-section in the state space generally differs from the isochron, which is defined as a curved cross-section comprising isophase points. The inference error involved with this discrepancy is expected to be larger as the trajectory is farther from the limit cycle orbit.

Finally, we discuss some extensions of our theory. Although our proposed methods are for two coupled oscillators, their extensions to the case of a population of oscillators with all-to-all coupling is straightforward, as briefly explained in . In addition, it is important to consider the presence of common noise (33, 34) because it ubiquitously exists as, e.g., environmental noise. This extension is also straightforward, as shown in . A further extension of our methods to the systems subjected to colored noise (35) is a significant challenge. Such an extension would contribute to the inference of coupled chaotic oscillators (33, 36–38) as chaotic oscillators can be regarded as a periodic oscillator subjected to colored noise. These extensions are expected to enhance the applicability of our methods.