A mathematical model of dysfunction of the thalamo-cortical loop in schizophrenia

Abstract

Background

Recent experimental results suggest that impairment of auditory information processing in the thalamo-cortical loop is crucially related to schizophrenia. Large differences between schizophrenia patients and healthy controls were found in the cortical EEG signals.

Methods

We derive a phenomenological mathematical model, based on coupled phase oscillators with continuously distributed frequencies to describe the neural activity of the thalamo-cortical loop. We examine the influence of the bidirectional coupling strengths between the thalamic and the cortical area with regard to the phase-locking effects observed in the experiments. We extend this approach to a model consisting of a thalamic area coupled to two cortical areas, each comprising a set of nonidentical phase oscillators. In the investigations of our model, we applied the Ott-Antonsen theory and the Pikovsky-Rosenblum reduction methods to the original system.

Results

The results derived from our mathematical model satisfactorily reproduce the experimental data obtained by EEG measurements. Furthermore, they show that modifying the coupling strength from the thalamic region to a cortical region affects the duration of phase synchronization, while a change in the feedback to the thalamus affects the strength of synchronization in the cortex. In addition, our model provides an explanation in terms of nonlinear dynamics as to why brain waves desynchronize after a given phase reset.

Conclusion

Our model can explain functional differences seen between EEG records of healthy subjects and schizophrenia patients on a system theoretic basis. Because of this and its predictive character, the model may be considered to pave the way towards an early and reliable clinical detection of schizophrenia that is dependent on the interconnections between the thalamic and cortical regions. In particular, the model parameter that describes the strength of this connection can be used for a diagnostic classification of schizophrenia patients.

Introduction

Schizophrenia is a severe and complex mental illness causing disability [1–3]. It has been conceptualized as a disconnectivity syndrome concerning the interplay of the brain areas involved. As information on the activity of some of the deeply localized involved brain areas, such as the thalamus is not accessible to noninvasive electroencephalography (EEG) measurement, alternative methods, like mathematical models, need to be developed in order to deepen our understanding of the fundamental neural processes underlying schizophrenia, and to detect dysfunctions in the interactions between the participating brain areas. Such methods aim at deriving reliable criteria that indicate the progress of the disease at an early stage. The early recognition is considered to be of special importance in schizophrenia.

As schizophrenia has a very high degree of complexity, due to the large number of neuronal processes involved, there is no model that treats all aspects of the disease concurrently. The model in [4], for instance, focuses on the empirical dopamine hypothesis of schizophrenia and postulates that an imbalance between glutamate and dopamine activity plays a key role in schizophrenia disorder. In particular, the authors could show that both hypoglutamatergia and hyperdopaminergia result in reduced activation of the striatal complex and thus leads to schizophrenia [4].

Heiden et al. [5] on the other hand modeled the basic neural circuit underlying schizophrenia as a dynamical system on a microstructural level of pyramidal cells (see also Mackey et al. [6]). In their model, periodic firing patterns were associated with healthy behavior, whereas aperiodic/chaotic firing patterns were associated with schizophrenic states and the switch from periodic to aperiodic firing took place due to an increase in dopamine level. In [7], the same model was analyzed in a hypo-glutamatergic setting, and it also exhibited aperiodic firing.

A top-down approach in modeling the symptoms of schizophrenia is proposed in [8, 9]. These authors relate cognitive, negative and positive symptoms of schizophrenia to a reduced depth of attractor basins of the model regarded as a dynamical system. The model consisting of pools of leaky integrate-and-fire neurons for the involved pyramidal cells and inhibitory interneurons and AMPA, NMDA and GABA _A synapses either developed spontaneous firing or was attracted to a high-firing state. The decrease in NMDA-receptor conductance led to a decrease in attractor stability and therefore a decrease in memory and an increase in distractability [8, 9]. An additional decrease in the conductance of GABA-synapses led to jumping from spontaneous activity to attractors which could be identified with the increase of positive symptoms [8].

In [10], fMRI data of healthy and schizophrenic test persons obtained in a memory-task experiment were used to assess the connectivity between the visual, parietal and prefrontal regions using DCM as proposed by [11]. Each of these regions has a self-coupling and bidirectional all-to-all coupling. Moreover, the working memory modulates the coupling from visual to prefrontal and from prefrontal to parietal regions. The schizophrenic patients and the healthy subjects differed significantly in the self-coupling strengths, in the coupling between parietal and prefrontal regions in both directions and in the coupling from visual to prefrontral regions [10].

In contrast to the modeling studies described above, we focused our attention on the phenomena observed on the large-scale level of dysfunction of the thalamo-cortical loop in schizophrenia. We did therefore not include any explicit biophysical properties into our model. We constructed a mathematical model based on a study by [12]. In this study, the difference between healthy subjects and schizophrenia patients was investigated, using individual EEG recordings with respect to phase locking (PL) in the four frequency bands (θ, α, β and γ). The two groups of test persons showed significant difference in the strength of PL in the θ- and α-frequency band, while no significant difference in PL was observed in the other two frequency bands. The duration of PL, i.e. the time until the system desynchronize again, differed for each frequency band. However it did not differ significantly between healthy subjects and schizophrenia patients in each of the four frequency bands.

In our earlier model [13], we described each component of the thalamo-cortical loop as a single phase oscillator, where each of them operated at its natural frequency. We used phase oscillators for the description of the dynamics in each brain area, since previous studies suggested that the timing of the brain rhythms, i.e. the phases, were more important than their amplitudes [12]. The coupling between these oscillators was expressed in form of weighted phase differences with suitable coefficients, which were determined from the natural frequencies of the oscillators. Using this model, we were able to explain the difference in phase locking in the respective frequency bands of the two groups of test persons depending on the feedback from the cortex to the thalamus. In this model, however, the synchronization effect observed directly after the given auditory stimulus did not abate and vanish after a given period of time, contrary to what was seen in the experiments [12]. We had therefore to change our mathematical model in order to be able to account for the desynchronization effects, too.

In the present study we extended this mathematical model such that now each area of the thalamo-cortical loop is represented by a large population of phase oscillators. The coupling between populations is driven by a complex meanfield (definition see below). To reduce this high-dimensional model to a low dimensional system which still reflects the behavior observed in the EEG data and to allow its analysis, we use the reduction methods of Watanabe- Strogatz [14], Ott-Antonsen [15] and Pikovsky-Rosenblum [16]. The mathematical analysis of the model offers a conducive explanation for the underlying mechanisms leading to the differences observed between healthy subjects and schizophrenia patients, as seen in the experiments by [12]. Our results suggest that the differences are due to a decrease in strength of the coupling from the auditory cortex to the thalamus in schizophrenia patients. Even so, our model is a rather abstract description of the neural dynamics that take place in the thalamo-cortical loop. A decrease in coupling strength can occur due to changes in the dopamine, glutamate or serotonin concentrations. This means in any case a reduction in signal transduction from the auditory cortex to the thalamus. Furthermore, our analysis of the reduced system reveals that the mechanism underlying the abolition of synchrony observed in all four brain wave bands is based on a fold limit cycle bifurcation that takes place when the coupling between the auditory cortex and the thalamus is changed (in either direction). Our model additionally predicts that a change in coupling strength from the thalamus to the auditory cortex, however, affects the duration of phase synchrony.

