Gaussian process regression with additive periodic kernels for two-body interaction analysis in coupled phase oscillators

Let us consider a system comprising $N$ interacting dynamical systems:

Here, the function $F_{i}$ represents the spontaneous dynamics of the element $i$, and $G_{j\to i}$ represents a weak action from the element $i$ to the element $j$. We assume that each element has a stable limit cycle (i.e., each element is an oscillator). It is well known that such a system can be reduced to coupled phase oscillators [1, 2].

The functions $\{\Gamma_{ij}\}$ are coupling functions, which represent the action from the oscillator $j$ to the oscillator $i$.

The advantage of employing phase-oscillator model lies in the significant simplification of the system. While the original system (1) is generally multidimensional, the phase oscillators (2) assign only one variable $\phi_{i}$ to each limit cycle. Consequently, the focus can be directed solely towards rhythm synchronization characteristics.

In this paper, we estimate all phase coupling functions $\{\Gamma_{ij}\}$ from phase time series $\mathcal{D}=\{(\bm{x}_{j},y_{j})\}_{j=1}^{n_{\mathrm{data}}}$ with $\bm{x}_{j}\in\mathbb{T}^{N-1}$ and $y_{j}\in\mathbb{R}$ for $j=1,2,\ldots,n_{\mathrm{data}}$. Here, $\mathbb{T}$ represents a torus. Since the phase coupling functions determine synchronization characteristics of the oscillators, it is desirable to estimate these functions from observed time series.

However, obtaining time-series data that comprehensively spans a high-dimensional state space is exceedingly challenging. Even if such data were acquired, the amount of data would be enormous, necessitating significant computational resources for estimation. To avoid this difficulty, we focus on the fact that the phase coupling functions are represented as the sum of two-body interactions, which enables us to use an additive kernel. See section 3.1 for details.

The covariance function $k$ is a fundamental component of Gaussian process regression, as it encodes our assumptions about the function which we wish to predict (see [11] for a detailed discussion). In our problem setting, since the coupling function have a domain of torus $\mathbb{T}^{N-1}$, which means the function is $2\pi$-periodic, it is appropriate to use a periodic kernel as the covariance function. In addition, since the coupling function is additive (i.e., it is represented as the sum of two-body interactions), it can be seen that the kernel is also additive [12]. In conclusion, we propose to use following additive periodic kernel to estimate the coupling functions:

Here, $\{\theta_{0}^{(j)}\}_{j=1}^{N-1}$ and $\{\theta_{1}^{(j)}\}_{j=1}^{N-1}$ are positive hyperparameters. $\theta_{0}^{(j)}$ determines the amplitude of the covariance function. $\theta_{1}^{(j)}$ can be said as an inverse lengthscale, which specifies the width of the $j$-axis direction of the covariance function.

In the Gaussian process regression described so far, hyperparameters remain in the kernel $\theta_{0,1}$ and the noise strength $\sigma$. For brevity of notation, we will write these parameters as $\bm{\theta}$.

To estimate $\bm{\theta}$, we consider the maximum likelihood estimation, which makes inferences about the population that is most likely to have generated the data $\mathcal{D}$. The (marginal) log-likelihood for parameters $\bm{\theta}$ is

where $K_{\bm{\theta}}=k_{XX}+\sigma^{2}I_{n_{\mathrm{data}}}$ is the covariance matrix for the noisy data $\mathcal{D}$. The most probable parameters $\bm{\theta}$ are, therefore, estimated by finding the maximum of the log-likelihood function $\mathcal{L}(\bm{\theta})$.

The problem of finding the point that maximizes the marginal likelihood while changing hyperparameters is formulated as an optimization problem. In this paper, we use the gradient descent method is to find the maximum point of the marginal likelihood. The hyperparameters $\bm{\theta}$ is updated using the gradient with the following manner:

where $\alpha$ denotes the learning rate.

Stochastic gradient decent (SGD) method, regarded as a stochastic approximation of gradient descent optimization methods, is another way of optimization. In the SGD method, the gradient is determined by different minibatches at each step, and parameter updating proceeds accordingly. The greatest benefit is that it requires less computation time to converge than the gradient method, which uses all data at each step. It is also thought to be less likely to be trapped in local solutions, thanks to stochastic fluctuations in the gradient direction for each minibatch. Note that we can calculate the derivative of the covariance function analytically with respect to their hyperparameters. This enables efficient optimisation using automatic differentiation. The implementation was carried out using GPflow [13].

