Gaussian Process Phase Interpolation for estimating the asymptotic phase of a limit cycle oscillator from time series data
Taichi Yamamoto Hiroya Nakao Ryota Kobayashi

















































Tokyo1]Taichi Yamamoto

[Tokyo1]organization=Graduate School of Frontier Sciences, The University of Tokyo,   

city=Kashiwa, postcode=277-8561, 
country=Japan

TokyoTec1,TokyoTec2]Hiroya Nakao
[TokyoTec1]organization=Department of Systems and Control Engineering, Tokyo Institute of Technology,  
            
            city=Tokyo,  postcode=152-8552,  
            country=Japan  
[TokyoTec2]organization=Research Center for Autonomous Systems Materialogy, Tokyo Institute of Technology,  
            
            city=Yokohama,  postcode=226-8501,  
            country=Japan  

Tokyo1,Tokyo2]Ryota Kobayashi
r-koba@k.u-tokyo.ac.jp

[Tokyo2]organization=Mathematics and Informatics Center, The University of Tokyo,   

city=Tokyo, postcode=113-8656, 
country=Japan            

§ ABSTRACT 


Rhythmic activity commonly observed in biological systems, occurring from the cellular level to the organismic level, is typically modeled as limit cycle oscillators. 
The phase reduction theory serves as a useful analytical framework for elucidating the synchronization mechanism of these oscillators. Essentially, this theory describes the dynamics of a multi-dimensional nonlinear oscillator using a single variable phase model. In order to understand and control the rhythmic phenomena in the real world, it is crucial to estimate the asymptotic phase from the observed data. 

In this study, we propose a new method, Gaussian Process Phase Interpolation (GPPI), for estimating the asymptotic phase from time series data. The GPPI method first evaluates the asymptotic phase on the limit cycle and subsequently estimates the asymptotic phase outside the limit cycle employing Gaussian process regression. Thanks to the high expressive power of Gaussian processes, the GPPI is capable of capturing a variety of functions. Notably, the GPPI is easily applicable even when the dimension of the system  increases. The performance of the GPPI is tested by using simulation data from the Stuart-Landau oscillator and the Hodgkin-Huxley oscillator. The results demonstrate that the GPPI can accurately estimate the asymptotic phase even in the presence of high observation noise and strong nonlinearity. Additionally, the GPPI is demonstrated as an effective tool for data-driven phase control of a Hodgkin-Huxley oscillator. 
Thus, the proposed GPPI will facilitate the data-driven modeling of the limit cycle oscillators. 















Synchronization Limit cycle oscillators Phase reduction Machine learning   Gaussian process regression

















§ INTRODUCTION




Rhythmic activity is a ubiquitous phenomenon observed in a wide range of biological systems. 
Examples include cortical networks in the brain <cit.>, circadian rhythms in mammals <cit.>, the human heart and respiratory system <cit.>, and animal gait <cit.>. Most mathematical models of rhythmic activity are based on limit cycle oscillators <cit.>. 

Phase reduction theory <cit.> is a valuable tool for analyzing the synchronization of limit cycle oscillators. This theory represents the state of a multi-dimensional limit cycle oscillator using a single variable, called the phase, and describes their dynamics in a reduced phase model. Theoretical studies based on the phase model have elucidated the key factors underlying synchronization phenomena, including the periodic external force, the coupling between oscillators, and the common inputs <cit.>. 


A fundamental challenge in the field of complex systems is to identify the mechanisms by which a real-world system achieves synchronization <cit.>. 
While these theoretical studies offer explanations for synchronization phenomena, they require the precise knowledge of the mathematical models. 
Consequently, previous studies have developed methods for identifying the phase model from real data (for a review, see <cit.>).

The most common approach for identifying the phase model from data is to estimate the phase sensitivity function (also known as the  infinitesimal phase response curve), which characterizes the linear response property of the oscillator phase <cit.>. 
However, when a limit cycle oscillator is subjected to strong perturbations, such as a strong impulse input, the approximation using the phase sensitivity function may be inaccurate. 
Therefore, it is essential to obtain the asymptotic phase distant from the limit cycle in order to understand the synchronization property and to control the phase of the limit cycle oscillator subjected to strong perturbations. Recently, <cit.> proposed a method to obtain the asymptotic phase in a data-driven manner. However, it is still challenging to estimate the asymptotic phase of high-dimensional systems. 


In this study, we propose the Gaussian Process Phase Interpolation (GPPI) method, which estimates the asymptotic phase from time series obtained from a limit cycle oscillator without assuming the mathematical model of the system. The GPPI method can be readily applied to systems with more than three dimensions, as the same algorithm is used regardless of the dimensionality of the dynamical system. 
To validate the GPPI, we apply it to two limit cycle oscillators: the Stuart-Landau oscillator and the Hodgkin-Huxley oscillator. 
Furthermore, we apply the GPPI to control the phase of the Hodgkin-Huxley oscillator by impulse inputs, which demonstrates the effectiveness of the GPPI for data-driven control. 