The paper is organized as follows. In section “The experimental setup and results”, we review the experimental setup and the results of the study by [12]. In section “Mathematical model”, we present the general structure of the thalamo-cortical loop and set up a mathematical model which we use to analyze first the behavior of two coupled brain regions, the thalamus and one cortical region, and then that of three coupled ones, the thalamus and two cortical regions.

The experimental setup and results

The mathematical model presented below is based on experimental results (for details of the experiment and methods used see [12]). In the following, we outline the experiment performed in [12]. Two groups of participants were investigated during the experiment: the first group consisted of 32 schizophrenia patients and the second of 32 healthy subjects. The experiment was based on the well-established paired click paradigm [17]. It consisted of 96 paired clicks (S₁ and S₂). Each click had a duration of 1 ms. The interstimulus interval between the two clicks within a pair lasted 500 ms, and the inter-trial interval between pairs of clicks 10 s. The EEG was continuously recorded using 32 electrodes during the whole experiment. Data from the vertex electrode Cz were taken for the analysis, because the cortex around the location of this electrode performs sensory and motor functions, see [18]. The recorded data have been divided into epochs of 1500 ms (500 ms prior to S₁ and 500 ms following stimulus S₂). The occurence of stimulus S₁ in each segment was set to t=0, hence the stimulus of each segment appeared at t=0. To obtain detailed information on the temporal and spectral properties of the EEG, a single-trial analysis was applied to the epochs. Thus a complex Morlet wavelet transformation in the frequency range from 3 Hz to 60 Hz in 1 Hz steps was performed to compute the phases of the single-trial data. A typical result is displayed in Figure 1 (adapted from [13]) where the cosine of the single trial phase after the wavelet transformation for a fixed frequency 54 Hz is shown. It includes 82 superimposed segments. Uniform distribution of the phases prior stimulus onset, i.e. for t∈[ -50,0] is clearly visible, while the so called phase locking effect after the stimulus, i.e. for t∈[ 0,75], and the effect of desynchronisation after t=75 ms can also be clearly discerned.

Figure 1

Experimental data. The cosine of the phases, derived from EEG data of a single participant, expressing a frequency of 54 Hz is shown (from [13]).

The stimulus locking index

$$\mathit{\text{SLI}}\left(t\right)=\left|\frac{1}{n}\sum _{k=1}^n{e}^{2\pi i{\varphi }_k\left(t\right)}\right|$$

(1)

defined in [19], can be used to measure the degree of phase locking for a certain frequency at time t. We denote the number of repetitions of the auditory double clicks by n and the phase of the k-th oscillation at time t by ϕ _k (t). Values near 0 indicate a uniform distribution of phases and values near 1 nearly synchronized phases. It was found that schizophrenia patients produced significantly less phase locking in lower frequency bands after the first stimulus than healthy subjects [12].

Table 1 shows the SLI and the duration of synchronization for the two groups of subjects for the θ- and α-frequency band for which significant differences were found in the SLI (see [12]). Based on these experimental results, we constructed a mathematical model of the thalamo-cortical loop and used it to explain the observed differences between the neural activities of schizophrenia patients and healthy subjects.

Table 1 Example behavior of cortical regions

Mathematical model

Our model of the thalamo-cortical loop is based on the results of [4, 12, 13, 20]. According to the experimental findings in these sources, we assume that essentially three main brain areas are actively involved in auditory signal processing, i.e. the thalamic auditory relay nucleus (here for the sake of simplicity, they are referred to as thalamus), the thalamic reticular nucleus (here named TRN) and areas of the auditory cortex. An auditory input signal reaches the thalamus and then propagates to the auditory cortex. From the auditory cortex, the signal propagates to higher cerebral regions such as the prefrontal cortex and back to the TRN, which inhibits the thalamus. Furthermore, backpropagation from higher regions such as the prefrontal cortex modulates the activity of the thalamus. These inhibitory and modulating influences lead to a reduced response of the thalamus to the second of the two clicks [21–23]. Since we are only interested in the dynamics after the first and before the second stimulus, we neglected the impact of the TRN. This means that only the thalamus and different regions of the auditory cortex are present in our model. It can be assumed that the different cortical regions act in different oscillatory frequency ranges, which correspond to the θ, α, β and γ ranges. The structure of the thalamo-cortical loop used in our model is shown in Figure 2.

Figure 2

Thalamo-cortical loop. The general structure of the thalamo-cortical loop used in our modeling study. K _T and K _C denote coupling strengths between the thalamus and the auditory cortex. Both the thalamus and the cortex are represented by one population of oscillators each. I denotes the stimulation strength.

In our earlier work [13], we described each part of the thalamo-cortical loop by a single phase oscillator. Each oscillator had a natural frequency, which was chosen according to biological/experimental observations. We assumed, that coupling between all phases is a form of weighted phase difference with suitable coefficients to be determined from the natural frequencies. With this model, it was possible to reproduce the effects that correspond to phase locking as observed in the EEG data of the two groups of subjects. However, in contrast to what the data showed (see Figure 1, t >75 ms), the oscillators remained in the synchronized state perpetually, and their phases did not desynchronize again. To overcome this major drawback of our old model, we now describe each element of the thalamo-cortical loop as a large ensemble of nonidentical phase oscillators. Each oscillator in the population has a natural frequency ω, which is chosen from a Lorentz distribution n(x). The coupling between populations is driven by a complex meanfield. In the experiments, repeated stimulations were used in order to obtain stimulation moments at different phases. In our mathematical model we use 1000 oscillators with distributed phases in each population and stimulate each oscillator at t=0. Since we choose the initial conditions for each oscillator to be different, the results obtained by stimulating them at only one point in time are comparable to the experimental conditions.

First, we will consider a minimal mathematical model, which consists of only two populations of oscillators, one for the thalamus and one for the θ-frequency band of the auditory cortex. We will use this simplified model to understand the mechanism behind the transition from the synchronized to the desynchronized state after stimulation (as seen in Figure 1).

Minimal mathematical model (two populations)

In the minimal model, two populations of oscillators are coupled via their complex mean fields as shown in Figure 3. One of them describes brain wave activity in the thalamus and the other one in the cortex, in this case in the θ-band. In the course of this work, we will refer to these populations as thalamus population and cortex population, respectively.

Figure 3

Mean field coupling. Visualization of the complex mean field coupling between the thalamic region (T) and a single cortical region (C). Black represents the individual oscillators and red their complex mean fields. The oscillators of one population are coupled through the mean phases Θ and the degree of synchronization R of the oscillators in the other population. See text for details.

The system describing the neural activities of the two populations reads as follows

