Ethics statement

The scalp EEG experiment was approved by the ethics committee of the Unit of the Integrated Studies of the Human Mind, Kyoto University (24-p-19). Participants provided written informed consent according to the Declaration of Helsinki and were paid for their participation.

Estimation of phase coupling functions for cross-frequency synchronization

Neural oscillations can be observed easily from EEG data, and many EEG studies have reported various types of synchronization [8, 38], which can be roughly divided into same- and cross-frequency synchronization. A few experimental results suggest that same-frequency phase-phase synchronization plays a role in modulating neuronal interaction [6, 7]. In contrast, cross-frequency synchronization is considered to play a role in the integration of activities over different time scales [18]. However, it is unclear how these synchronizations, particularly 1:p phase synchronization, are achieved by a dynamical system. Therefore, we developed an effective method to identify the dynamical system that performs these synchronizations.

In general, synchronization of neural oscillation is thought to be described by a network of limit-cycle oscillators, which can be described generally by the multidimensional differential equation , where X_i denotes the multidimensional state of the i-th oscillator, such as membrane voltages and gate variables of ionic channel. We assume that a system X_i can generate a limit-cycle oscillation by itself without external interaction. An EEG signal is thought to be generated by some neuronal system consisting of many interacting neurons. In this context, it is plausible that the neuronal system of an EEG signal can be represented by the X_i system. According to the phase reduction method, the limit-cycle oscillator can be characterized theoretically by a phase ϕ as a simple dynamical system with one degree of freedom. If the oscillators are weakly coupled, the dynamics of the networks among N oscillatory systems can be described by [27, 39]: (1) where ω_i is the natural frequency of the oscillator and Γ_i,j is a phase coupling function representing the influence from the j-th oscillator to the i-th oscillator. It is known theoretically that the phase coupling function depends only on the phase difference ϕ_j−ϕ_i. When the phase difference is constant over time, these oscillators are said to be synchronized. Specifically, the synchronization of same-frequency oscillators is referred to as 1:1 phase locking. Eq (1) can describe the 1:1 phase-locking state between rhythms in real systems.

Various synchronizations between slow and fast oscillators, e.g., theta (4–8 Hz) and gamma (>30 Hz) EEG activities, have been observed ubiquitously, and they appear to play an important role in brain function [8–11]. In fact, 1:p phase locking has been observed ubiquitously in human EEG experimental studies during the resting state, mental arithmetic tasks, and working memory tasks [12–15]. However, 1:p phase synchronizations cannot be described by the model expressed by Eq (1). Therefore, we consider the 1:p phase-locking state among slow and fast oscillators. If 1:p phase locking occurs, the value of the phase difference ϕ₂−pϕ₁ is constant over time, where ϕ₁ and ϕ₂ are slow and fast phases, respectively. Using the phase reduction theory, we found that 1:p phase locking can be described as follows [39–41]: (2) (3) (4) Here, we explain a simple case of two coupled oscillators described by Eqs (2)–(4). Note that many real rhythmic systems generally consist of many oscillators. We assume that the ratio of the natural frequencies of the two oscillators is close to some integer p. Note that, in this situation, the coupling function Γ_1,2 depends on only the phase difference ϕ₂−pϕ₁.

In our approach, to investigate the nature of interactions between neuronal rhythms, we directly estimate both the natural frequencies ω_i and the phase coupling functions Γ_i,j from experimental time-series data. In addition, considering unavoidable sources of uncertainty (e.g., observational error or additional unknown disturbance to the system), we introduce independent Gaussian white noise η_i(t) into the phase oscillator model as follows: (5) (6) Here, we assume that the independent Gaussian white noise η_i(t) satisfies ⟨η_i(t)⟩ = 0,⟨η_i(t)η_j(t′)⟩ = 2D_iδ_ijδ(t − t′), where δ_ij and δ(t) are the Kronecker delta and the Dirac delta functions, respectively. D_i indicates the noise strength and p_iϕ_j−p_jϕ_i denotes the phase difference, where the p values are integers. Note that this phase oscillator model can explain p_i:p_j synchronization (e.g., 2:3 phase synchronization and 2:7 phase synchronization). We estimate the phase oscillator model (Eq (6)) using almost periodic time-series data. In the following, we employ a straightforward extended version of a previously proposed method [28] and explain the outline of our method.