This paper is organized as follows: Section 2 outlines the phase reduction theory and the asymptotic phase. Section 3 describes the proposed method, the GPPI, for estimating the asymptotic phase. In Section 4, the GPPI is validated using synthetic data obtained from the two limit cycle oscillators. In Section 5, we demonstrate data-driven control of a Hodgkin-Huxley oscillator by impulse input using the estimated asymptotic phase. Finally, Section 6 summarizes the results and discusses the conclusions. 














§ PHASE REDUCTION OF LIMIT-CYCLE OSCILLATORS





We consider a limit-cycle oscillator described by

    d/dt(t) = 𝐅((t)),

where (t)∈^d is the state vector of the system at time t and 𝐅:^d→^d is a smooth vector field representing dynamics. 
We assume that the system has a stable limit cycle orbit (t) with the natural period and frequency, T and ω= 2π/T, respectively, which satisfies the condition: (t+T) = (t). 

For any point in the basin of the limit cycle, an asymptotic phase can be determined <cit.>. The basin ℬ is defined as the set of initial conditions that converge to the limit cycle. 
In the following, the asymptotic phase is simply written as the `phase'. 

First, we assign the phase θ∈ [0,2π) on the limit cycle, where the phase θ= 0 and 2 π are considered to be identical. 
We choose a point on the limit cycle, which is set to be the phase origin, θ=0, and we define the phase to increase at a constant rate ω, i.e., θ(t) = ω t (mod 2π). 
As a result, the phase of the state (t) on the limit cycle follows 

    d/dtΘ((t)) = ω,
 
where the phase function Θ() gives the phase θ of the state  on the limit cycle. 

Moreover, the phase function can be extended so that Eq. (<ref>) holds for any orbit (t) in the basin of the limit cycle ℬ:

    d/dtΘ((t)) = ω.

Here, we extended the domain of the phase function as Θ: ℬ→ [0,2π).  
The above procedure reduces the d-dimensional differential equation (<ref>) to a simple one-dimensional equation (<ref>) for any (t) ∈ℬ. By integrating Eq.(<ref>), we obtain 

    Θ((t+τ)) - Θ((t)) = ωτ   
        ( mod 2π).



We introduce the phase response function (PRF, also known as the phase response or resetting curve, PRC) <cit.>, which is essential for the analysis and control of synchronization phenomena. 
The PRF is a function that describes the effect of an impulse perturbation applied to an oscillator in phase θ and is defined as

    g(θ,   ς  ) = Θ((θ)+ ς ) - Θ((θ)) = Θ((θ)+ ς ) - θ,


where ς∈^d denotes the intensity and direction of the impulse applied to the system (<ref>) <cit.>, and 
𝒳_0(θ) denotes the point with phase θ on the limit cycle. 
If the impulse intensity ς is sufficiently small, the PRF can be approximated as 

    g(θ, ς ) = 𝐙(θ)·ς +O(ς^2)

where

    𝐙(θ) = ∇Θ((θ))

is called the phase sensitivity function (PSF, also known as the infinitesimal phase resetting curve, iPRC). The PSF characterizes the linear response to a given weak perturbation for the state of phase θ.







§ PROPOSED METHOD


We propose a method, Gaussian Process Phase Interpolation (GPPI), for estimating the phase function Θ() from time series data of multiple orbits converging to the limit cycle. 
The main idea is to obtain the phase function from the time series data using the phase equation (<ref>) and interpolate them by Gaussian process regression. 
The procedure can be divided into the following three steps (Fig. <ref>):   

    
    
  1. Determine the phase function on the limit cycle. 
    
  2. Obtain the phase function outside the limit cycle. 
    
  3. Interpolate the phase function using Gaussian process regression.

  





 §.§ Step 1: Determine the phase function on the limit cycle
  


In the first step, we obtain the value of the phase function on the limit cycle. 
Let us assume that a sufficiently long time series, _0(t), is obtained from the limit cycle oscillator (<ref>) and that at time t=0, the orbit has converged sufficiently to the limit cycle.  

We first estimate the natural frequency ω. 

Let us consider a Poincaré section that intersects the limit cycle once.  
The period and frequency of the limit cycle can be estimated using the following formulae: T̂=(s_n+1-s_1)/n and by = 2π / T̂, where s_1 and s_n+1 are the times of the first and (n+1)-th passage through the Poincaré section, respectively. 

In order to determine the phase, a point _0(t_0) (t_0 ≥ 0) on the time series is taken as the origin of the phase, i.e., Θ(_0(t_0))=0. Subsequently, the phase on the limit cycle can be obtained using the phase equation (<ref>): Θ(_0(t)) = (t-t_0) for any t ≥ 0. 




 §.§ Step 2: Obtain the phase function outside the limit cycle
 



In the second step, we obtain the phase function on each data point outside the limit cycle. 
Let us assume that we have collected n time series, represented by {_1(t), ⋯, _n(t) }, which converge to the limit cycle from its basin. 
It should be noted that such time series can be collected, for example, by applying perturbations to a limit cycle oscillator and observing the convergence to the limit cycle. 


We describe the procedure for obtaining the phase function on a data point that is outside the limit cycle. Let us consider a time series _i(t) whose last point _i(t_ end) is close enough to the limit cycle. 
We can calculate the phase of the last point, Θ(_i(t_ end)), by applying 

the linear interpolation of the data points on the limit cycle, Θ(_0(t)), obtained in Step 1. 
Then, the phase at any time t on the time series can be calculated from the last point using the phase equation (<ref>) as

    Θ(_i(t)) = Θ(_i(t_ end)) + ( t- t_ end ).


In this way, the phase of all points on the time series _i(t) are obtained. Similarly, the phase of each data point can be obtained by repeating this procedure for all time series (i= 1, 2, ⋯, n). 
In this study, the entire time series is discarded if the last point is not sufficiently close to the limit cycle. 





 §.§ Step 3: Interpolate the phase function using Gaussian process regression

  


In the third step, we estimate the global phase function Θ() by employing Gaussian process regression. 
The proposed method does not directly estimate the phase function, Θ(), but rather estimates two real-valued functions 

    s()    := sinΘ(), 
    
        c()    := cosΘ(),

and obtains the phase function from these functions. 

First, we generate two datasets {(_i,s(_i))} and {(_i,c(_i))} (i=1, 2, ⋯, N) from the data {(_i,Θ(_i))} obtained in step 2 (Sec. <ref>). 

Next, we estimate the two functions (s(), c()) from the datasets using Gaussian process regression. 
Finally, the phase function is estimated as 

    Θ̂()=arctan(ŝ()/ĉ()),

where ŝ() and ĉ() are the estimates of Gaussian process regression. 


Gaussian process regression is a non-parametric Bayesian approach for estimating a real-valued function f:ℝ^d→ℝ <cit.>. 
Gaussian process regression has two main advantages: 
1) it captures a wide variety of functions by leveraging a theoretically infinite number of parameters, and 
2) it evaluates the uncertainty of the function by predicting the distribution of the function from the data.

From the data { (_1,y_1), ⋯, (_N,y_N) }, Gaussian process regression predicts the mean m(_*) and variance v(_*) of the function f(_*) for an unknown input _* as follows: 



    m(_*)    = _̨*^⊤(K+σ^2I)^-1, 
    
        v(_*)    = k(_*,_*) - _̨*^⊤(K+σ^2I)^-1_̨*,

where 

    K   ∈ℝ^N× N,  
        K_ij = k(_i,_j), 
    _̨*   ∈ℝ^N,  _̨* = ( k(_1,_*), k(_2,_*), ⋯ , k(_N,_*) )^⊤, 
    
           ∈ℝ^N,   = (y_1, y_2, ⋯ , y_N)^⊤,

I ∈ℝ^N× N is the identity matrix, 
and σ^2 is the variance of the observation noise, respectively. 

The kernel function k(_i,_j) represents the correlation between the outputs (f(_i) and f(_j)). 
Here, we adopt the Matern-5/2 kernel for the kernel function: 


    k(_i,_j)    = σ_f^2
        (1+√(3)r/σ_l)
        exp(-√(3)r/σ_l),

where r= √((_i-_j)^⊤(_i-_j)) is the distance between the data points, and σ_f and σ_l are hyperparameters that determine the height and width of the kernel function, respectively. 
The hyperparameters are determined by maximizing the log-likelihood function given by 

    logp(|_1,...,_N)
    
        = -1/2^⊤(K+σ^2I)^-1
        -1/2log(K+σ^2I)-d/2log2π.

It should be noted that the noise variance is set to σ^2=0.01. 
In this study, we used Gaussian process regression algorithm  implemented in MATLAB Statistical and Machine Learning Toolbox <cit.>. 







§ TEST OF THE PROPOSED METHOD



In this section, we test the validity of the proposed method, GPPI, by applying it to two limit cycle oscillators, the Stuart-Landau oscillator and the Hodgkin-Huxley oscillator.







 §.§ Evaluation methods



We test the proposed method by calculating the true phase function through exhaustive numerical simulation based on the mathematical model (1) and comparing it with the estimation results. 
The true phase function is calculated as follows: First, we obtain the phase function on the limit cycle in the same way as in the proposed method (Step 1, Sec. <ref>). Then, for each state  in a region in the basin, we simulate the dynamical system (<ref>) from  until the orbit converges to the limit cycle. Finally, the phase function Θ() is obtained using the phase equation (<ref>) in the same way as in the proposed method (Step 2, Sec.<ref>). 
In this study, we use a two-dimensional grid in the state space, and calculate the phase function at all grid points. It is noteworthy that we use noiseless simulation to calculate the true phase function.



Moreover, we test the proposed method by comparing the normalized phase response function (nPRF) with its estimate.  
The nPRF for the x direction is defined as G_x, κ(θ) := g(θ, κ𝐞_x) / κ, where 𝐞_x is the unit vector in the x direction and κ is the intensity of the impulse. 
The nPRF can be calculated from the estimate of the phase function Θ̂ as follows: 

    Ĝ_x,κ(θ) :=
        Θ̂((θ)+ κ𝐞_x) - Θ((θ))/κ,

where (θ) denotes a state of the phase θ on the limit cycle.   


The estimation accuracy of the nPRFs is evaluated using the coefficient of determination defined as follows:

    R_Y^2 = 1 - ∑_i(Y(θ_i)-Ŷ(θ_i))^2/∑_i(Y(θ_i)-Y̅)^2,

where Y and Ŷ are the true and estimated values of the nPRF, respectively, and Y̅ is the average of Y over all the phase θ_i. The closer the coefficient of determination is to 1, the higher the estimation accuracy is. 







 §.§ Stuart-Landau oscillator















We evaluate the proposed method using the Stuart-Landau (SL) oscillator and compare it with an existing method, which we refer to as the Derivative Phase Regression method (DPR method) <cit.>. We also examine the robustness of these methods against observation noise. 


The SL oscillator is described as follows:		

    ẋ   = x - α y - (x - β y)(x^2+y^2),		
    ẏ   = α x + y - (β x + y)(x^2+y^2),
 
where x, y are the state variables and α=2, β=1 are parameters. The numerical simulation was conducted using the fourth-order Runge-Kutta method with a time step of Δ t=0.005. 


We generate synthetic data for validation as follows: an orbit is simulated based on Eqs. (<ref>) and (<ref>) for the time length of 2.5 from an initial state randomly sampled from the uniform distribution of [-1.6, 1.6]^2. It should be noted that the period of the limit cycle is 2 π. The time series data { x(t), y(t) } is collected if the final state is sufficiently close to the limit cycle. This procedure is repeated until we obtain n= 100 time series. 
Finally, observation noise is added to the time series for testing the robustness against noise. 
The observation noise is given by independent Gaussian random variables of mean 0 and standard deviation of η=0.005. 




We first estimate the natural frequency (Step1, Sec. <ref>). To reduce the impact of observation noise, a moving average with a window width of 0.07 was calculated from the time series data.
The natural frequency is estimated as = 1.000, which is in agreement with the theoretical value of ω=1.0 derived from the model. The estimated frequency = 1.000 is also used in the DPR method. 

Next, we estimate the phase function Θ from the time series using the proposed method, the GPPI, as well as the DPR. 
To evaluate the performance of the proposed method for a small dataset, 1,000 data points are sampled with a coarse interval Δ t = 0.25, and are used as a training data set for Gaussian process regression. 
In contrast, the DPR uses all simulation data, 50,000 points with sampling interval of Δ t = 0.005. The other parameters of the DPR were the same as those used in <cit.> (see Appendix A for details).



Figure <ref> compares the true phase function with its estimate by the proposed GPPI method and the DPR method from the time series. 
While the both methods accurately estimate the phase function near the limit cycle, they cannot estimate the phase function around the origin (x,y)=(0,0) due to the lack of data. 

In particular, in the boundary region, the GPPI provides a more accurate estimate of the phase function than the DPR (see Fig. <ref>a in Appendix B).  
One potential explanation for the decline in the estimation performance of the DPR is the use of a polynomial function, which tends to diverge the absolute value of the function in the boundary region. Instead of the polynomial function, the GPPI employs Gaussian process regression for estimating the phase function. 
Overall, the proposed GPPI method provides an accurate estimation of the phase function in the global region with smaller data compared to the DPR method (1,000 vs 50,000 data points). 


We also compared the normalized phase response function (nPRF) calculated from the estimated phase functions with the true values. Figure <ref>a shows the true and estimated nPRFs obtained by applying impulses in four directions (+x,-x,+y,-y) at an intensity of 0.2.
Both methods yield accurate estimates of the phase response. However, the present GPPI method outperforms to the DPR method, with a higher coefficient of determination averaged over the four directions for the GPPI and the DPR (0.998 and 0.952). 
It should be noted that the accuracy of the DPR is worse than the results in <cit.> (Figure 3 and Table 1). This discrepancy may be attributed to the small number of time series used for the estimation. 



Furthermore, we examine the robustness of the methods against observation noise by generating new data sets with increasing noise levels, specifically, η=0.01, 0.05. The estimation procedure was repeated for both methods using the new data, including the estimation of the natural frequency. 

Figures <ref>b and c show the estimation results for nPRFs. As the noise strength increases, the estimation error of the DPR method also increases. The DPR is unable to estimate the phase response from noisy time series (η=0.05). 
In contrast, the GPPI is capable of estimating the phase response even in the presence of the substantial observation noise (η=0.05).  
Table <ref> shows a comparison of the performance of the GPPI with that of the DPR based on the average of the  coefficient of determination R_G^2. 
These results demonstrate that the proposed GPPI method can estimate the phase function with less data than the DPR method and it is also robust against noise. 











 §.§ Hodgkin-Huxley oscillator




As a second example, we test the proposed method using the Hodgkin-Huxley model (HH model), which is a four-dimensional dynamical system. The HH model has a higher dimension than the SL oscillator and also exhibits stronger nonlinearity than the SL oscillator. 



The HH model is one of the most important neuron models in neuroscience, which has been used to describe neural systems mathematically <cit.>. This model is written as a differential equation of four variables, namely, V, h, m, and n, as follows <cit.>:




    V̇ =   1/C( -Na(V-Na)m^3h-K(V-K)n^4 
        -L(V-L) + I_b ),  
    ṁ =   α_m(V)(1-m) - β_m(V)m,	  
    ḣ =   α_h(V)(1-h) - β_h(V)h,	  
    ṅ =   α_n(V)(1-n) - β_n(V)n,

where the functions α_m,h,n, β_m,h,n are written as

    α_m(V)    = 0.1(V+40)/1-e^-(V+40)/10,
           β_m(V)    = 4e^-(V+65)/18,
    α_h(V)    = 0.07e^-(V+65)/20,
           β_h(V)    = 1/1+e^-(V+35)/10,
    α_n(V)    = 0.01(V+55)/1-e^-(V+55)/10,
           β_n(V)    = 0.125e^-(V+65)/80,

and the parameters are set as follows:  C= 1 [μ F/cm^2],  Na  = 120 [mS/cm^2],   K= 36 [mS/cm^2],   L= 0.3 [mS/cm^2], Na = 50 [mV],  K = -77 [mV], L = -54.4 [mV],  and I_b= 10 [μ A/cm^2].  

The numerical simulation was conducted using the fourth-order Runge-Kutta (RK4) method with a time step Δ t=0.01 [ms].





 



A limit cycle that corresponds to periodic firing can be obtained from the HH model by setting the input current I_b above a certain threshold <cit.>. In this study, we focus on the oscillatory regime, where the HH model has a limit cycle.  
Figure <ref> illustrates a single period of the limit cycle obtained from the HH model. 
Following the literature <cit.>, the origin of the phase is defined as the peak of the voltage V. 
The period and the natural frequency are estimated from a time series of 1,000 periods, which results in T̂=14.6 [ms] and =0. 429 [rad/ms]. It should be noted that these values are also used in the calculation of the true phase function, as we do not consider the observation noise in this example. 



Here, we consider two strategies for collecting the time series data used to estimate the phase function. 

The first strategy is to simulate the system from the initial state, which is randomly sampled from a uniform distribution of a hyperrectangle, as in the case of the SL oscillator. For the HH oscillator, we employ a hyperrectangle formed by V ∈ [-100, 50], m ∈ [0,1], h ∈ [0, 0.6], n ∈ [0.3, 0.8]. 
We refer to this strategy as uniform sampling. The uniform sampling strategy is able to collect data over a global region of the state space. However, as the dimension of the state space increases, the data points become increasingly sparse, thereby rendering the estimation of the phase function more difficult. 

The second strategy is to simulate the system from the initial state generated by adding a small perturbation to the state on the limit cycle. For the perturbation of the HH oscillator, we employ a uniform random number in the range [-σ_x/10, σ_x/10], where σ_x is the standard deviation of the variable x ∈{V, m, h, n } along the limit cycle. 
We refer to this strategy as vicinity sampling. While the vicinity sampling is an effective strategy for accurately estimating the phase function of regions in proximity to the limit cycle, it is difficult to apply this strategy for estimating the phase function of regions distant from the limit cycle. 

We generate the orbits of time length 100 [ms] from 100 initial states obtained for each sampling strategy by simulating the HH model. Subsequently, the time series of each state variable is collected. It should be noted that, in uniform sampling, several time series that do not converge to the limit cycle are excluded from the analysis. 




The number of data points thus prepared is approximately 10^6. 
Given the high computational cost associated with Gaussian process regression, it is necessary to reduce the number of data points to a range of 10^3 to 10^4. Accordingly, the training data used for Gaussian process regression is extracted from the data set in the following manner. 

First, we extract the initial part of the time series, which is relatively distant from the limit cycle. 
Specifically, the initial 20 [ms] of the time series with a sampling interval of 1 [ms] are extracted in the uniform sampling, and the initial  6 [ms] of the time series with a sampling interval of 0.3 [ms] are extracted in the vicinity sampling.  

Furthermore, due to the rapid change in the state variables before and after a spike, the data density in the region corresponding to the spike tends to be relatively low. To ensure sufficient data collection in the spike region, we employ a smaller sampling interval. Specifically, the sampling interval is decreased to 1/10 during the spike period, t_f-4 ≤ t ≤ t_f+3 [ns], where t_f is the spike time defined as the time when the voltage V exceeds 0 [mV] <cit.>. 

Finally, 100 points from the limit cycle are added to the training data to supplement the phase information. 














Figure <ref> compares the true phase function obtained from the HH model with its estimate by the proposed GPPI method using the two sampling strategies: uniform and vicinity sampling. The limit cycle of the HH model approximately lies on a surface that can be written as follows <cit.>:

    m    = m_∞(V) := α_m(V) / {α_m(V) + β_m(V)},
    
        h    = 0.8 - n.
 
We therefore compare the phase function in the region V∈[-100,50], n∈[0.3, 0.8] on this surface. 
Although the data obtained from the vicinity sampling is distributed over a small region in proximity to the limit cycle, this strategy still worked well for estimating the phase functions of a more global region (see Fig. <ref>b in Appendix B). 
It should be noted that neither method is able to estimate the phase function in the vicinity of the equilibrium (V∼-60 [mV], n∼0.5), where the phase varies in a complex way <cit.>.






Next, we examine the estimation performance of normalized phase response functions (nPRFs). Figure <ref> compares the estimate of the nPRF from the HH oscillators, which shows that the vicinity sampling provides better results than the uniform sampling. 
Additionally, the difference between the nPRF for the positive and negative impulses reflects the nonlinearity of the HH oscillator. While the HH oscillator has a strong nonlinearity in the n direction, this nonlinearity can be captured by the vicinity sampling strategy. It should be noted that the intensity of the impulse used to calculate the nPRF is sufficiently low to fall within the range of the neighborhood sampling. 
Furthermore, we evaluate the estimation accuracy of the nPRF using the average of the coefficients of determination (Table <ref>). 
As can be seen in Fig. <ref>, the coefficients of determination of the vicinity sampling are much higher than those of the uniform sampling, except for the direction in h; both sampling methods provide accurate estimates of the nPRF in the h direction.  






§ DATA-DRIVEN PHASE CONTROL USING IMPULSE INPUT



As an application of the estimated phase function, we examine the control of the phase of a limit cycle oscillator using an impulse input. 
As in the previous section, we consider the Hodgkin-Huxley (HH) oscillator, but we consider a scenario that is more feasible in the experiments. 
Let us suppose that a time series (V(t),m(t),h(t),n(t)) can be measured from the HH oscillator that is stimulated by an impulse input with the direction of +V. 
Our objective is to control the phase of the oscillator from the time series without assuming any knowledge of the mathematical model (Eqs. <ref>, <ref>, <ref>, and <ref>). 
While the kinetic variables (m(t),h(t),n(t)) cannot be recorded in experiments, we can estimate them from the experimental data V(t) <cit.>. 


A total of 90 time series are generated by applying impulses of intensity Δ V =1, 2, and 3 [mV], respectively, to 30 points on the limit cycle of the HH oscillator. It should be noted that comparable time series data can be obtained in experiments where neurons are stimulated with repeated impulses.
We first estimate the phase function Θ() using the proposed method (GPPI) from the time series data. The training data for the Gaussian process regression was prepared by sampling at 0.1 [ms] intervals from the initial 3 [ms] of the time series. We then calculate the normalized phase response functions (nPRFs) in the +V direction from the estimated phase function. 

Figure <ref>a compares the true and estimated nPRFs for impulse intensities of Δ V = 1.5 and 2.5 [mV], as well as the phase sensitivity function Z_V:= ∂Θ/∂ V((θ)). 
It should be noted that Figure <ref>a shows the prediction performance of the proposed method for the impulses that were not included in the training data set. 

In addition, the phase sensitivity function Z_V was not estimated from the data, and it was obtained from the numerical simulation of the HH model with an impulse intensity of Δ V=0.1 [mV].
For a strong impulse (Δ V =2.5 [mV]), the true nPRF deviates from Z_V due to nonlinearity of the oscillator, whereas the GPPI is able to estimate the phase response accurately from data. 








Next, we synchronize the phases of two uncoupled HH oscillators (oscillator A and B) by using impulse inputs in the +V direction to the oscillator A. 
The phases of the two oscillators (θ_A(t) and θ_B(t), respectively) are obtained from the observed data. The phase information is used to determine the intensity and timing of the impulse in order to align the phases of the two oscillators. 
Due to the limitation of the impulse intensity (0 ≤Δ V ≤ 3 [mV]), multiple impulses are required if the phase difference between the oscillators is too large. In such cases, the next impulse was added to the oscillator after sufficient time has elapsed to converge to the limit cycle. 

We examine a strategy for controlling the phase of the oscillator based on the phase response function (PRF): g(θ, Δ V 𝐞_V). The strategy is as follows: 
First, if the PRF indicates that it is possible to align the phases between the oscillators with a single impulse, we adopt the minimum strength of the impulse that achieves the complete synchronization. 
Second, if the PRF indicates that the phases cannot be aligned with a single impulse, we adopt the impulse that brings the phase difference closest to zero. In other words, if the phase of oscillator A is too far ahead (behind), the impulse that delays (advances) the phase the most is selected.


This control strategy based on the PRF obtained from the estimated phase function: ĝ(θ,𝐞_VΔ V):= Θ̂((θ)+𝐞_VΔ V) - θ, was applied to the HH oscillator. 
Figures <ref>b and c compare the performance of the phase control based on the estimated phase function with that based on the phase sensitivity function, where the PRF is approximated with the phase sensitivity function: ĝ(θ,𝐞_VΔ V) ≈ Z_V(θ)Δ V. Again, the phase sensitivity function was not estimated from the data, and it was obtained from the HH model. 

Given that the maximum phase change due to an impulse is approximately 0.5, it is expected that a single impulse is sufficient when the initial phase difference is 0.3 (Fig. <ref>b and c: left). It can be also expected that two impulse inputs are required when the initial phase difference is 0.8 (Fig. <ref>b and c: right).

Figures <ref>b and c confirm that the phase of the HH oscillator can be controlled as expected from observation data alone, without the knowledge of the HH model.  

While the phase can be also controlled using the phase sensitivity function Z_V, an additional  impulse is required to achieve the complete synchronization. 
In particular, the method based on the phase sensitivity function becomes inefficient for delaying the phase (Fig. <ref>b). 
This is because the phase sensitivity function Z_V overestimates the PRF when it takes negative values (Fig. <ref>a), yielding too strong input for the control. 

This result implies that the data-driven control can be achieved by using the proposed method (GPPI), and 
the GPPI method achieves superior performance to the method based on the phase sensitivity function obtained from the HH model. 





§ DISCUSSION


In this study, we have proposed a method, Gaussian Process Phase Interpolation (GPPI), for estimating the asymptotic phase (i.e., the phase function) using time series data obtained from the limit cycle oscillator. This method is based on Gaussian process regression and is applicable even when the dimension of the dynamical system is increased. 
We have applied the method to the Stuart-Landau (SL) oscillator and the Hodgkin-Huxley (HH) oscillator. Our results demonstrate that the GPPI method can accurately estimate the phase function even in the presence of substantial observation noise and strong nonlinearity from relatively small dataset. 
Furthermore, we have demonstrated that the phase function estimated by the GPPI method is effective for data-driven phase control with impulse input.



As previously stated in the Introduction, the DPR method proposed in <cit.> is an exiting method with the same objective as the proposed method. We have shown that the GPPI method can estimate the phase function with fewer data points (Fig. <ref>) and is more robust to the observation noise than the DPR (Fig. <ref> and Table <ref>). 
The GPPI differs from the DPR mainly in two respects, which contribute to the improvement. The first difference is in the regression method. The GPPI employs Gaussian process regression, whereas the DPR is based on the polynomial regression. Gaussian processes have a higher expressive power than polynomial functions.
The second difference is in the object of the regression. The GPPI fits the functions of the asymptotic phase, i.e., sinΘ(), cosΘ(), whereas the DPR fits a differential equation that the phase satisfies (see <ref>). 
One advantage of the GPPI is that it can be applied to time series with a large sampling interval, because it does not require the derivative of the state variables.



The Phase AutoEncoder (PAE), recently proposed in <cit.>, is comparable to the GPPI method in that it estimates the phase function from the time series. Moreover, it can be applied to limit cycle oscillators from high-dimensional dynamical systems. 
One advantage of the PAE is its capacity to learn the phase function from a large amount of data. However, it requires the tuning of a large number of hyperparameters, including the number of neurons and layers, as well as the weights of individual loss functions, which is typically done through a process of trial and error. 
In contrast, the proposed method has only two hyperparameters (σ_f and σ_l), which are optimized by maximizing the likelihood function. 

It has been shown that the asymptotic phase can also be interpreted as the argument of the Koopman eigenfunction <cit.>. 
The existing methods for estimating Koopman eigenfunctions, including Dynamic Mode Decomposition (DMD) and its extensions  <cit.> are applicable for estimating the asymptotic phase.  
However, there remain challenges for estimating the asymptotic phase, 
as an accurate estimate of the Koopman eigenvalue (i.e., the natural frequency) is required. 
In contrast, we separately estimated the accurate natural frequency by using the Poincaré section in the present method.  





As an application of the phase function estimation, we have considered the problem of controlling the phase of an oscillator by impulse input. 
Specifically, we consider the HH oscillator and control the phase by using the pulse in the +V direction, assuming that the time series {V(t), h(t), m(t), n(t) } is available. 

While it is possible to stimulate the neuron using such an impulse in an experiment <cit.>, the variables (m, h, and n), referred to as kinetic variables, cannot be measured directly. 
Nevertheless, it is possible to infer the kinetic variables from the measured data <cit.>. 
Consequently, the control of limit cycle oscillators in the real world, such as neurons and the human heart and respiratory system, using the proposed method would be one of the important directions of future research. 




A limitation of this study is the computational cost required for Gaussian process regression. A potential way to improve the estimation accuracy is to increase the amount of data. 
However, the computational complexity of Gaussian process regression is O(N^3), where N is the number of data points, rendering it an impractical method when data size increases. 
Indeed, the HH oscillator required approximately 10 times as many data points as the SL oscillator for the phase function estimation, which resulted in a computation time more than 100 times longer. 

A promising avenue for future research would be to improve the Gaussian process regression algorithm for more efficient phase function estimation. For example, we could use fast approximation algorithms for Gaussian process regression <cit.>, including the FIC method <cit.>.    





For high-dimensional limit cycle oscillators, the data distribution in the state space may be non-uniform for the following reasons: 1) The dynamical system may possess multiple time constants, resulting in fast and slow parts of the limit cycle orbit, 2) The shape of the limit cycle may become complex in state space. 
It is essential to obtain the data points so that their distribution is close to uniform distribution when using the GPPI method. This is because the Gaussian process regression implicitly assumes a uniform distribution of the data. 
While the GPPI can accurately estimate the phase function of the HH oscillator (Fig. <ref> and <ref>), it requires proper data sampling. 
When we collect all the time series, the data tends to be concentrated in the vicinity of the limit cycle. Additionally, the data density is relatively low in regions where the state variable changes rapidly, such as the spike in the HH oscillator. 
It would be an interesting topic of future research to extend the GPPI so that it can be easily applied to a variety of limit cycle oscillators. 



In this study, we have presented the Gaussian Process Phase Interpolation, a method for estimating the phase function of the limit cycle oscillator from time series data. We have demonstrated the utility of the proposed GPPI method by applying it to simulated limit cycle oscillators. Our future goal is to apply the GPPI to limit cycle oscillators in the real world and develop a methodology for controlling their synchronization dynamics. To this end, we will refine the phase estimation method so that it can handle more complex oscillatory systems.






§ ACKNOWLEDGEMENTS


We thank Norihisa Namura for sharing us the code used in  <cit.>. We thank Hiroshi Kori and Massimiliano Tamborrino for fruitful discussion. 
This study was supported by
the World-leading Innovative Graduate Study Program in Proactive Environmental Studies (WINGS-PES), the University of Tokyo, to T.Y., 
JSPS KAKENHI (Nos. JP22K11919 and JP22H00516) and Japan Science and Technology Agency CREST (No. JPMJCR1913) to H.N., and 
JSPS KAKENHI (Nos. JP18K11560, JP21H03559, JP21H04571, JP22H03695, and JP23K24950) and AMED (Nos. JP21wm0525004 and  JP223fa627001) to R.K. 
















§ DERIVATIVE PHASE REGRESSION METHOD




We briefly describe the derivative phase regression method (DPR method) proposed by <cit.>. 
The DPR method estimates two functions, s():=sinΘ() and c():=cosΘ(), as the proposed method (the GPPI method). 
In contrast to the GPPI, which regressed the functions (s() and c()) directly using Gaussian processes, 
the DPR regresses the time derivatives of these functions using the polynomial function. 

Here, we describe an overview of the DPR method (see <cit.> for details). 
Taking the time derivatives of the the functions (s() and c()) using Eq. (<ref>), we obtain

    ∇ s((t)) ·(t)    = ω c((t)), 
    ∇ c((t)) ·(t)    = -ω s((t)),

where ∇ denotes the gradient of a function, and · denotes the inner product. 
Let us approximate these functions with a polynomial function, 

    s()    = ∑_ia_kf_k(),
    
        c()    = ∑_ib_kf_k(),

where f_i() is a polynomial function of the state variables (x_1, x_2, ⋯ ,x_d). 
The parameters { a_k, b_k} are determined to satisfy Eqs. (<ref>) and (<ref>), which is achieved by minimizing the error function E( {a_k}, { b_k } ), 

    E:= ∑_i(∑_ka_k∇ f_k(_i)·_i - ω∑_kb_kf_k(_i))^2
        
     + ∑_i(∑_kb_k∇ f_k(_i)·_i + ω∑_ka_kf_k(_i))^2
        .

We also require a constraint to fix the origin of the estimated phase function ( Θ(_0)= 0 ), i.e., 

    s(_0)    = ∑_ia_kf_k(_0) = 0,
    
        c(_0)    = ∑_ib_kf_k(_0) = 1.

The minimization problem of the error function E( {a_k}, { b_k } ) with the constraint (Eqs. <ref> and <ref>) is a quadratic programming. Thus, we can determine the parameters { a_k, b_k} using an efficient solver, such as quadprog in MATLAB. 

Finally, the phase function Θ() is obtained from the estimate of the functions (s(), c()) using Eq. (<ref>).  





§ ESTIMATION ERROR OF THE PHASE FUNCTION
 

We evaluate the estimation performance of the asymptotic phase function by calculating the error: E()= |Θ̂() - Θ()|, where Θ̂() and Θ() represent the estimated and true phase functions, respectively.  
Figure <ref>a shows the data used for estimation and the estimation error in the phase function for the Stuart-Landau oscillator, which corresponds to Fig. <ref> (Sec.<ref>). While the DPR method used all the 50,000 data points for estimating the phase function, the proposed method (GPPI) used only 1000 data points sampled from the data. The strength of the observation noise was η=0.005. 

Figure <ref>b shows the data used for estimation and the estimation error in the phase function for the Hodgkin-Huxley oscillator, which corresponds to Fig. <ref> (Sec.<ref>).







elsarticle-harv