$$\begin{array}{lcr}\frac{d{\varphi }_T\left({\omega }_T\right)}{\mathit{\text{dt}}}& =& {\omega }_T+K_C{R}_{C}sin\left({\theta }_C-{\varphi }_T\left({\omega }_T\right)\right)+I\left(t\right)cos\left({\varphi }_T\right({\omega }_T\left)\right),\end{array}$$

(2)

$$\begin{array}{lcr}\frac{d{\varphi }_C\left({\omega }_C\right)}{\mathit{\text{dt}}}& =& {\omega }_C+K_T{R}_{T}sin({\theta }_T-{\varphi }_C({\omega }_C\left)\right),\end{array}$$

(3)

where ω _T and ω _C are continuous parameters distributed in each ensemble of oscillators as

$$n_a\left({\omega }_a\right)=\frac{2}{\pi (1+4{({\omega }_a-{\tilde{\omega }}_a)}^2)},a=\{T,C\}$$

(4)

and represent the natural frequencies of the oscillators. ϕ _T and ϕ _C denote the phases of the oscillators of the thalamus and cortex populations, respectively. The state of each population can be described by the distribution density W(x,ϕ,t)=n(x)w(x,ϕ,t), with the conditional distribution density of oscillators denoted by w(x,ϕ,t) [16].

Each oscillator in the cortical population is coupled to the complex mean field

$$Y_T=R_T{e}^{i{\theta }_T}=\int n\left(x\right)\underset{-\pi }{\overset{\pi }{\int }}{e}^{i{\varphi }_T}w(x,{\varphi }_T,t)d\varphi dx$$

(5)

of the thalamic population and each oscillator in the thalamic population is coupled to the complex mean field

$$Y_C=R_C{e}^{i{\theta }_C}=\int n\left(x\right)\underset{-\pi }{\overset{\pi }{\int }}{e}^{i{\varphi }_C}w(x,{\varphi }_C,t)d\varphi dx$$

(6)

of the cortical populations (see Figure 3). Coupling strengths are denoted by K _T and K _C , respectively.

The thalamus population is stimulated by an external stimulus that acts directly on it (see Figure 2). This stimulus is represented by the term I(t) cos(ϕ _T ) where