First, we transformed the experimentally-recorded signal s(t) into the phase time-series θ(t) by computing the analytic signal as follows: (7) where denotes the Hilbert transformation of the recorded signal s_i(t) [42], and θ(t) is the phase of the analytic signal. However, the variable θ is generally different from the phase ϕ in Eq (1) because, according to phase reduction theory, ϕ evolves linearly over time without interaction and noise. It is therefore necessary to transform θ into ϕ, as follows [30, 31]: (8) where f(θ) denotes the probability density distribution of θ.

Second, Bayesian linear regression [43, 44] is applied to estimate the parameters of the phase oscillator model given by Eq (6). Because the coupling function is periodic, we consider Fourier series expansion of the coupling function as: (9) where ψ_i,j(t_τ) is the extended version of the phase difference p_iϕ_j(t_τ) − p_jϕ_i(t_τ) at time t_τ. The times t_τ are discrete points, t_τ = t₁ + (τ − 1)Δt for τ = 1,2,⋯,T, and Δt is the sampling interval. In this expansion, M_ij denotes the Fourier series order for each coupling function, and the parameters M_ij control the complexity of the coupling function. The parameters M_ij can be determined by a model selection method based on logarithmic evidence, as explained later.

Finally, the proposed method estimates the model as follows: (10) where . Thus, the unknown model parameters are and D_i. Here, denotes . The phase velocity is a dependent variable in a standard linear regression problem that is computed from phase time-series data as {(ϕ_i(t_τ+1) − ϕ_i(t_τ))/Δt}. Furthermore, an independent variable is computed by the phase difference as {cos(mψ_ij(t_τ)),sin(mψ_ij(t_τ))}_i. Here, η_i is an independent and identically distributed random variable. This linear regression problem corresponds to maximizing the following likelihood function: (11) where equals , and N(x|μ,σ²) denotes Gaussian distribution, where μ and σ² are the mean and variance of x, respectively. Using Bayesian theory, the product of the likelihood function and the prior distribution is proportional to the posterior distribution : (12) Here, we adopt a Gaussian inverse gamma distribution for the prior distribution . This prior distribution is a conjugate to the likelihood function (the so-called conjugate prior): (13) where IG(x|α,β) denotes the inverse gamma distribution with shape parameter α and scale parameter β. and are the mean and covariance of model parameters , and b^(m), respectively. Note that the prior parameters , and are referred to as hyperparameters. We can easily calculate the posterior distribution parameters from the conjugate prior distribution and the likelihood function in Eq (12) (S3 File). The posterior distribution with the updated parameters takes the form of a Gaussian inverse gamma distribution: , and . The estimated model parameters are the mean of . The estimated noise level is the mean of . Using the posterior, prior, and likelihood functions, we can compute the logarithmic evidence log p({ϕ_i(t_τ)}): (14) Thus, we find M_ij with the largest logarithmic evidence among all models.

Electronic circuit of van der Pol oscillator

To test the proposed method, we conducted an experiment in which an electronic circuit was used to implement two coupled van der Pol oscillators [45]. The coupling function Γ_i,j between the oscillators can be obtained theoretically from the corresponding differential equations.

In this experiment, we recorded the rhythmic signals of the electronic circuit. Each oscillator consisted of two multipliers U₁ (AD633, Low Cost Analog Multiplier) and three operational amplifiers U₂ (TL082, ½ Dual BiFET Op Amp) (Fig 1A). We conducted two experiments under different conditions. In the first experiment, two same-frequency oscillators were coupled directly. In this experiment, 1:1 phase locking was expected to occur. The parameters of the electronic component were set to R₁ = 100kΩ, R₂ = 1kΩ, R_coupling = 1MΩ, C₁ = C₂ = 0.01μF, V₁ = 0.115V, and V₂ = 0.12V. R_k and C_i are the parameters of the resistor and capacitor, respectively. Voltages V_i were monitored using a digital voltmeter. In the second experiment, a slow oscillator was coupled to a fast oscillator. In this experiment, 1:2 phase locking was expected to occur. The natural frequencies of the slow and fast oscillators were set to satisfy a nearly 1:2 ratio. The parameters of the electronic components were the same as those in the first experiment, except that C₁ was changed to C₁ = 0.02μF to reduce the natural frequency by one-half. In both cases, the electric potentials x_i and y_i were recorded using an I/O device (NI SCB-68, National Instruments, US). The sampling rate of the electric potential was 15,000 Hz, and the data size was 180 s.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Electronic circuit of a pair of van der Pol oscillators and recorded electric potential.

