A Weighted Spike-Triggered Average of a Fluctuating Stimulus Yielding the Phase Response Curve

We demonstrate that the phase response curve (PRC) can be reconstructed using a weighted spike-triggered average of an injected fluctuating input. The key idea is to choose the weight to be proportional to the magnitude of the fluctuation of the oscillatory period. Particularly, when a neuron exhibits random switching behavior between two bursting modes, two corresponding PRCs can be simultaneously reconstructed, even from the data of a single trial. This method offers an efficient alternative to the experimental investigation of oscillatory systems, without the need for detailed modeling.

We demonstrate that the phase response curve (PRC) can be reconstructed using a weighted spike-triggered average of an injected fluctuating input. The key idea is to choose the weight to be proportional to the magnitude of the fluctuation of the oscillatory period. Particularly, when a neuron exhibits random switching behavior between two bursting modes, two corresponding PRCs can be simultaneously reconstructed, even from the data of a single trial. This method offers an efficient alternative to the experimental investigation of oscillatory systems, without the need for detailed modeling. Synchronization phenomena in networks consisting of interacting nonlinear dynamical elements exhibiting limit-cycle oscillations have been the subject of intensive study in physical, biological and social systems [1]. In general, such a system is described by ordinary differential equations of the form dX dt = F (X), where X denotes a multi-dimensional state of the system. Owing to nonlinearity and complicated interactions, however, we often encounter difficulties when attempting to analyze synchronization properties of such systems (see arrow A in Fig.1). One of the most successful approaches for dealing with such difficulties is the phase reduction method (arrow B), in which the evolution of each active element is described by only one degree of freedom, the phase [2]. With this method, the dynamics of N coupled oscillatory systems can generally be described by a set of equations of the form dφi dt = ω i + N j=1 Γ ij (φ i − φ j ), where ω i is the natural frequency without coupling. The phase of the i-th system representing the timing of its limit-cycle oscillation is denoted by φ i . This description is valid if the perturbation due to the interactions is so weak that its effect only changes the phase asymptotically. The coupling function Γ ij (φ) can then be derived from the original dynamical system as represent the interaction exerted by the j-th system on the i-th system and the limitcycle solution for the uncoupled i-th system, respectively [3]. Here, Z i represents the phase response curve (PRC) of the i-th system. Therefore, if we can obtain the PRC, then we can predict the synchronization properties of the coupled system for any weak interaction (arrows C).
It is well known that the PRC can be calculated from the original model equations with the adjoint method [4]. Owing to technical difficulties limiting the applicability of experimental studies, however, it is often difficult to construct system-specific detailed models including the 0 0 essential properties for synchronization (arrow D). For instance, if we wish to determine how a given neuromodulator causes the synchronization properties of a coupled neuronal system to change, we cannot evaluate all the effects of the neuromodulator on the dynamics of the various ion channels. As an alternative approach, there have been some recent attempts to directly measure the PRC experimentally (arrow E). The most important advantage of this approach is that, once the exact PRC is known, the dynamics of a set of weakly coupled oscillators are fully determined, whether or not the original detailed dynamics are known. In the example considered above, therefore, all the effects of the neuromodulator are reflected in the measured PRC, and with this, the dynamical behaviors can be correctly predicted. The standard method for obtaining PRCs through direct experimental measurements is to measure the phase shift induced by stimulating the system with a brief pulse-like perturbation at various phases of the period [5]. Unfortunately, in such a study we often face the difficulty that the information from which we can derive the PRC is lost in the noise arising from the uncontrollable nature of the environment, though there are several statistical fitting methods which have been applied to efficiently extract the PRC from noisy data [6].