$$I\left(t\right)=\left\{\begin{array}{ll}I& \text{during stimulus},\\ 0& \text{otherwise}.\end{array}\right.$$

(7)

A complete analysis of the model can be performed by means of the Watanabe-Strogatz (WS) theory [14]. With this theory an N-dimensional system of identical oscillators can be reduced to a three-dimensional system with the global variables ρ, ϕ and ψ. Here ρ is the global amplitute, ϕ and ψ are global phases. The original phase variables can be reconstructed from the obtained WS variables by means of the time-dependent transformation. The theory is described in more detail in the Appendix.

Following [14–16, 24], we transform system (2)-(3) using Equation 36 (see Appendix) and obtain a reduced system of WS equations with the new variable ρ _T , ϕ _T , ψ _T and ρ _C , ϕ _C , ψ _C . By additionally introducing $z_a\left({\omega }_a\right)={\rho }_a{e}^{i{\varphi }_a}$ and the phase shift α _a (ω _a )=ϕ _a -ψ _a , (a=T,C), we obtain

$$\begin{array}{lcr}\frac{{\mathit{\text{dz}}}_T\left({\omega }_T\right)}{\mathit{\text{dt}}}& =& i{\omega }_T{z}_T\left({\omega }_T\right)+\frac{K_C}{2}\left(Y_C-Y_C^{\ast }{z}_T{\left({\omega }_T\right)}^2\right)+\frac{I\left(t\right)i}{2}\left(1+z_T{\left({\omega }_T\right)}^2\right),\end{array}$$

(8)

$$\begin{array}{lcr}\frac{d{\alpha }_T\left({\omega }_T\right)}{\mathit{\text{dt}}}& =& {\omega }_T+\text{Im}\left(z_T{\left({\omega }_T\right)}^{\ast }\left(K_C{Y}_C+I\left(t\right)i\right)\right),\end{array}$$

(9)

$$\begin{array}{lcr}\frac{{\mathit{\text{dz}}}_C\left({\omega }_C\right)}{\mathit{\text{dt}}}& =& i{\omega }_C{z}_C\left({\omega }_C\right)+\frac{K_T}{2}\left(Y_T-Y_T^{\ast }{z}_C{\left({\omega }_C\right)}^2\right),\end{array}$$

(10)

$$\begin{array}{lcr}\frac{d{\alpha }_C\left({\omega }_C\right)}{\mathit{\text{dt}}}& =& {\omega }_C+\text{Im}\left(z_C{\left({\omega }_C\right)}^{\ast }{K}_T{Y}_T\right).\end{array}$$

(11)

Here, $i=\sqrt{-1}$ and A^∗ denotes the conjugate complex of A.

Now we consider this reduced set of equations with respect to the Ott-Antonsen manifold [25]. In this case, z(ω) no longer depends on α(ω), and the mean fields Y _C and Y _T can therefore be written as

$$Y_a=R_a{e}^{i{\theta }_a}=\int n_a\left(x\right)z_a\left(x\right)\mathit{\text{dx}},(a=T,C).$$

(12)

Following the work by Ott and Antonsen for a similar distribution [25], the integrals in Equation (12) can be calculated by applying the residue theorem, under an additional assumption that z _a (ω) (a=T,C) is analytic in the upper half-plane. This calculation yields

$$\begin{array}{lcr}{Y}_T& =& z_T({\tilde{\omega }}_T+i/2),\end{array}$$

(13)

$$\begin{array}{lcr}{Y}_C& =& z_C({\tilde{\omega }}_C+i/2).\end{array}$$

(14)

Thus Equation (8) for ${\omega }_T={\tilde{\omega }}_T+i/2$ and Equation (10) for ${\omega }_C={\tilde{\omega }}_C+i/2$ provide a 2-dimensional system of complex ODEs that describe the behavior of the order parameter of the thalamic and the cortex population, respectively:

$$\begin{array}{l}\frac{{\mathit{\text{dY}}}_T}{\mathit{\text{dt}}}=\left({\tilde{\omega }}_{T}i-\frac{1}{2}\right)Y_T+\frac{1}{2}\left(K_C{Y}_C+I\left(t\right)i-Y_T^2\left(K_C{Y}_C^{\ast }-I\left(t\right)i\right)\right)\end{array}$$

(15)

$$\begin{array}{l}\frac{{\mathit{\text{dY}}}_C}{\mathit{\text{dt}}}=\left({\tilde{\omega }}_{C}i-\frac{1}{2}\right)Y_C+\frac{K_T}{2}\left(Y_T-Y_C^2{Y}_T^{\ast }\right)\end{array}$$

(16)

In the following, we will investigate this system of two complex differential equations, i.e. its dynamics during the post-stimulus interval.

Analysis of the model behavior in the post-stimulation interval

For the analysis of eqs. (15)-(16) in the post-stimulus interval, i.e. when I(t)=0, we transform them to a 4-dimensional system of real ODEs via Y _T =x _T +i y _T and Y _C =x _C +i y _C . This leads to

$$\begin{array}{lcr}\frac{{\mathit{\text{dx}}}_T}{\mathit{\text{dt}}}& =& -\frac{x_T}{2}-{\omega }_T{y}_T+\frac{K_C}{2}\left(x_C-\left(x_T^2-y_T^2\right)x_C-2x_T{y}_T{y}_C\right)\end{array}$$

(17)

$$\begin{array}{lcr}\frac{{\mathit{\text{dy}}}_T}{\mathit{\text{dt}}}& =& {\omega }_T{x}_T-\frac{1}{2}{y}_T+\frac{1}{2}{K}_C\left(y_C-2x_T{y}_T{x}_C+\left(x_T^2-y_T^2\right)y_C\right)\end{array}$$

(18)

$$\begin{array}{lcr}\frac{{\mathit{\text{dx}}}_C}{\mathit{\text{dt}}}& =& -\frac{1}{2}{x}_C-{\omega }_C{y}_C+\frac{1}{2}{K}_T\left(x_T-\left(x_C^2-y_C^2\right)x_T-2x_C{y}_C{y}_T\right)\end{array}$$

(19)

$$\begin{array}{lcr}\frac{{\mathit{\text{dy}}}_C}{\mathit{\text{dt}}}& =& {\omega }_C{x}_C-\frac{1}{2}{y}_C+\frac{1}{2}{K}_T\left(y_T-2x_C{y}_C{x}_T+\left(x_C^2-y_C^2\right)y_T\right).\end{array}$$

(20)

In a next step, we linearize this new system about its fixed point x _F =(0,0,0,0) and investigate the stability of this fixed point with the coupling strengths K _C and K _T as parameters. The linearized system reads:

$$\begin{array}{lcr}\left(\begin{array}{l}\frac{{\mathit{\text{dx}}}_T}{\mathit{\text{dt}}}\\ \frac{{\mathit{\text{dy}}}_T}{\mathit{\text{dt}}}\\ \frac{{\mathit{\text{dx}}}_C}{\mathit{\text{dt}}}\\ \frac{{\mathit{\text{dy}}}_C}{\mathit{\text{dt}}}\end{array}\right)& =& A\left(\begin{array}{l}{x}_T\\ y_T\\ x_C\\ y_C\end{array}\right)\text{with}\end{array}$$

(21)

$$\begin{array}{lcr}A& =& \left(\begin{array}{llll}-\frac{1}{2}& -{\omega }_T& \frac{K_C}{2}& 0\\ {\omega }_T& -\frac{1}{2}& 0& \frac{K_C}{2}\\ \frac{K_T}{2}& 0& -\frac{1}{2}& -{\omega }_C\\ 0& \frac{K_T}{2}& {\omega }_C& -\frac{1}{2}\end{array}\right).\end{array}$$

(22)

A one-dimensional bifurcation diagram is displayed in Figure 4. For the calculation of this diagram, we fix one of the coupling parameters, K _C =1.2, and show the dependence of one of the system variables (x _C ) on the second coupling parameter K _T . The system has a fixed point x _F =(0,0,0,0), which is stable for $K_T<K_T^{\mathit{\text{HB}}}$. At $K_T=K_T^{\mathit{\text{HB}}}=\frac{{({\omega }_T-{\omega }_C)}^2+1}{K_C}$ a Hopf bifurcation (HB) occurs in the system, i.e. a complex conjugate pair of eigenvalues of A passes through the imaginary axis (see Figure 4). At this point, the branch of stable fixed points (x _F , red line in Figure 4) loses its stability because it collapses into a branch of unstable periodic orbits (x _PU , blue circles). Additionally, the system exibits a fold limit cycle bifurcation (FLB) at $K_T^{\mathit{\text{FLB}}}<K_T^{\mathit{\text{HB}}}$. At this bifurcation point, two periodic orbits, a stable (x _PS , green discs) and an unstable one (blue circles) are born. The bifurcation diagram reveals three parameter regions in which the system displays different behavior. In region I $\left(0<K_T<K_T^{\mathit{\text{FLB}}}\right)$ we have a stable fixed point x _F =(0,0,0,0), which corresponds to the state of full desynchronization in the non-reduced system (2)-(3). In region III $\left(K_T>K_T^{\mathit{\text{HB}}}\right)$ the fixed point x _F has lost its stability and all trajectories are attracted to the stable periodic orbit x _PS (filled green circles in Figure 4). This corresponds to a state near perfect synchronization of the non-reduced system. In region II $\left(K_T^{\mathit{\text{FLB}}}<K_T<K_T^{\mathit{\text{HB}}}\right)$, the system is bistable. It can exhibit fixed point solutions as well as periodic ones. Both behaviors are separated by an unstable periodic orbit x _PU . Depending on the initial conditions of the system, the trajectory will stay in the region of attraction of the fixed point x _F or is attracted by the stable periodic orbit x _PS .

Figure 4

Bifurcation diagram. Bifurcation diagram in the (K _T -x _C )-plane. The coupling strength from the cortex to the thalamus has been set to K _C =1.2 and the coupling strength K _T from the thalamus to the cortex is varied; ${\tilde{\omega }}_C=3$ Hz and ${\tilde{\omega }}_T=7$ Hz. Red line: stable fixed points (x _F ), black line: unstable fixed points; green filled circles represent the maximum and minimum of stable periodic orbits (x _PS ) while the blue circles represent the maximum and minimum of unstable periodic orbits (x _PU ). Bifurcation points, i.e. the fold limit cycle bifurcation (FLB) at $K_T^{\mathit{\text{FLB}}}$ and the Hopf-bifurcation (HB) at $K_T^{\mathit{\text{HB}}}$ are marked with vertical blue lines. The light blue box represents the parameter region where the initially synchronized system desynchronizes in finite time (like in Figure 1).

Let us now focus on the blue region surrounding the bifurcation points FLB. It is possible to choose a value of the parameter K _T inside this region such that the trajectory is resetted to a state near the periodic orbit and drops back to the stable fixed point after a certain amount of time.

Figure 5 shows the behavior of the reduced system (left) and the corresponding behavior of the non-reduced system (right) for the three parameter regions described above. Red lines in the figures of the left column indicate the maximum amplitude of the solutions. Black thick lines in the figures of the right column indicate the SLIs. In region I (see Figure 4), i.e. for low value of K _T , e.g. K _T =1, the reduced system exhibits a stable fix point solution (Figure 5(a), left), and the corresponding non-reduced system is in the desynchronization regime (Figure 5(a), right), hence SLI = 0. For K _T in region III (see Figure 4), e.g. K _T =16, the stable periodic orbit of the reduced system is shown in Figure 5(b) on the left and the corresponding synchronization regime of the non-reduced system on the right with SLI = 0.6. In Figure 5(c) K _T is fixed near $K_T^{\mathit{\text{FLB}}}$ (K _T =5.5, left edge of blue region in Figure 4). With this parameter choice, we observe the same dynamics as seen in the experiments: before stimulation we have desynchronization in the non-reduced system and a fix point in the reduced one. After the stimulation interval (marked with vertical dash lines), we see a phase reset, and the phases of the oscillators of the non-reduced system are now synchronized (thin wave). After some time (t >700 ms) they desynchronize again. The length of the synchronization state can be modulated by changing the distance of K _T to $K_T^{\mathit{\text{FLB}}}$. The closer K _T is set to $K_T^{\mathit{\text{FLB}}}$ the longer the trajectory will stay in the state of synchronization before it desynchronizes again.

Figure 5

Bifurcation behaviors. The real part of the solution of the reduced system (left) and the corresponding cosine of the solution of the non-reduced system (right) for (a) K _T =1 in region I, (b) K _T =16 in region III and (c) K _T =5.5 in region II of Figure 4 are shown. In (c) the behavior of both systems is shown during and after stimulation. For the reduced system, red lines indicate the maximum amplitude of the solutions (panels in the left column). For the non-reduced system, the oscillators’ (only 200 out of N=1000 displayed) activities are shown as blue curves, and the corresponding max SLI is plotted in black (panels in the right column). K _C =1.2, ${\tilde{\omega }}_C=3$ Hz, ${\tilde{\omega }}_T=7$ Hz. The stimulus intensity was set to I=100 and had a duration of 50 ms for all simuli.

In Figure 6, top left, the maximum SLI of the cortex population for different pairs of coupling strengths (K _C ,K _T ) is shown. For this we stimulated the system of two populations (N=1000 oscillators each) with a stimulation strength of I=100 for 50 ms and calculated $\underset{t}{max}\frac{1}{1000}\sum _{j=1}^{1000}{e}^{2\pi i{\phi }_{C_j}\left(t\right)}$ for each fixed pair (K _C ,K _T ). The dependence of the peak of the frequency distribution on these parameters was calculated using the fast Fourier transformation (MATLAB function FFT). It is shown in the bottom left panel of Figure 6. The coupling strengths K _T and K _C were varied independently from 0 to 10 with a step size of 0.05. The grey and black curves in the K _C - K _T plane represent the branches of fold limit cycle bifurcations (FLBs) of the periodic orbits and of Hopf-bifurcations (HBs) of the fixed points, respectively. As seen in Figure 4, only parameter values to the left and close to the grey curve and parameter values which yield S L I≥0.3 guarantee a strong synchronized population which desynchronizes in finite time. The closer the coupling parameters are set to the branch of FLBs the longer the population will stay in the synchronized state. When the value of the coupling parameter K _T is decreased the system moves away from the curve of FLBs, thus the synchronization of the system becomes weaker and shorter. This happens even for small changes in K _T (direction denoted by B in Figure 6, top left). Changing the value of the coupling parameter K _C , however, has a much weaker effect on the strength of synchronization (direction denoted by A in Figure 6, top left). An enlargement of this region of interest is shown in Figure 6, top right.

Figure 6

Coupling parameter dependence. Left: Maximal SLI (top) and mean frequency (bottom) of the oscillators in the cortical region (for ${\tilde{\omega }}_C=3$ Hz, ${\tilde{\omega }}_T=7$ Hz and N=1000) as functions of the coupling strengths K _C and K _T . Right: enlargments of the marked areas in the corresponding panels. Solid curves mark Hopf-Bifurcations (black) and fold-limit cycle bifurcations (gray) of the reduced system. The arrows A and B show possible transitions from a level of high synchronization (K _C =2,K _T =5) to a level of lower synchronization (tip of arrow A: K _C =1.1,K _T =5; tip of arrow B: K _C =2,K _T =4.5).

We know from the experimental data of [12] that schizophrenia patients show a lower synchronization in the θ-band than healthy subjects (max SLI = 0.3 and 0.37, respectively). If we now change the parameters K _C and K _T such that the SLI changes from 0.37 to 0.3 (in the A or B direction or in a direction representing a linear combination of the two), the mean frequency of the θ-oscillators does only minimally change, i.e. it remains approximately 6 Hz. A major change in the mean frequency should of course not happen when the conditions of the system are changed from the healthy to the schizophrenic state. Figure 6, bottom left shows the mean frequency of the cortical oscillators for different pairs of coupling parameters (K _C , K _T ). An enlargement of the region of interest is shown in Figure 6, bottom right.

Summing up, our model simulations nicely show that the neural dynamics observed in EEG data of schizophrenia patients and healthy subjects strongly depend on the strength of the coupling between the thalamus and the cortex: decreasing one or both of the coupling parameters K _T or K _C in an appropriate manner decreases the max SLI of the system and changes the time the system stays synchronized but leaves its mean frequency nearly unchanged.

Model with three populations

In this section, we consider the extension of our minimal model to three populations of coupled phase oscillators, where one of them describes brain activity in the thalamus and the other two θ- and α-brain waves in the auditory cortex. We only consider the θ- and α-frequency band here since the experimental data of [12] showed significant differences between schizophrenia patients and healthy subjects in these frequency bands, only.

Our extended system has the form

$$\begin{array}{lcr}\frac{d{\varphi }_T\left({\omega }_T\right)}{\mathit{\text{dt}}}& =& {\omega }_T+\sum _{j=1}^2{K}_{C_j}{R}_{C_j}sin\left({\theta }_{C_j}-{\varphi }_T\left({\omega }_T\right)\right)+I\left(t\right)cos\left({\varphi }_T\left({\omega }_T\right)\right),\end{array}$$

(23)

$$\begin{array}{lcr}\frac{d{\varphi }_{C_1}\left({\omega }_{C_1}\right)}{\mathit{\text{dt}}}& =& {\omega }_{C_1}+K_{T_1}{R}_{T}sin\left({\theta }_T-{\varphi }_{C_1}\left({\omega }_{C_1}\right)\right),\end{array}$$

(24)

$$\begin{array}{lcr}\frac{d{\varphi }_{C_2}\left({\omega }_{C_2}\right)}{\mathit{\text{dt}}}& =& {\omega }_{C_2}+K_{T_2}{R}_{T}sin\left({\theta }_T-{\varphi }_{C_2}\left({\omega }_{C_2}\right)\right).\end{array}$$

(25)

All notations are the same as in the minimal model. ϕ _T and ${\varphi }_{C_1}$, ${\varphi }_{C_2}$ are the phases corresponding to the oscillators of the thalamus population and the oscillators of the θ cortex and α cortex populations, respectively. ω _T , ${\omega }_{C_1}$ and ${\omega }_{C_2}$ represent the natural frequencies of the corresponding oscillators chosen from the distribution described in Equation 4.

Each oscillator of both cortical populations is coupled to the complex mean field

$$Y_T=R_T{e}^{i{\theta }_T}=\int n\left(x\right)\underset{-\pi }{\overset{\pi }{\int }}{e}^{i{\varphi }_T}w(x,{\varphi }_T,t)d\varphi dx$$

of the thalamic population with the corresponding coupling parameters $K_{T_j}$, j=1,2. Each oscillator in the thalamic population is coupled to both complex mean fields

$$Y_{C_j}=R_{C_j}{e}^{i{\theta }_{C_j}}=\int n\left(x\right)\underset{-\pi }{\overset{\pi }{\int }}{e}^{i{\varphi }_{C_j}}w(x,{\varphi }_{C_j},t)d\varphi dx,j=1,2$$

of the θ and α cortical populations with the corresponding coupling parameter $K_{C_j}$, j=1,2.For the sake of simplicity and as a first approximation, we assume, that there is no direct connection between the two cortical populations. Note, that these populations are considered as functionally distinct groups in the cortex, not as anatomically distinct ones. They, however, influence each other indirectly through the feedback they receive from the thalamic population. Figure 7 shows a schematic illustration of the thalamo-cortical loop with two cortex populations. The thalamus population is stimulated with the same external stimulus as in the case of the minimal model (see Equation 7).

Figure 7

Three population structure. Scheme of the three-population model. Fixed couplings are shown in black and variable ones in red.

Again, we apply the Watanabe-Strogatz ansatz [14] and its extension by Pikovsky and Rosenblum [16] to the system (23)-(25) and obtain

$$\begin{array}{lcr}\frac{{\mathit{\text{dz}}}_T\left({\omega }_T\right)}{\mathit{\text{dt}}}& =& i{\omega }_T{z}_T\left({\omega }_T\right)+\sum _{j=1}^2\frac{K_{C_j}}{2}\left(Y_{C_j}-Y_{C_j}^{\ast }{z}_T{\left({\omega }_T\right)}^2\right)+\frac{I\left(t\right)i}{2}\left(1+z_T{\left({\omega }_T\right)}^2\right),\end{array}$$

(26)

$$\begin{array}{lcr}\frac{d{\alpha }_T\left({\omega }_T\right)}{\mathit{\text{dt}}}& =& {\omega }_T+Im\left(z_T{\left({\omega }_T\right)}^{\ast }\left(K_{C_1}{Y}_{C_1}+K_{C_2}{Y}_{C_2}+I\left(t\right)i\right)\right),\end{array}$$

(27)

$$\begin{array}{lcr}\frac{{\mathit{\text{dz}}}_{C_1}\left({\omega }_{C_1}\right)}{\mathit{\text{dt}}}& =& i{\omega }_{C_1}{z}_{C_1}\left({\omega }_{C_1}\right)+\frac{K_{T_1}}{2}\left(Y_T-Y_T^{\ast }{z}_{C_1}{\left({\omega }_{C_1}\right)}^2\right),\end{array}$$

(28)

$$\begin{array}{lcr}\frac{d{\alpha }_{C_1}\left({\omega }_{C_1}\right)}{\mathit{\text{dt}}}& =& {\omega }_{C_1}+Im\left(z_{C_1}{\left({\omega }_{C_1}\right)}^{\ast }{K}_{T_1}{Y}_T\right),\end{array}$$

(29)

$$\begin{array}{lcr}\frac{{\mathit{\text{dz}}}_{C_2}\left({\omega }_{C_2}\right)}{\mathit{\text{dt}}}& =& i{\omega }_{C_2}{z}_{C_2}\left({\omega }_{C_2}\right)+\frac{K_{T_2}}{2}\left(Y_T-Y_T^{\ast }{z}_{C_2}{\left({\omega }_{C_2}\right)}^2\right),\end{array}$$

(30)

$$\begin{array}{lcr}\frac{d{\alpha }_{C_2}\left({\omega }_{C_2}\right)}{\mathit{\text{dt}}}& =& {\omega }_{C_2}+Im\left(z_{C_2}{\left({\omega }_{C_2}\right)}^{\ast }{K}_{T_2}{Y}_T\right).\end{array}$$

(31)

This system can again be reduced to a 3-dimensional system of complex ODEs representing the dynamics of the order parameters. The reduced system has the form:

$$\begin{array}{lcr}\frac{{\mathit{\text{dY}}}_T}{\mathit{\text{dt}}}& =& \left({\tilde{\omega }}_{T}i-\frac{1}{2}\right)Y_T+\frac{1}{2}\left(I\left(t\right)i+Y_T^{2}I\left(t\right)i\right)+\\ +\frac{1}{2}\sum _{j=1}^2\left(K_{C_j}{Y}_{C_j}-Y_T^2{K}_{C_j}{Y}_{C_j}^{\ast }\right),\end{array}$$

(32)

$$\begin{array}{lcr}\frac{{\mathit{\text{dY}}}_{C_1}}{\mathit{\text{dt}}}& =& \left({\tilde{\omega }}_{C_1}i-\frac{1}{2}\right)Y_{C_1}+\frac{K_{T_1}}{2}\left(Y_T-Y_{C_1}^2{Y}_T^{\ast }\right),\end{array}$$

(33)

$$\begin{array}{lcr}\frac{{\mathit{\text{dY}}}_{C_2}}{\mathit{\text{dt}}}& =& \left({\tilde{\omega }}_{C_2}i-\frac{1}{2}\right)Y_{C_2}+\frac{K_{T_2}}{2}\left(Y_T-Y_{C_2}^2{Y}_T^{\ast }\right).\end{array}$$

(34)

The main results of the analysis of this model are illustrated in Figure 8. The panels of the figure on the left show the maximum SLI and the frequency for both cortex populations over the plane of the coupling parameters $K_{C_1}$ and $K_{C_2}$. The right column shows enlargements of a region of interest of each of the figures shown in the left column. In these simulations we fixed the coupling from the thalamus population to the θ- and α-populations, i.e. $K_{T_1}=5.5$ and $K_{T_2}=7$, respectively. For each parameter pair $(K_{C_2},K_{C_1})$, we stimulated the system (N=1000 oscillators in each of the three populations) with a stimulus intensity of I=100 for 50 ms and calculated the maximum SLI as $\underset{t}{max}\frac{1}{1000}\sum _{j=1}^{1000}{e}^{2\pi i{\phi }_{C_j}\left(t\right)}$. The results are shown in the first and third row of Figure 8. The figures in the second and fourth row of Figure 8 show the mean frequency of the θ- and α-populations, respectively, as functions of the coupling strengths $K_{C_1}$ and $K_{C_2}$. The coupling strengths $K_{C_1}$ and $K_{C_2}$ were varied independently from 0 to 10 with a step size of 0.05. The grey and black curves again denote the branches of fold limit cycle bifurcations (FLBs) of the periodic orbits and of Hopf bifurcations (HBs) of the fixed points, respectively.

Figure 8

Coupling parameter dependence for three populations. Left: maximum SLI and frequency of the θ-population (${\tilde{\omega }}_{C_1}=3$ Hz; top two) and the α-population (${\tilde{\omega }}_{C_1}=14$ Hz; bottom two) as functions of the coupling parameters $K_{C_2}$ and $K_{C_1}$. Right: enlargement of the marked areas. The points in the corresponding left panels C ($K_{C_2}=3.9,K_{C_1}=1.2$) and P ($K_{C_2}=3.2,K_{C_1}=1.2$) mark parameter choices that yield the max SLIs and frequencies observed in healthy subjects (controls) and schizophrenia patients, repectively. The coupling parameters from the thalamus (${\tilde{\omega }}_T=7$ Hz) to the cortex are fixed at $K_{T_1}=5.5$ and $K_{T_2}=7.0$. Solid curves mark branches of Hopf bifurcations (black) and of fold limit cycle bifurcations (gray) of the reduced system.

Using these simulations, we can now navigate through the $(K_{C_2},K_{C_1})$ parameter plane and find pairs of coupling parameter values at which our model exhibits brain dynamics as observed in healthy subjects or in schizophrenia patients. To change the dynamics of the model from the healthy to the schizophrenic state, the parameter pairs need to be chosen such that: i) the max SLI of the θ-brain waves changes from 0.37 to 0.3; ii) the one of the α-waves from 0.26 to 0.19; iii) the frequencies of the θ- and α-populations change only minimally, i.e. they remain approximately 5-7 Hz for the θ- and 9-12 Hz for the α-population.