(a) Schematic of electronic circuit of two coupled van der Pol oscillators, where x_i and y_i are positions for recording electric potential, R_k denotes resistors, and C_i denotes condensers. Electronic units U₁ and U₂ represent the multiplier and operational amplifiers, respectively. R_coupling is a resistor whose resistance is the parameter of the strength of connectivity. (b) Experimental data of electric potentials x₁ and y₁ show the limit-cycle oscillator under the same-frequency (129.1 Hz) coupling condition (gray dots and line). The black trajectory shows the theoretical value computed by the van der Pol oscillator Eqs (15–18). Here, the frequency is 142.1 Hz. Blue dots represent the zero-phase reference points on the experimental data, which were determined automatically via Hilbert transformation. Green crosses represent the theoretical zero-phase reference points defined as the peak points of x_i. Red dots denote the adjusted zero-phase reference points. (c) x₂ and y₂ show the oscillators under the same-frequency oscillator condition. The frequency of the experimental data is 132.5 Hz and that of the theoretical trajectory is 146.4 Hz. (d) Recorded electric potentials show the slow limit-cycle oscillator under cross-frequency coupling conditions (experimental frequency, 64.1 Hz; theoretical frequency, 71.1 Hz). (e) x₂ and y₂ denote the fast oscillator (experimental frequency, 131.1 Hz; theoretical frequency, 146.4 Hz).

https://doi.org/10.1371/journal.pcbi.1005928.g001

The corresponding theoretical equation of the electronic circuit is given as: (15) (16) (17) (18) where x_i and y_i are the corresponding theoretical electric potentials of the i-th oscillator. The trajectories of the van der Pol oscillator are shown by the x_i and y_i signals in Fig 1B–1E. Note that only x_i was used to estimate the coupling function. The experimental parameters (R_k and C_i) and those in Eqs (15–18) were the same as the parameters of the electronic components. Note that the term b does not exist in the original van der Pol oscillator. In the case of b = 0, the theoretical trajectory and coupling function do not agree with the experimental data and the estimated coupling functions, respectively. Note that the original van der Pol oscillator generates a symmetrical trajectory. However, in the electronic circuit experiment, the trajectory was not exactly symmetrical due to small additional disturbances in the system or the uneven quality of the electronic components. By introducing parameter b (b = 0.4), the simulated trajectory can be adjusted to the experimental data (Fig 1B–1E).

The corresponding theoretical coupling functions were calculated based on Eqs (15–18) using the adjoint method [46]. In this method, the zero phase is defined as the peak point of x_i (green crosses in Fig 1B–1E represent the zero-phase reference points on the theoretical orbits). From the experimental data, the zero-phase reference points were determined based on electric potentials x_i using Hilbert transformation (blue dots denote the point of zero-phase reference on the experimental data). As shown in Fig 1B–1E, a small gap exists between the zero-phase reference points of the theoretical and experimental orbits. In principle, an arbitrary point on the limit-cycle orbit can be defined as the zero-phase point. However, to compare theoretical and estimated results, the zero-phase reference point of the theoretical model should be consistent with that of the experimental data. Therefore, the zero-phase reference points of the experimental orbits were adjusted to coincide with those of the theoretical model by shifting the experimental reference points to the theoretical points as (red dots denote the revised zero-phase reference points), where ϕ₁ is the experimental phase based on the electric potential x₁ (Fig 1B and 1D) and ϕ₂ is the experimental phase based on the electric potential x₂ (Fig 1C and 1E). We estimated the coupling function from the shifted phase data and compared the theoretical coupling function with the estimated function.

Scalp EEG experiment

We applied the proposed method to scalp EEG data. The method can estimate the coupling functions for same-frequency and cross-frequency synchronization assuming that EEG activities can be considered weakly coupled oscillators. However, it is unclear whether EEG activity can be considered a weakly coupled oscillator system. Thus, we must confirm that the proposed method can estimate the dynamical system from the EEG data successfully.