One of the main merits of Gaussian process regression is its ability to provide uncertainty estimates for the regression parameters, which can be useful in situations where the data are noisy or sparse. It also allows for the modeling of complex functions without requiring the specification of a fixed functional form, which can be difficult to determine in many cases. Additionally, the use of a kernel function allows for the incorporation of prior knowledge about the function, such as its smoothness or periodicity, which can improve the accuracy of the estimates.

Gaussian process regression has the advantage that all calculations can be done in a closed form using only matrix operations. However, as the number of data $n_{\mathrm{data}}$ increases, inverse matrix calculations for $n_{\mathrm{data}}\times n_{\mathrm{data}}$ matrices are required, and the amount of memory and computation is enormous. Many sparse approximations have been proposed to overcome the computational complexity of Gaussian process regression [14, 15, 16, 17, 18].

To confirm the validity of our approach, we first consider a simple 3-body system of coupled phase oscillators:

where $\phi_{1},\phi_{2},\phi_{3}$ are the phases of the oscillators. We prepare the time series by solving the ODEs (11) with forward Euler method from $t_{0}=0.0$ to $t_{1}=15.0$ using a randomly selected initial value and sample the data with sampling rate 0.2. We repeat the above calculation 50 times. Estimation results are shown in Figure 1.

Next, we apply the proposed method to coupled Van der Pol oscillators:

with $\langle\xi_{\alpha}(s)\xi_{\beta}(t)\rangle=\sigma^{2}\delta_{\alpha,\beta}\delta(s-t)$ for $\alpha,\beta\in\{x_{1},y_{1},x_{2},y_{2}\}$. The values of parameters are $\varepsilon_{1}=0.3,\varepsilon_{2}=0.7,K=0.01,\sigma=0.03$.

We reconstruct the phase time series from the obtained time-series data using the following method. First, we calculate the geometric angle $\vartheta_{i}(t)$ as follows:

where $\hat{x}_{i}$ and $\hat{y}_{i}(t)$ are centered time series of $x_{i}(t)$ and $y_{i}(t)$, respectively. Even though $\vartheta_{i}(t)\in[0,2\pi)$ for all $t$, it is inappropriate to consider $\vartheta_{i}(t)$ as the phase time series of the oscillator $i$. This is because $\vartheta_{i}$ does not vary monotonically in the absence of interactions and noise. To overcome this difficulty, the following transformation from $\vartheta_{i}$ to $\theta_{i}$ should suffice [19, 20]:

where $\sigma_{i}$ denotes the empirical distribution of $\vartheta_{i}$. We consider $\theta_{i}(t)$ as the phase of oscillator $i$ and use them in the regression.

Estimation results are shown in Figure 2. The results are highly consistent with the true coupling functions derived using the adjoint method. For comparison purposes, we also present the estimated coupling functions with the existing method [9].

Here is another example which demonstrates the effectiveness of the proposed method. Consider the following coupled FitzHugh-Nagumo oscillators [21, 22]:

with $\langle\xi_{\alpha}(s)\xi_{\beta}(t)\rangle=\sigma^{2}\delta_{\alpha,\beta}\delta(s-t)$ for $\alpha,\beta\in\{v_{1},w_{1},v_{2},w_{2}\}$. Here, $a=1,b=0.8,I_{\mathrm{ext}}=0.8,K=0.0087,\tau_{1}=0.1^{-1},\tau_{2}=0.09^{-1}$.

We estimated the coupling functions for various numbers of data (Figure 3). Specifically, we varied the number of initial values used to create the data (denoted by $N_{\text{init}}$) and the number of oscillator rounds per initial value used for estimation (denoted by $N_{\text{cycle}}$). We set the sampling rate as 0.2. The phase time series are reconstructed using the same method as described in section 4.2.

From Figure 3, it is evident that when the data sample size is small, our proposed method accurately aligns with the theoretical result in estimating the coupling function, whereas the existing method fails in estimation. This implies that our method is superior to the other methods when a number of the data are small and the data are biased.

Finally, we consider a network of multiple coupled Hodgkin-Huxley models [23] and estimate the coupling function for all of them. There are seven neurons, five of which are excitatory and two of which are inhibitory. Details of the Hodgkin-Huxley model is summarized in A.

We employed a kernel function in the form of equation (8) with $N=7$. See the lower graphs of Figure 4 and Figure 5 for the estimation result compared to the theoretical coupling functions obtained by the phase reduction approach. We confirm that the coupling functions between each oscillators are qualitatively consistent with the results from Gaussian process regression and the theoretical results.