With the letters P and C, we label the positions in the $(K_{C_2},K_{C_1})$ parameter plane which correspond to the max SLIs and frequencies observed in schizophrenia patients and healthy subjects, respectively. For healthy subjects the coupling constants are $K_{C_2}=3.9,K_{C_1}=2.2$, for schizophrenia patients $K_{C_2}=3.2,K_{C_1}=1.2$. We can see that schizophrenia patients, compared to healthy subjects, have a reduced feedback from both cortex populations to the thalamus.

We therefore hypothesize based on our model that schizophrenia patients have deficits in signal transduction from the auditory cortex back to the thalamus.

Figure 9 shows the behaviour of the θ- and α-populations in the case of schizophrenia patients, i.e. the coupling parameters are $K_{T_1}=5.5$, $K_{T_2}=7$, $K_{C_2}=3.2,K_{C_1}=1.2$. The cosine of the phases of 200 oscillators is shown (in total we calculated N=1000 oscillators in each population). The system was stimulated at t=0. The black lines indicate the max SLI. Before the stimulation, i.e. t∈[ -100,0], we have a uniform distribution of the phases, which means that the oscillators in each population are desynchronized, i.e. max SLI = 0. Directly after the stimulation at t=0 a phase reset occurs and the phases of the oscillators synchronize (thin blue waves), hence now max SLI >0. After a certain time, they desynchronize again: the desynchronization in the α-population happens earlier (at ≈250 ms) followed by that in the θ-population (at ≈400 ms). These simulation results fully agree with the experimental data [12].