We used EEG data recorded during a speech recognition task. Note that detailed information is provided in our previous paper [37]. The participants categorized what they heard as a target or distractor as soon and as accurately as possible. Four-letter Japanese words were used as speech sounds, and the words were uttered within approximately 1 s. The sampling rate of the speech sound was 48,000 Hz. The speech envelopes on each frequency were high-pass filtered with a cutoff frequency of 3 Hz to avoid phase-resetting. Furthermore, the speech sounds were masked with pink noise. The noise volume was increased linearly over 0.5 s after onset to avoid the phase-resetting effect by noise onset. The speech sound always started 2 s after the onset of noise and lasted approximately 1 s. The noise sound was terminated 1.5 s after speech onset. The EEG experiment consisted of four sessions for each participant. Each session consisted of 100 trials.

A 32-channel EEG amplifier (Brain Amp MR, Brain Products, Germany) with an international 10% standard electrode cap with a sintered Ag/AgCl ring electrode (Easy Cap, Falk Minow Services, Germany) was used for the EEG recording (sample rate, 5 kHz). Four electrodes were used for the vertical and horizontal electrooculogram (VEOG and HEOG) channels. The VEOG and HEOG were used to remove ocular artifacts. The measurement reference was linked earlobes, and the ground was on the inion. The EEG signal was filtered using a 1-Hz high-pass software filter, a 250-Hz low-pass hardware filter, and a 60-Hz notch filter. In a preprocessing step, ocular artifacts were corrected using EEG analysis software (Brain Vision Analyzer, Brain Products, Germany) and the VEOG/HEOG signals [47]. The reference was changed to the average of all electrodes, except VEOG and HEOG. The preprocessed EEG data were then downsampled to 500 Hz.

The participants were 16 healthy adults (five females; 11 males; 21–32 years; mean age, 25 years). One participant was excluded due to a low response rate, and another participant was excluded due to an excessive artifact that could not be removed during preprocessing. Note that these participants were also excluded in our previous study [37].

We estimated the phase oscillator model between the 3–6 Hz EEG (theta oscillation) and the speech envelope. The theta oscillation is synchronized with the envelope of speech sound, and it plays an important role in speech processing [35, 48, 49]. Speech rhythm consists of linguistic components, e.g., syllabic and prosodic rhythms. A syllable is a unit of speech that separates a word into sound chunks. For example, the Japanese word “KaKuShiKi” (“formality” in English) is composed of four syllables “/Ka/Ku/Shi/Ki/,” and its sound envelope appears in the 4–5 Hz frequency range in the current speech stimulus. Prosody is the stress and intonation patterns of an utterance. The envelope of prosodic rhythms appears in the 1–3 Hz frequency range. We estimated the coupling function between the EEG oscillation and these speech rhythms. To obtain the theta oscillations, the preprocessed EEG signals were bandpass filtered within 3–6 Hz. The syllabic and prosodic rhythms were computed from the sound stimulus consisting of the noise and speech sounds (Fig 2A). The sound envelope was computed as the absolute value of the Hilbert-transformed sound data (Fig 2B). To compute the syllabic and prosodic rhythms, the envelope signal was bandpass filtered within 3–6 Hz and 1–3 Hz, respectively (Fig 2C and 2D). The syllabic and prosodic signals were downsampled to 500 Hz. The instantaneous phases of these rhythms were computed using both Hilbert transformation and correction (Eq 8). We estimated the coupling functions of 1:1 phase locking (syllable and theta oscillation) and 1:2 phase locking (prosody and theta oscillation).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Syllable and prosody rhythms in speech sound.

(a) Example of speech stimulus. The stimulus consisted of noise and a four-syllable Japanese word. The red line represents a speech wave. The blue line represents the presented sound wave, which consists of speech plus noise sounds. (b) Speech envelope was computed as the absolute value of Hilbert-transformed speech sound. (c) Syllabic rhythms were computed from the speech envelope through the bandpass filter within 3–6 Hz. (d) Prosodic rhythms were computed from the speech envelope through the bandpass filter within 1–3 Hz.

https://doi.org/10.1371/journal.pcbi.1005928.g002

Numerical simulation of phase oscillator model

First, we applied our Bayesian method to numerical simulation data which was generated from three cross-frequency oscillators with somewhat complicated connections, to evaluate the validity of the proposed method. Simulation data were generated from a network comprising one fast oscillator and two slow oscillators (Fig 3) based on the Euler–Maruyama method [50] using the following differential equations: (19) (20) (21)

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Estimated coupling function for numerical simulation data.