Real system
An alternative approach for obtaining the PRC which we consider here, is to design a different experimental procedure for the measurement, instead of the standard method described above. In connection to this, G. B. Ermentrout et al. have recently shown that the spiketriggered average (STA) of the injected noisy current is proportional to the derivative of the PRC under suitable conditions [7]. This result inspired us to develop a more useful, novel method for experimental protocols.
In this paper, we demonstrate that the PRC can be directly reconstructed using an appropriately weighted spike-triggered average of the injected fluctuating inputs.
The key idea in this method is to choose the weight to be proportional to the magnitude of the fluctuation of the oscillatory period. In the conventional method employing pulse-like perturbations, the spontaneous fluctuation of the period tends to make the estimation of the PRC quite rough. Interestingly, the disadvantage of the fluctuating period in the case of the conventional method becomes an advantage in our method. Although the proposed method can be applied to various real oscillatory systems, such as Belousov-Zhabotinsky reactions, for simplicity, we consider two neuronal systems to explain the experimental procedure in the following.
Experimental procedure -We consider the situation in which a neuron with some constant bias current exhibits spikes regularly with a period T . When an additional fluctuating current with zero mean is injected into this neuron, the inter-spike interval (ISI) generally fluctuates about the average value T , as shown in Fig.2(a). As the first step of the experimental procedure, we record the stimulus current I(t) and all the spike times over the entire spike train. Next, we divide the stimulus I(t) into individual stimuli I i (t) between successive spikes, whose interval duration is denoted by τ i (i = 1, . . . , N ; N denotes the number of the recorded spikes). Then ( Fig.2(b)), the slightly different duration times of the stimulus currents I i (t) are rescaled to the uniform period T as Finally, the weighted spike-triggered average (WSTA) of the rescaled stimulus current I i (t i ) is defined as where ∆ i = T −τi τi , and the angular brackets represent the average over all spikes, i.e. the index i (Fig.2(c)). As shown below, we can prove that the above-defined WSTA is proportional to the PRC, provided that the fluctuations in the stimulus current are not too strong.
Theory -When a neuron fires almost periodically under the influence of a zero-mean fluctuating stimulus current I(t), the corresponding phase dynamics can generally be described by the equation in which ω = 2π T and Z(φ) represents the PRC with respect to the input current. Let us consider the effect of the i-th stimulus, I i (t) ≡ µξ i (t), where the parameter µ is introduced to account for the effect of the magnitude of the fluctuating current. Hereafter, we omit the subscript i of the rescaled timet i for simplicity. In addition, we assume that ξ i (t) is a stationary stochastic process with unit strength. The auto-correlation of ξ i (t) is then defined as C(t) ≡ ξ i (t)ξ i (0) . Then, substituting I(t) = µξ i (t) into Eq.(3), and integrating from t = 0 to τ i , we obtain ∆ i = µ 2π T 0 Z(φ i (s))ξ i (s)ds, where Eq.(1) has been used. Next, substituting this into Eq.(2), we find We assume that we can expand the time development of the phase with respect to µ as φ i (t) = ωt + µφ i (t) + · · · . Using this and the Taylor expansion in Eq.(4), we get Assuming that the timescale on which the autocorrelation of the stimulus current decays is much shorter than the period T , we can approximate C(t) as δ(t), where δ is the Dirac delta function. Using the fact that µ is then given by I i (t)I i (0) = µ 2 δ(t), we finally obtain where φ = ωt. We have thus found that the PRC can be reconstructed from the spike-triggered average of the stimulus current using a weight proportional to the normalized difference between the ISI and the average interval.
Numerical examples -To confirm the validity of our proposed method, we compare the PRCs obtained from the WSTAs and the PRCs calculated using the adjoint method. In Fig.3, we first show typical results in the case of the Hodgkin-Huxley model [8]. Figure 3(a) displays a typical voltage trajectory in response to the fluctuating current. Figure 3(b) presents comparisons between the true PRC (black trace) and the PRC estimated from the WSTA (red trace) for three different measurement durations. We see that as the measurement duration increases, the WSTA-estimated PRC converges to the true one. This result suggests that the PRCs of real systems can be accurately estimated within practically reasonable recording times experimentally.
We next consider the more complicated situation in which a bursting neuron exhibits singlet or doublet firing randomly, under the influence of the injected fluctuating current, as shown in Fig.4(a). For the bursting neuron, here we adopt the chattering neuron model, which exhibits an increasing number of intra-burst spikes, such as from singlet to doublet, when the injected current is increased [9]. In the case considered in Fig.4, the injected current fluctuates about the mean level corresponding to the transition between singlet and doublet firings. Using the WSTA with an additional procedure, we can simultaneously reconstruct two PRCs corresponding to singlet and doublet states from a single trial. The key idea here is to separately calculate two conditional WSTAs, depending on whether a singlet or doublet firing occurs. In other words, the PRC for the singlet (doublet) firing can be reconstructed by calculating the WSTA over only the injected currents generating the singlet (doublet) firing. In Fig.4(b), closer investigation reveals that the multi-modal distribution of the ISI can be regarded as a mixture of four different unimodal distributions, which are specified by two successive firing modes. From each unimodal distribution, the average ISI is separately com- puted for calculation of the ∆ i in Eq.(4). Finally, two rescaled sets of data with the same final firing mode are used to obtain the two PRCs for the singlet and doublet modes. Figure 4(c) demonstrates that both of the PRCs obtained in this way are in reasonably good agreement with the exact ones. We emphasize that two different PRCs can be reconstructed from a single-trial data set, even when the naive conventional method employing a pulse-like perturbation is not practical. Noting the sharp difference between these two WSTA-estimated PRCs, we predict a transition between anti-synchronous and synchronous firing for a two-neuron system with excitatory synaptic couplings, as shown in a previous study [10] (Fig.4(d)).
Discussion -In actual applications to real systems, we have to choose an appropriate magnitude of the fluctuation of the injected current. Since some unpredictable and uncontrollable noise is inevitable in real systems, we examine the discrepancy between the estimated PRC and the true one in the presence of unknown background noise, as summarized in Fig. 5. We find that for two types of Morris-Lecar models near the threshold for firing [12], the WSTA estimates the PRC most precisely in the case that an intermediate magnitude of the fluctuation of the injected current is chosen. This is because the signal of the PRC becomes lost in the background noise if the injected fluctuations are too weak, while if the fluctuations are too strong, the approximation of the phase description becomes poor. Figures 5(c) and (d) suggest that the same results hold for all types of bifurcations and smoothing algorithms considered.
We now give some comments on three relevant studies. First, when the fluctuation of the ISIs vanishes, all the c1 weights of the WSTA are the same. This is equivalent to the situation in which the STA yields Z ′ (φ) instead of Z(φ) [7]. This superficial inconsistency can be resolved by considering the fact that the difference between the WSTA and the STA converges to zero as µ 2 . Second, c1 recent studies have pointed out that a correction term appears in the phase description (3) if the limit-cycle oscillation is perturbed by the noise with a correlation time shorter than the relaxation time of the limit-cycle attractor [13]. Under real conditions, a fluctuating noisy c2 signal can be regarded as a smooth signal in the limit of a small time scale. Therefore, Eq. (3) and our result are practically valid. Third, in another study [14], c1 the collective PRC to a macroscopic external force for globally coupld oscillators is investigated. In principle, c2 our method can be applied to measure such a collective PRC. Closer examinations of the above-mentioned issues c3 are beyond the scope of this paper, but they are of great importance and should be studied further.
In conclusion, we demonstrated that the PRC can be reconstructed from our proposed WSTA, in which the weight is proportional to the magnitude of the fluctuation of the period. In this study, considering limit-cycle oscillations observed ubiquitously in nonlinear dissipative systems far from equilibrium, we found a theoretical relation between the fluctuations in the system and its response to an external force. This might provide useful insight for further studies as the development of fluctuation dissipation theorem. In comparison with the standard method employing pulse-like perturbations, furthermore, the proposed WSTA method is experimentally reliable, fast and wide applicable, as shown by the fact that two different PRCs can be reproduced even from a single-trial data set for two-mode bursting neurons. We believe that this method will contribute to elucidating the nature of real dynamical system in a broad range of contexts.