Figure 9

Dynamic behavior of cortical regions. An example of the desynchronization in the two cortical areas of the non-reduced system (only 200 of N=1000 oscillators are shown); ${\tilde{\omega }}_{C_1}=3$ Hz ${\tilde{\omega }}_{C_2}=14$ Hz, ${\tilde{\omega }}_T=7$ Hz, $K_{T_1}=5.5$, $K_{T_2}=7.0$, $K_{C_2}=3.2$, $K_{C_1}=1.2$. The max SLI is shown in black. Frequency of coupled oscillators in the phase-locking interval: 6 Hz for θ-population and 9.4 Hz for α-population.

Discussion

We have constructed a mathematical model which describes some aspects of the dysfunction in neural activity of the thalamo-cortical loop during schizophrenia. This model is, in contrast to the models introduced at the beginning [5, 7, 8, 10], an abstract description of the neural activity of the brain regions involved. The description is based on synchronization phenomena as found in the EEG records of healthy subjects and schizophrenia patients. The model allows us to study the interaction of the brain areas involved in a systematic way and to detect the effects of changes of the couplings between them.

A major advantage of our model is that it provides insight into activities of brain areas not directly accessible by EEG measurements. A second advantage is its flexibility: since we did not include any explicit biophysical properties into our model and stayed on a phenomenological macrostructural level, it can be adjusted to fit similar effects in other diseases, e.g. Morbus Parkinson [5] or bipolar disorder [26]. Patients suffering from Morbus Parkinson show a decrease in post-movement synchronization in the β-band [27] in the motor cortex. To simulate the system dynamics underlying this disease, it is possible to extend our model by including oscillators with β-frequencies in the cortex population (as done by the extension of our model from initially 2 to 3 populations). In a graph-theoretical analysis, Kim et al. [26] could show that only synchronization in α frequencies is significantly lower in patients with bipolar disorder than in healthy controls. Thus it would be possible to adjust the coupling strengths in our model in a way that it could appropriately reproduce the symptoms of this disease.