Upper-left diagram shows the network structure. The estimated coupling functions (red lines) were nearly identical to the correct functions (dashed black line). The gray dots represent the phase time-series data. When the interaction did not exist, the estimated coupling function was identically zero. The proposed method estimated all coupling functions correctly for the simulation data.

https://doi.org/10.1371/journal.pcbi.1005928.g003

We set the natural frequencies to ω₁ = 0.9, ω₂ = 2.1, and ω₃ = 1.1. Here, ϕ₂ is the phase of the fast oscillator, and ϕ₁ and ϕ₃ are the phases of the slow oscillators. η is independent Gaussian white noise with zero mean and a standard deviation of 0.1. Using the proposed method, we estimated the coupling functions and the natural frequencies from the phase time-series data generated in the experiment.

Fig 3 shows the estimated coupling functions and the correct coupling functions. In this case, the correct coupling functions were defined explicitly by Eqs (19)–(21). Despite the complicated connections, the results indicate that the estimated and correct coupling functions agree fairly well. Furthermore, the complexity parameter of the coupling function was selected correctly by optimizing the logarithmic evidence. Therefore, the proposed method works quite well at estimating a nontrivial network of phase oscillators comprising oscillators with different natural frequencies.

Electronic circuit experiment

Before applying the proposed method to EEG data, we recorded the electric potential of the van der Pol electronic circuit and tested the ability of the estimation method using the experimental data. Since the electronic circuit can be explained by the corresponding theoretical differential equations (Eqs 15–18), we can derive the correct coupling function theoretically using the adjoint method. We conducted two experiments. One involved coupling oscillators of the same frequency (Fig 4A), and the other involved coupling slow and fast oscillators (Fig 4E). We transformed the x₁ and x₂ electric potentials of the first and second oscillators, respectively, to phase time-series data and estimated the coupling functions from these data.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Estimated coupling function of electronic circuit.

(a) The diagram shows the coupling direction between oscillators of the same frequency. The first oscillator was coupled to the second oscillator. (b) The red line shows the estimated phase coupling function with the natural frequency in the same-frequency coupling case. The dashed black line shows the theoretical coupling function. The coupling function from the second to first oscillator Γ₁₂ is identically zero. When there is no interaction, the coupling function is nearly zero. The gray dots show the experimental data points. (c) The coupling functions from the first to second oscillator Γ₂₁. (d) The blue line shows the phase difference histogram of the experimental data in the case of 1:1 phase locking (experimental histogram). The red line shows the simulated histogram calculated in the phase oscillator model estimated from the experimental data (estimated histogram). The dashed black line shows the simulated histogram calculated in the phase oscillator model using the theoretical natural frequencies and coupling functions (theoretical histogram). (e) In the cross-frequency coupling case, the slow oscillator was coupled to the fast oscillator. (f) The coupling function from the fast to slow oscillator is identically zero. (g) The coupling function from the slow to fast oscillator. (h) The experimental, estimated, and theoretical histogram in the 1:2 phase-locking case.

https://doi.org/10.1371/journal.pcbi.1005928.g004

The estimated coupling functions with no interaction from the second to first oscillator were identically zero in the experiments involving oscillators of the same frequency (Fig 4B) and cross frequency (Fig 4F). When coupling existed, the estimated coupling function was the same as the theoretical function under the same-frequency (Fig 4C) and cross-frequency conditions (Fig 4G). Furthermore, to confirm whether the estimated phase oscillator model can explain the experimental data, we compared the phase difference histograms for the experimental phase data and the two types of simulated phase data. The experimental phase data were computed from electric potentials by both Hilbert transformation and Eq (8). One set of simulated phase data was computed in the phase oscillator model based on the estimated parameters using the Euler–Maruyama method. The other was computed based on the theoretical coupling functions and natural frequency, which were determined by Eqs (15)–(18). We computed the experimental, estimated, and theoretical histograms from these phase time-series data. The experimental and estimated histograms of the phase difference were nearly the same under each condition (Fig 4D and 4H). However, the theoretical histograms differed. The difference among these histograms was caused by the difference of the natural frequency between the theoretical and experimental oscillators (Fig 1B–1E) due to uncontrollable experimental conditions. In other words, the electronic circuit experimental data did not follow the theoretical equations exactly; however, the coupling functions between the experimental van der Pol oscillators were the same as the theoretical coupling functions. In fact, when the theoretical histograms were computed based on the theoretical coupling functions and the estimated natural frequencies rather than the theoretical natural frequencies, the experimental and estimated histograms coincided relatively well with the theoretical histograms. These results indicate that the proposed method works well with real data, even if the data contain observational errors or additional unknown disturbances.