We have seen that the description of each area by a population of phase oscillators with distributed frequencies instead of only a single one per region allows the model to exhibit phase desynchronization in addition to phase-resetting after a given auditory stimulus. With this model it is now possible to get insights into the brain dynamics of the thalamic region which are not accessible to EEG measurements and to investigate the impact of different connection topologies between the thalamic and cortical brain regions on the duration and strength of synchronization in the respective brain frequency bands.

By analyzing our model, we have seen that it supports the current view [9] that the coupling between the thalamic and the cortical regions is responsible for dysfunction of the thalamo-cortical loop in schizophrenia. In particular, our model shows that a reduction in the strength of coupling from the thalamic to a cortical region, i.e. decreasing the coupling parameters $K_{T_i}$, shortens the time interval of phase synchronization after stimulus. The closer $K_{T_i}$ is to the fold limit cycle bifurcation the longer the state of high synchrony lasts. A reduction in the strength of the coupling from the cortical regions back to the thalamus, i.e. decreasing the coupling parameters $K_{C_i}$, however reduces the strength of synchronization. Comparing this with the findings in [12] yields that patients differ mostly in the θ- and α-frequency bands in that they have a significantly lower coupling $K_{C_i}$ from the cortex to the thalamus and only a little weaker coupling $K_{T_i}$ from the thalamus to the cortex, than healthy subjects.

The bifurcation analysis of our model was only possible, because we were able to reduce this large-scale model to systems of only two or three dimensions by using the Ott-Antonsen theory and the reduction methods by Pikovsky and Rosenblum. Thanks to these reduction methods, we could do so without loosing the desynchronization phenomena observed in experiments and exhibited by the large-scale model. The bifurcation analysis helped us to understand the mechanisms underlying the desynchronization which follows the initial phase reset after the auditory stimulus. As far as we are aware of, our model can show for the first time which bifurcations underlie the changes in the simulated brain dynamics.We have, furthermore, seen that the mechanism of desynchronization is preserved if the model is extended from initially 2 to 3 regions to describe the neural activity in the thalamo-cortical loop. In this model, however, both parameter regions which correspond to data of schizophrenia patients and of healthy subjects are rather small (see little island in Figure 8). For the sake of simplicity and as a first approximation, we did not include any direct couplings between the cortical regions. Preliminary results (simulation results not shown) suggest that the size of this parameter region will increase, if additional intra-cortical connections are included in the model, and the mechanisms of desynchronization will still be preserved. However, further investigations will be needed to gather firm evidence that supports the preliminary results.