Human EEG experiment

Finally, we applied the proposed method to the EEG data. We estimated the coupling functions between the theta oscillation and the envelope of speech stimulus. The theta oscillation was observed to be synchronized with the speech envelope [35] and is considered to play an important role in speech perception. Generally, the speech envelope consists of syllabic (3–6 Hz) and prosodic (1–3 Hz) rhythms. Thus, we estimated the phase oscillator model under the same-frequency and cross-frequency conditions, i.e., theta oscillation and syllabic rhythm (Fig 5), and theta oscillation and prosodic rhythm (Fig 6). In these experiments, the phase-difference histograms between the EEG and the speech sound showed 1:1 and 1:2 phase locking (Figs 5A and 6A). Note that there is obviously no interaction from EEG activity to speech sound. Therefore, we can use this fact to confirm the validity of the estimation results.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Estimated distribution of phase difference between EEG data and syllable envelope.

(a) Experimental histogram of phase difference between the theta oscillation on the Cz electrode and the syllabic rhythm (the histograms show phase locking). The gray lines represent histograms of individual participants and the blue line represents the histogram averaged over all participants. (b) Histograms obtained from the simulated data in the estimated phase oscillator model. The averaged histogram is similar to the experimental histogram. The gray lines represent the phase difference histograms of individual participants. The red line represents the average of the simulated histograms. (c) Blue lines represent the averaged experimental histogram and the standard error of mean (SEM). Red lines represent the averaged simulated histograms and the SEM. (d) Estimated coupling functions Γ_θ,s from syllabic rhythm to theta oscillation. The gray and red lines represent the results of individual participants and the average results of all participants, respectively. (e) Estimated coupling functions Γ_s,θ are considerably smaller than the opposite directional coupling functions. (f) Simulated histograms where the coupling functions Γ_s,θ are removed. The effect on the original phase-locking state was negligible. (g) Histograms where Γ_θ,s were removed are nearly flat.

https://doi.org/10.1371/journal.pcbi.1005928.g005

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Estimated distribution of phase difference between EEG data and prosody envelope.

(a) Experimental phase difference histograms for 1:2 phase locking. (b) Simulated histograms based on the estimated phase oscillator model. (c) Blue lines represent the average and SEM of experimental phase difference histograms. Red lines represent the average and SEM of simulated histograms. (d) Estimated coupling functions Γ_θ,p. (e) Estimated coupling functions Γ_p,θ. (f) Simulated histograms where coupling functions Γ_p,θ are removed. (g) Simulated histograms where Γ_θ,p is removed are uniform.

https://doi.org/10.1371/journal.pcbi.1005928.g006

In the theta oscillation and syllable data, we assumed that the syllabic rhythm modulated the theta oscillation. The instantaneous phase of the theta oscillation is denoted ϕ_θ, and the phase of the syllable is denoted ϕ_s. The phase difference between the theta oscillation and the syllable is defined as ϕ_s − ϕ_θ. The phase-difference histogram of each participant showed 1:1 phase locking (Fig 5A). We estimated the phase oscillator models of 1:1 phase locking as follows: (22) (23) To confirm that the estimated dynamical system can explain the experimental data, we simulated phase synchronization based on the estimated phase oscillator model. Fig 5B shows the phase difference histograms obtained from the simulated data. Note that, for the simulated data, phase differences were calculated by numerical simulation performed using the estimated phase oscillator model (Eqs (22 and 23)). The histograms of the simulation data were similar to those of the experimental data (Fig 5C), which indicates that the estimated phase oscillator model can explain the experimental data. Furthermore, the estimated coupling functions Γ_θ,s were consistent among all participants (Fig 5D). In contrast, the estimated coupling functions Γ_s,θ had smaller amplitudes than Γ_θ,s (Fig 5E) and were not consistent among all participants. These results are reasonable in terms of the relationship between EEG and speech sound because direct interaction from theta oscillation to speech sound never exists. To examine the effects of each coupling function on the dynamics, we computed the simulated histogram under the condition that either Γ_s,θ or Γ_θ,s was set to identically zero. In the case of Γ_s,θ = 0 (Fig 5F), the resultant histogram shows that, compared to the original dynamics in Fig 5B, the synchronized state is almost maintained. This implies that the coupling function Γ_s,θ does not contribute to the realization of 1:1 phase locking. In the case of Γ_θ,s = 0, the synchronized state disappeared, as shown in the flatter histograms (Fig 5G). Consequently, the results indicate that the coupling function Γ_θ,s primarily contributed to 1:1 phase locking.

In the theta oscillation and prosody data, we assumed that the prosodic rhythms modulated the theta oscillation. Here, let ϕ_p denote the prosody phase. Considering the 1:2 phase-locking state, the phase difference between theta oscillation and prosody is reasonably defined as 2ϕ_p − ϕ_θ. The phase difference histograms of each participant showed 1:2 phase locking (Fig 6A). Next, we considered the phase oscillator model for 1:2 phase locking as follows: (24) (25) We confirmed that the estimation result can explain the experimental data by calculating the simulated phase difference histograms under the 1:2 phase-locking condition (Fig 6B). The simulated histograms were similar to the experimental histograms (Fig 6C). The results indicate that the estimation phase interaction functions can explain the experimental data, as well as the 1:1 phase-locking condition. The estimated coupling functions of all participants were consistent (Fig 6D and 6E). Furthermore, the estimated coupling functions Γ_p,θ showed small amplitude or were identically zero (Fig 6E). These results clearly show that there were no coupling function from EEG to speech sound, which is reasonable given the relationship between EEG and speech sounds in the experiments. In the case of Γ_p,θ = 0, the simulated phase difference histograms also showed phase locking (Fig 6F), as well as Fig 6B. In contrast, phase locking disappeared for Γ_θ,p = 0 (Fig 6G). These results indicate that the coupling function Γ_θ,p contributed to 1:2 phase locking.

In both cases, our method could estimate whether there was a relationship between the EEG activity and speech sound even though there was some variance due to estimation inaccuracies. In the case of the theta oscillation and prosody data, the estimated coupling functions Γ_p,θ showed small amplitudes or were identically zero, clearly demonstrating the asymmetry of the relationship between the EEG activity and speech sound. In contrast, for the theta oscillation and syllable data, the estimated coupling functions Γ_s,θ showed somewhat larger amplitudes than Γ_p,θ, giving no clear indication of an asymmetric relationship. Therefore, in order to determine whether there was a asymmetry relationship, we estimated the coupling functions using surrogate data (S1 Fig). The surrogate data consisted of randomly time-shifted speech-sound phase data and original EEG data that had not been time-shifted. Note that there was no temporal relationship between the time-series, speech-sound phase data and the original EEG signals. This random shifting process was repeated 100 times for each of the 14 participants. Then, we estimated the coupling functions using these 1,400 surrogate datasets and computed histograms of the model selection results for the appropriate Fourier modes based on the logarithmic evidence and the coupling function powers . For the coupling functions Γ_θ,s and Γ_θ,p, the model selection histograms showed that the M = 0 model was selected more often than the other models (S1A and S1E Fig), while the original data results showed that none of the participants selected the M = 0 model. Furthermore, the integrated values of coupling function power showed that the original data results did not follow the same histograms as the surrogate data results (S1B and S1F Fig). For the coupling functions Γ_s,θ and Γ_p,θ, the model selection results for the original data were not largely different from those that for the surrogate data (S1C and S1G Fig), and the integrated coupling function values were relatively small (S1D and S1H Fig). Consequently, these results suggest that the coupling functions, Γ_s,θ and Γ_p,θ, showed no relationship between the EEGs and speech sounds.

Validity of estimation method for cross-frequency synchronization

We can obtain the theoretical coupling function from the numerical simulation and data from the electronic circuit experiment. To confirm the reliability of the proposed method, we compared the estimated results to theoretical results. In the simulation, time-series data were generated by numerial simulation used in the given phase oscillator model. In this case, we knew the true instantaneous phase of the time-series data and the correct coupling functions. The proposed method worked well with the simulation data (Fig 3), and the complexity parameter for the coupling function M_ij was selected correctly.