Simulation results obtained with our model support the notion that schizophrenia is not caused by focal brain abnormalities, but results from pathological interactions between brain regions [28]. In contrast to other studies that found abnormal functional connectivity between temporal and frontal regions, as measured by PET and fMRI [29–31], our model hypothesizes that it is the feedback from the thalamus to the auditory cortex that causes the disabilities, even in the absense of direct coupling between cortical frequencies.

Of course, the validity of our model needs to be further tested by comparing the simulation results with experimental ones obtained in new measurements. The “disconnection hypothesis” suggests that the core pathology of schizophrenia is an impaired neuromodulation of synaptic plasticity that leads to abnormal functional integration of neural systems, i.e., “dysconnectivity” [32, 33]. A possible next step would therefore be to investigate the influence of medications like eicosapentaenoic acid [34], pregnenolone [35, 36] or antipsychotics [37] on the functional signal transduction by calculating individual SLIs from the EEG data or in a DCM study [38]. This would provide information on the coupling between the thalamus and the auditory cortex in schizophrenia patients. One could then compare whether a regeneration of functions can be linked to a restored connection strength between thalamus and cortex.

Moreover, EEG measurements on individual patients and the calculation of their individual coupling strength between the thalamus and the cortex could then be used as a tool to diagnostically classify the different types of schizophrenia. The advantage of this approach would be that the grouping of patients would not simply be driven by data but would be constrained by a well-founded and carefully specified theory.

In order to explain the experimental observations we consider in this paper a simple and plausible model which inherits properties of the coupling of the Kuramoto model. Despite the simplicity of the model, we observed good agreement between experimental data and numerical simulations. This justifies our chosen degree of approximation. Of course more sophisticated phase models including the nonresonance case [39, 40] as well as detailed networks of the brain regions of interest could improve the quality of that results. Particularly, if we extended our model to include other brain regions with distinct natural frequencies, such as β- or γ-bands, we should definitely consider models that incorporate the nonresonance property.

Conclusion

Our model can explain functional differences seen between EEG records of healthy subjects and schizophrenia patients on a system theoretic basis. Because of this and its predictive character, the model may be considered to pave the way towards an early and reliable clinical detection of schizophrenia that is dependent on the interconnections between the thalamic and cortical regions. In particular, the model parameter that describes the strength of this connection can be used for a diagnostic classification of schizophrenia patients.

Appendix

Mathematical theory

In this section, we will present the main reduction methods that are used to analyze our systems in section “Mathematical model”. The Pikovsky-Rosenblum ansatz [16] is an extension of the reduction method by Watanabe and Strogatz [14] which covers infinitely large systems of nonidentical phase oscillators. The general model treated in this theory reads:

$$\frac{d\varphi (x,t)}{\mathit{\text{dt}}}=\omega (x,t)+\text{Im}\left[H(x,t)e^{-i\varphi }\right].$$

(35)

The natural frequencies of the oscillators are denoted by ω(x,t). They depend on a continuous parameter x. The oscillators are coupled via a complex field H(x,t). The state of the system can be described by the distribution density W(x,ϕ,t) which is determined by the distribution density of the parameter n(x) and the conditional distribution density of oscillators w(x,ϕ,t), written as

$$W(x,\varphi ,t)=n\left(x\right)w(x,\varphi ,t).$$

According to the idea presented in [14], three new variables ρ(x,t), Φ(x,t), Ψ(x,t) and constants of motion ψ(x) are introduced to the original system of equations via the transformation

$$e^{i\varphi }=e^{\mathrm{i\Phi }}\frac{\rho +e^{i(\psi -\Psi )}}{\rho e^{i(\psi -\Psi )}+1},$$

(36)

which transforms the time-dependent density w(x,ϕ,t) to a stationary density σ(x,ψ) with the new variables Φ(x,t), Ψ(x,t) and ρ(x,t) which satisfy the following

$$\begin{array}{lcr}\frac{\mathrm{\partial \rho }(x,t)}{\mathrm{\partial t}}& =& \frac{1-{\rho }^2}{2}\text{Re}\left(H(x,t)e^{-\mathrm{i\Phi }}\right),\end{array}$$

(37)

$$\begin{array}{lcr}\frac{\mathrm{\partial \Phi }(x,t)}{\mathrm{\partial t}}& =& \omega (x,t)+\frac{1+{\rho }^2}{2\rho }\text{Im}\left(H(x,t)e^{-\mathrm{i\Phi }}\right),\end{array}$$

(38)

$$\begin{array}{lcr}\frac{\mathrm{\partial \Psi }(x,t)}{\mathrm{\partial t}}& =& \frac{1-{\rho }^2}{2\rho }\text{Im}\left(H(x,t)e^{-\mathrm{i\Phi }}\right).\end{array}$$

(39)

As three new variables are added to the system, additional constraints have to be defined to guarantee that the transformation (36) determines ρ, Φ, Ψ uniquely. These constraints are defined in [14, 16] as follows:

$$\begin{array}{l}\underset{-\pi }{\overset{\pi }{\int }}\sigma (\psi ,x)e^{i\psi }d\psi =0\\ \text{and}\\ \mathit{\text{Re}}\underset{-\pi }{\overset{\pi }{\int }}\sigma (\psi ,x)e^{i2\psi }d\psi =0.\end{array}$$

(40)

In section “Mathematical model”, we apply this ansatz to populations of oscillators with complex mean field coupling. It is therefore important to determine the order parameters of each subpopulation. In [16] it was shown that

$$Y=R{e}^{\mathrm{i\Theta }}=\int n\left(x\right)Z\left(x\right)\mathit{\text{dx}}=\int n\left(x\right)\gamma \left(x\right)z\left(x\right)\mathit{\text{dx}}$$

holds. This simplifies further if the constants of motion are distributed uniformly in which case γ(x)=1 (see [16]).

In [15, 25] the same general model (35) is treated with a different ansatz. Their idea was to extend the density function W(x,ϕ,t) to a Fourier series:

$$W(x,\varphi ,t)=\frac{n\left(x\right)}{2\pi }\left\{1+\left[\sum _{m=1}^{\infty }{f}_m(x,t)e^{-im\varphi }+\text{c.c.}\right]\right\},$$

where c.c. denotes complex conjugate. Ott and Antonsen have shown that the continuity equation that expresses the conservation of the number of oscillators is fullfilled if the Fourier coefficients f _m can be written in terms of a single function F^m. This set of solutions is the so-called OA-manifold. Pikovsky and Rosenblum argued in [16] that the OA-manifold corresponds to the case of uniformly distributed constants of motion in the Watanabe-Strogatz Theory. Ott and Antonsen discussed in [25] that the OA-manifold is the only attractive region in terms of long-time evolution under the assumption that the parameter distribution n(x) is continuous. This makes it possible to reduce the systems to the OA-manifold for purposes of a long-time analysis.