We also estimated a dynamical system using electric potential data, i.e., real time-series data. In this situation, the corresponding theoretical coupling function was computed using the adjoint method from which we constructed a theoretical model of real electronic circuits. The results demonstrated that the proposed method can correctly estimate the coupling functions (Fig 4). Furthermore, to confirm that the estimated phase oscillator model can reproduce real time-series data, we compared the phase difference histogram of real time-series data with those of the simulated data. The results demonstrate that the estimated coupling functions and the noise strengths can explain the real data, including any additional disturbance imposed on the system.

In the EEG data, the correct coupling function to be compared to the estimated function is unknown. Therefore, we must consider an alternative procedure to examine the validity of the estimation results. To this end, we considered the following three steps. In the first step, we focused on the coupling functions from the EEG activity to the speech stimulus (Figs 5E and 6E). Under this EEG experimental condition, EEG activity did not influence speech sound because the timing of the external speech sound was given by a recorded sound. Our estimated dynamics showed that the coupling function from EEG to speech sound had smaller amplitude than the coupling function in the opposite direction (Figs 5D and 6D) and did not influence the phase difference histograms (Figs 5F and 6F). Therefore, the estimated network structure is consistent with the real EEG and speech system under this experimental condition. In the second step, we confirmed whether the simulated phase difference histograms were similar to the experimental histograms. Our estimation results and the experimental data both showed phase locking (Figs 5A, 5B, 6A and 6B). In addition to results of the first step, these results suggest that the phase oscillator model can explain the EEG and speech sound data. In the final step, we confirmed that the estimated coupling functions were consistent across all participants to check whether the above results occurred by chance. Our results indicated that the coupling functions and phase difference histograms were similar across all participants on the Cz electrode (Figs 5 and 6). Furthermore, the estimation results for neighbor electrodes (e.g., FCz, Pz, CP, and CP2) showed results similar to those obtained on the Cz electrode. Based on the results obtained by performing these three steps, the mechanism between EEG and speech sound can be explained by the dynamical phase oscillator system.

Remarks on estimation method

To apply the proposed method to EEG data, it is necessary to consider whether the systems to be estimated can be considered a weakly coupled oscillator system. It is well known that EEG phases are often locked by external input among trials [51, 52]. In addition, it is possible that phase locking is generated by phase-resetting through strong external inputs. If phase-resetting depends on the timing of external inputs rather than the phase difference, the interaction between the EEG phase and the stimulus cannot be explained by the phase oscillator model. Therefore, it is necessary to avoid phase resetting caused by sudden and strong external inputs. Event related phase locking is often induced at stimulus onset. Therefore, to prevent such phase-resetting, we presented noise (increased linearly over 0.5 s) prior to presenting speech sound. Furthermore, to avoid a situation where a strong external input induces phase-resetting at speech onset, we employed a bandpass filter to decrease the strong periodic speech sound signals.

Note that the proposed method cannot estimate the coupling function if the EEG phase is completely synchronized. Under such synchronization conditions, each phase difference between the two oscillators takes only a specific value. Therefore, except for this specific value, there is no information about the coupling function on the other phase difference value. To obtain the full range of the coupling function, the phase differences in the data must be distributed in the range 0 to 2π, as shown in Fig 3.

This study focused primarily on 1:p phase synchronization; however, other types of cross-frequency synchronization exist [8, 10], e.g., phase-amplitude, amplitude-amplitude, and phase-frequency synchronization. These synchronizations are also important from a cognitive function perspective; however, the proposed method cannot be applied to such dynamical systems. In future, we plan to construct a method that is applicable to the experimental data of these synchronizations.

Methods to quantify the causality between different frequency rhythms [53–56] have been proposed. Such methods, including transfer entropy and Granger causality, may reveal more general causality than the proposed method, which can only estimate the coupling function related to the p_i:p_j phase synchronization. However, the proposed method can reveal causality and quantify the phase interaction function; thus, it can examine the role of connection in phase synchronization from a dynamical system perspective.

The proposed method can estimate the coupling functions of simulation data and experimental electronic circuit data accurately. In the EEG experiment, we estimated the dynamical system of EEG and speech sound. Note that we examined the dynamics of phase synchronization between a single EEG activity and speech sound rather than between two EEG activities. We applied the proposed method to synchronization between EEG and speech because the direction of causality between EEG and speech sound is clear, whereas that of EEG phases is unknown. Therefore, EEG and speech data were used to verify the estimation results. It is expected that the proposed method can serve as a useful tool to reveal the role of connectivity and causality in neural oscillations.