Network physiology is a new field initiated by Ivanov et al. in 2012, aiming to reveal the relationship between network topology and physiologic function (Bashan et al., 2012; Ivanov and Bartsch, 2014; Bartsch et al., 2015; Ivanov et al., 2016). Generally, a physiological network is consist of different physiologic organ systems, such as the dynamical network of six physiologic systems with nodes being brain activity (EEG waves: δ, θ, α, σ, β), cardiac (HR), respiratory (Resp), chin muscle tone, leg, and eye movements. In contrast to a static complex network, the topology (i.e., number and strength of network links) of physiological network usually vary with time, resulting in different physiologic states such as different sleep stages [deep, light, rapid eye movement (REM) sleep, and quite wake]. That is, each physiological state of physiological network corresponds to a specific network structure. The transition from one physiologic state to another is associated with fast reorganization of physiological interactions. Network physiology can be successfully applied to explain phase synchronization between organ systems (Chen et al., 2006; Xu et al., 2006; Ivanov et al., 2009; Bartsch et al., 2012; Bartsch and Ivanov, 2014). A typical system of network physiology is the brain network, which is the focus of this review.

One of the most challenging and long standing problems is the understanding of the powerful brain functions such as data processing, function approximation, and pattern recognition etc., which has been considered as a black box for a long time. To solve this problem, it is necessary to figure out the brain network first. So far, it is well known that the human cerebral cortex is a huge network, consisting of 10¹¹ neurons and 10¹⁴ links. Thus, it is almost impossible to figure out the structure of this huge network. Fortunately, owing to the modern physical detections such as electroencephalogram (EEG), magnetic resonance imaging (MRI), magnetoencephalogram (MEG), and diffusion tensor imaging (DTI) etc., we can now conveniently obtain a coarse-grained brain network with a finite size of hundreds to thousands nodes. It has been found that the obtained brain network is of a modular structure with a complex connectivity (Hilgetag and Kaiser, 2004), which provides a base for the inherent parallel nature of brain computations (Meunier et al., 2010; Zamora-Lopez et al., 2010). The quantification of complexity in brain networks can be also measured by the multifractality and other factors (Liu et al., 2015; Lin et al., 2016; Xue and Bogdan, 2017; Racz et al., 2018; Gupta et al., 2019; Yang et al., 2019). A larger or more precise brain network can be obtained by reconstructing the connectivity under partial observability assumptions, which is common to many real world settings or experiments (Gupta et al., 2018; Xue and Bogdan, 2019).

By the physical detections, it is revealed that synchronization of neuronal ensembles in the network of cerebral cortex is the base of various neurobiological processes. For example, the alpha synchronization is task-related and is also associated to top-down processing (Benedek et al., 2011). Further, it is pointed out that synchronization and desynchronization of neural activity are closely related to both the normal functions of brain and its disorders, such as epileptic seizures and Parkinson's disease etc. (Rothkegel and Lehnertz, 2014). For instance, during an epileptic seizure, some regions of brain are strongly synchronized but the others are desynchronized (Ayala et al., 1973). While in Parkinson's disease, synchronized activity is absent in the brain regimes of damaged neurons (Levy et al., 2000).

Moreover, it was found that the synchronized and desynchronized behaviors are usually co-existed in the time series of brains. For example, the EEG data showed that during sleep of dolphin, its two brain hemispheres have independently synchronized and desynchronized behaviors at the same time (Mukhametov et al., 1977), i.e., one hemisphere is in sleep and another remains awake. This phenomenon is called as unihemispheric slow-wave sleep and has been also found in other aquatic animals and migrated birds (Rattenborg et al., 2000). On the other hand, the similar phenomenon was reported in the context of neuroscience, called bump state, where the firing rate is higher at some spatial locations but a constant at other spatial positions (Laing and Chow, 2001). This bump state takes an important role for feature selectivity in models of the visual system (Somers et al., 1995), the head direction system (Zhang, 1996), and working memory (Camperi and Wang, 1998).

To understand the mechanisms of these coexisted behaviors in brains, numerous efforts have been paid to the synchronization of coupled oscillators. One of its recent progresses is chimera state, which is closely related to the phenomenon of unihemispheric sleep. Chimera state represents the coexistence of coherent and incoherent dynamics. It is surprising that this behavior occurs in symmetrically coupled identical oscillators. This counterintuitive phenomenon was first discovered in 2002 (Kuramoto and Battogtokh, 2002) and then named in 2004 (Abrams and Strogatz, 2004). Since then, chimera state has become a hot topic and different kinds of chimera states have been revealed in different systems such as the chaotic dynamical systems (Omelchenko et al., 2012), time-delayed system Sethia et al. (2008), and systems with regular topology (Ko and Ermentrout, 2008; Yao et al., 2013; Tian et al., 2017) etc. Initially, chimera states were shown to emerge in systems of nonlocal coupling. Recently, it has been extended to the system of globally coupled oscillators (Chandrasekar et al., 2014), and even in complex networks (Zhu et al., 2014). Further, to explain the alternating activity patterns between the hemispheres over time (Mukhametov et al., 1977), Ma et al. considered the effect of time-delay in two coupled populations and found that the synchronous and desynchronous behavior do alternate between the two groups over time (Ma et al., 2010). At the same time, chimera states have also been implemented in several experiments such as on chemical oscillators (Tinsley et al., 2012), mechanical oscillators (Martens et al., 2013), electronic oscillators (Larger et al., 2013), electrochemical oscillators (Wickramasinghe and Kiss, 2013; Schmidt et al., 2014), and optoelectronic oscillators (Larger et al., 2015). See reviews Panaggio and Abrams (2015) and Majhi et al. (2019) for details.

In sum, chimera states are mainly studied on phase oscillators. As chimera state may represent the mechanism of unihemispheric sleep where the neurons in the sleepy hemisphere are synchronized and the neurons in the awake hemisphere are desynchronized, it is necessary to pay more attention on the chimera state in neural systems. Fortunately, some interesting results have already been obtained in this line, which involve most typical neural models. For examples, chimera states have been studied in leaky integrate-and-fire neurons (Olmi et al., 2010), Morris-Lecar neurons (Calim et al., 2018), FitzHugh-Nagumo neurons (Omelchenko et al., 2013, 2015), Hindmarsh-Rose neurons (Hizanidis et al., 2014, 2016), and Hodgkin-Huxley neurons (Sakaguchi, 2006; Glaze et al., 2016). To accelerate the study of chimera state in brain networks, it is necessary to systematically summarize the measures of chimera state and its recent progress in empirical brain systems, which has not been paid enough attention in the previous reviews (Panaggio and Abrams, 2015; Majhi et al., 2019). Thus, we here briefly summarize the main results of chimera state in neuronal systems and pay a special attention to the network of cerebral cortex.

To characterize the chimera state, three statistical measures have been proposed so far, by using the time series of network. The first measure is the order parameter by (1)ρ(t)eiΦ(t)=1N∑j=1Neiθj(t), where ρ represents the phase coherence of oscillators, Φ(t) is the average phase of all oscillators, θ_j(t) is the phase variable of the j-th oscillator. The system is a complete synchronization when ρ = 1 and a complete desynchronization when ρ = 0. This measure can be conveniently used to the system consisting of two groups. It is a chimera state when one group has an order parameter ρ ≈ 1 and the other ρ ≈ 0.

A key element of this measure is to calculate the phase variable θ_i(t). However, there is not such an explicit variable in all the neuronal models, in contrast to the Kuramoto phase oscillator. Generally, there are two ways to solve this problem. The first one can be used to a general nonlinear oscillator not necessarily having a well-defined rotational center. Let u_i and v_i represent two variables of the i-th neuron, respectively. The phase θ_i(t) can be calculated as (Osipov et al., 2003; Liu et al., 2009) (2)θi(t)=arctan(v.iu.i), where v.i and u.i denote the derivatives of u_i and v_i, respectively. The second one is to calculate the phase θ_i(t) by (3)θi(t)=arctan(viui), provided that u_i and v_i move around the origin. An alternative way of Equation (3) is to calculate its instantaneous angular frequency (Pereira et al., 2007) by (4)θ.i=v.iui-u.iviui2+vi2. Then, the phase θ_i can be integrated from Equation (4).

The second measure is based on the first one but only for the networks having more than two communities (Shanahan, 2010). Consider a network with m communities. We let ρ_i(t)(i = 1, …, m) represent the order parameter for each community i at time t. A chimera state means that the values of ρ_i(t) for different communities are not the same. Based on this feature, two indices are introduced to measure the chimera state, i.e., the index of the metastability λ and the chimera-like index χ. For the former, we first calculate ρ_i(t) for T points, i.e., t ∈ {1 … T}. Then, their variance can be obtained as (5)σmet(i)=1T∑t=1T(ρi(t)-〈ρi(t)〉T)2.

σ_met(i) gives the fluctuation of synchronization in the community-i, i.e., the metastability. The index of the metastability for the entire network is (6)λ=〈σmet(i)〉m. Similarly, an instantaneous variance over all the m communities can be introduced as (Shanahan, 2010) (7)σchi(t)=1m∑i=1m(ρi(t)-〈ρi(t)〉m)2. Then, the chimera-like index χ is (8)χ=〈σchi〉T. Thus, the chimera-like χ and metastability λ indices quantify the degree of synchronization along time and among communities, respectively.

The third measure is for a general network with or without clear communities (Kemeth et al., 2016). Its idea is based on the local curvature for the spatial coherence, represented by the second derivative in the case of one spatial dimension. In this case of one dimension, the local curvature D can be calculated as (9)Df≡f(x+Δx,t)-2f(x,t)+f(x-Δx,t), where f represents the spatial data on a snapshot at time t. Figure 1 shows its schematic figure where (Figure 1A) is a typical snapshot of chimera state and (Figure 1B) is the mapped D, with D_m being the maximal value of |Df|. D_m represents the curvature of the oscillator with its two neighbors being shifted 180^o in phase (Kemeth et al., 2016). From Equation (9) we see that |Df| equals to zero in the synchronous regime and finite values with pronounced fluctuations in the desynchronized regime.

From Figure 1B, we see that |Df| is distributed between 0 and D_m. Letting g be the normalized probability density of |Df|, g(|Df| = 0) is the fraction of spatially synchronous regimes. Thus, g(|Df| = 0) is unity for a complete synchronized state, zero for a complete desynchronized state, and a value between zero and unity for a chimera state. As numerical simulations have fluctuation, it was suggested that those points with |Df| < 0.01D_m should be considered as synchronized, and otherwise desynchronized (Kemeth et al., 2016). That is, the fraction of coherent regions can be calculated by (10)g0(t)=∫0δg(t,|Df|)d|Df|, with δ = 0.01D_m. Therefore, we have g₀ = 1 for a complete coherent, g₀ ≈ 0 for desynchronized, and 0 < g₀ < 1 for chimera states.

Except these three main measures, Gopal et al. introduced another approach to characterize the chimera and multichimera states (Gopal et al., 2014). Their approach is based on a transformation from the framework of x_i to a new framework z_i, i = 1, 2, ⋯ , N, where z_i = x_i − x_i+1. In their approach, chimera states can be measured by two indices, i.e., the strength of incoherence and discontinuity measure. Although these measures are different each other, they are all robust in numerical simulations and can be applied in different systems.

Although chimera state is mainly studied in phase oscillators, it has also been observed in other oscillators with amplitudes, including the neuronal systems. We here make a brief review on the results of chimera state in neuronal systems, mainly focusing on the FitzHugh-Nagumo (FHN) neurons and Hindmarsh-Rose (HR) neurons.

We first introduce the studies of chimera state in FHN neurons. In this aspect, Omelchenko et al. considered the case of FHN neurons by a stronger coupling and found a multi-chimera state (Omelchenko et al., 2013). For the sake of simplicity, they introduced a rotational coupling matrix and found that it is possible for the system to show both chimera and multichimera states. Their discussion is as follows. Consider N FHN neurons coupled nonlocally on a ring (11)ϵdukdt=uk-uk33-vk+σ2R∑j=k-Rk+R[buu(uj-uk)          +buv(vj-vk)],dvkdt=uk+ak+σ2R∑j=k-Rk+R[bvu(uj-uk)          +bvv(vj-vk)], k=1,2,⋯,N where u_k and v_k are the activator and inhibitor variables, respectively. ϵ is taken as ϵ = 0.05. A neuron is excitable for |a| > 1 and oscillatory otherwise (Omelchenko et al., 2013). For simplicity, Omelchenko et al. (2013) considered the case of identical neurons with a_k ≡ a ∈ (−1, 1). In the framework of Equation (11), the coupling strength σ is a constant within the R nearest neighbors from both sides but zero otherwise, marking the feature of nonlocal coupling. To observe chimera states in neuronal systems, Omelchenko et al. introduced the rotational coupling matrix (Omelchenko et al., 2013) (12)B=(buubuvbvubvv)=(cosϕsinϕ-sinϕcosϕ), where ϕ is a parameter of coupling phase, representing the relative phase difference of interacting oscillators. This is a kind of cross-coupling. The value of chosen ϕ determines the property of coupling, i.e., attractive or repulsive.

For convenience, the parameter r = R/N is used to represent the coupling radius. The system of Equation (11) has different behaviors for different r and σ. For example, for fixed a = 0.5 and ϕ = π/2−0.1, a chimera state can be observed for small coupling strength σ and multichimeras state for larger σ. Figure 2 shows the case of small σ where Figure 2A represents a snapshot of u_k at the moment of t = 5, 000. It is a typical chimera state consisting of coherent and incoherent parts. Figure 2B shows a further confirmation where the desynchronized part is distributed along the limit cycle.

On the other hand, we may also use the mean phase velocity of oscillators to characterize the chimera state, which is defined as (13)ωk=2πMkΔT, k=1,2,⋯,N where ΔT is the considered time interval and M_k is the total firing number of node-k in ΔT. Figure 2C represents the results. ω_k is a constant for coherent part but lie on a continuous curve for incoherent part. Moreover, we can also use the local order parameter to represent chimera state, defined as (14)Zk=|12δ∑|j-k|≤δeiθj|, k=1,2,⋯,N where θ_j is defined by Equation (3). The k-th unit is in the coherent part of the chimera state when Z_k = 1 and incoherent parts otherwise. Figure 2D shows the evolution of Z_k by choosing the window size δ = 25, where the yellow parts denote the coherent regions.

For a larger coupling σ, we may observe a multichimera state where the incoherent part is divided into several independent domains, in contrast to the single incoherent domain in chimera state. Figure 3 shows the phase diagram of chimera states on the r − σ plane, where the red region represents the chimera state with one incoherent domain, while the green and blue regions represent the multichimera states with two or three incoherent domains, respectively.

Further, Omelchenko et al. extended chimera state to the case of nonidentical FHN units (Omelchenko et al., 2015). They randomly choose the parameters a_k from a normal (Gaussian) distribution with mean value a_mean and variance δa, so that the neurons have different frequencies. For fixed a_mean = 0.5, chimera state changes with δa. Figure 4 represents a few typical snapshots and corresponding ω_k for chimera states with one, two, and three incoherent regions, as δa increases. It is clear that these chimera states have some robustness to the nonidentity of oscillators (Omelchenko et al., 2015).

Semenova et al. considered the case with noise (Semenova et al., 2016), as it is unavoidable in real systems. They focus on the question whether noise is beneficial for chimera states. In fact, this consideration is related to the phenomenon of coherence resonance (Hu et al., 1993; Pikovsky and Kurths, 1997; Liu and Lai, 2001), where an optimal noise intensity can result in the counterintuitive increase of temporal coherence. Semenova et al. found that noise is essential for chimera behavior and call it coherence-resonance chimera (Semenova et al., 2016). That is, an optimal noise can even induce a spatial chimera state, rather than purely temporal coherence. Their model is a slight modification of Equation (11) and can be written as (15)ϵduidt=ui-ui33-vi+σ2R∑j=i-Ri+R[buu(uj-ui)          +buv(vj-vi)],dvidt=ui+ai+σ2R∑j=i-Ri+R[bvu(uj-ui)          +bvv(vj-vi)]+2Dξi(t), where ξ_i(t) represents the Gaussian white noise with 〈ξ_i(t)〉 = 0 and 〈ξi(t)ξj(t′)〉=δijδ(t-t′), and D denotes the noise intensity. The rotational coupling matrix and phase parameter are fixed as the same as in Equation (11) except a_i ≡ a = 1.001, indicating that all neurons are excitable but close to the threshold.

The presence of noise may induce a spike for a single FHN neuron in the excitable regime of |a| > 1. When the noise intensity is too small, it is not enough to induce a spike in the system of Equation (15). While a strong enough noise may induce too many spikes and results in irregular behaviors. Thus, there may exist an optimum intermediate noise to induce a chimera state. To confirm this analysis, Figure 5 shows four distinct regimes by u_i and Z_i from Equation (14). Figure 5A shows the case of D = 0, with no spikes. Figure 5B shows the case of intermediate noise with D = 0.0002. The values of Z_i show the characteristic feature of chimera state, i.e., the coexistence of coherent and incoherent spiking. Thus, Semenova et al. called it as coherence-resonance chimera (Semenova et al., 2016). Figure 5C shows the case of strong noise with D = 0.0004, where the coherence-resonance chimera is destroyed. Figure 5D shows the case of much stronger noise with D = 0.1, where the spiking is incoherent in both time and space.

To see how the coupling parameters r and σ influence the coherence resonance chimera, Figure 6 shows the phase diagram by fixing the parameters ϵ, a, D, N. Region (a) represents a homogeneous steady state, region (b) shows the state of spiking patterns with temporal periodicity and spatial incoherence, regions (c–e) represent the coherence-resonance chimeras with one, two, and three incoherent domains, respectively. Therefore, except the noise intensity D, the coherence-resonance chimeras are also influenced by the two parameters r and σ.

Further, the values of Z_i in Figure 5B show a periodic switching between the coherent and incoherent regimes of chimera state. This feature may be helpful for the understanding of unihemispheric sleep, where the coherent and incoherent behaviors are known to switch between the two hemispheres of brain (Mukhametov et al., 1977). In fact, this kind of alternating chimera behavior has been previously addressed in a phase model by a sinusoidal signal (Ma et al., 2010), which can be also considered as an external perturbation.

Tian et al. considered the case of time-delay and electromagnetic induction in FHN neurons and found that either the time delay or electromagnetic induction can induce chimera states (Tian et al., 2018). By considering both the effect of time-delay and electromagnetic induction, Equation (11) becomes (16)ϵu.i=ui-ui33-vi+σ2R∑j=i-Rj=i+R[uj(t-τ)-ui]          +kρ(φi)ui,v.i=ui+a,φ.i=k1ui-k2φi, where φ represents the magnetic flux and τ denotes the time-delay. The term kρ(φ_i)u_i is the induction current (Tian et al., 2018). The nonlinear function ρ(φ) is taken as ρ(φ) = α + 3βφ² (Ma et al., 2017), where α and β are two parameters.

In numerical simulations, the parameters are fixed as N = 256, k₁ = 0.1, k₂ = 1.0, α = 0.1, β = 0.1, τ = 1.0, σ = 0.02 and r = 0.35. Figure 7 shows the results where Figures 7A–D denote the cases of k = 0.8, 1.2, 5.2, and 7.4, respectively. It is clear from Figures 7A–D that with the increase of k, the multi-chimera state gradually becomes a chimera state.

Now, we turn to the studies of chimera state in HR neurons. Different from the FHN model of spiking neurons, HR model may represent the bursting behavior of neurons. We here concern about how this bursting behavior influences the chimera state.

Bera et al. studied a network of bursting HR neurons with global coupling as follows (Bera et al., 2016) (17)ẋi=axi2-xi3-yi-zi+kN-1(vs-xi)∑j=1NcijΓ(xj),ẏi=(a+α)xi2-yi,żi=c(bxi-zi+e), i=1,2,…,N where k is the coupling strength, x is the membrane potential, and y and z are the transport of ions across the membrane through the fast and slow channels, respectively. The HR neuron is excitatory for x_i(t) < v_s, where v_s = 2 is the reversal potential. (c_ij) is the adjacent matrix with c_ij = 1 if i ≠ j and c_ii = 0. The coupling function Γ(x) is assumed to be the sigmoidal nonlinear function (18)Γ(x)=11+e-λ′(x-Θs), where λ′ determines the slope of the function and Θ_s is the synaptic threshold. These two parameters are taken as Θ_s = −0.25 and λ′ = 10.

To characterize the chimera state, a new transformed variable w_1,i = x_i − x_i+1 is introduced (Bera et al., 2016). Figure 8 shows the results for snapshots of the state variables x_i and w_1,i by black and red color dotted points, respectively. Figure 8A shows the case of a weak coupling k = 1.0. It is a disordered state. Figure 8B shows the case of a middle coupling k = 1.2. It is a multichimera state with two domains of disordered oscillators. The inset shows a typical time series of x_i (blue color line). Its behavior changes between the square-wave and plateau bursting. Figure 8C shows the case of a strong coupling k = 1.28. Its chimera state has only one incoherent domain, in contrast to the two incoherent domains in Figure 8B. Figure 8D shows the case of a stronger coupling k = 1.3, where all the neurons become coherent.

Hizanidis et al. considered a modular network of HR neurons (Hizanidis et al., 2016), which is from the C. elegans network and consists of six communities. They let the neurons be connected by two types of synapses: electrical and chemical. The former is for the connections within each community and the latter for the connections across the communities. Their model can be described as follows

where p_i represents the membrane potential, q_i and n_i are associated with the fast and slow currents, respectively. The parameters are chosen as r = 0.005, a = 1, V_syn = 2, b = 3, d = 5, c = 1, s = 4, p₀ = −1.6 and I_ext = 3.25 so that each neuron has the spiking-bursting behavior (Hizanidis et al., 2016). L is the Laplacian matrix with L_ij = E_ij − δ_ijk_i, where δ_ij = 1 if i = j, and δ_ij = 0 otherwise. E is the adjacency matrix with E_ij = 1 if there is an electrical synapse connecting the neurons i and j, and E_ij = 0 otherwise. g_el is the strength of the electrical coupling and its functionality is governed by the linear function H(p) = p. The adjacency matrix T is T_ij = 1 if there is a chemical synapse between neurons i and j, and T_ij = 0 otherwise. g_ch is the strength of the chemical coupling and its functionality is defined by the sigmoidal function Γ(p) from Equation (18).

To study the chimera state of Equation (19), Hizanidis et al. (2016) used ρ to represent the order parameter of Equation (1) where θ is calculated through Equation (4) by p and q. Figures 9A–G show the phase diagrams of ρ on the (g_ch, g_el) parameter space for each of the six communities and for the entire network, respectively. It is clear that ρ is not homogeneously distributed in the (g_ch, g_el) plane but with higher ρ in some region and lower ρ in other regions. For example, the red regions have ρ ≈ 1 and the yellow regions have 0 < ρ < 1.

The regions with 0 < ρ < 1 in Figures 9A–G may represent the chimera states. To confirm it, the two measures of metastability index λ of Equation (6) and chimera-like index χ of Equation (8) are used. Figures 9H,I show the results on the (g_ch, g_el) parameter space. From Figure 9H we see that λ reaches higher values in some regions, implying that the system often changes between coherent and incoherent states. From Figure 9I we see that χ reaches its highest values in the two synchronization “islands” of communities 3 and 6. For more detailed information on chimera state, Hizanidis et al. chose 3 interest points on the (g_ch, g_el) phase diagram, marked by letters A, B, and C (Hizanidis et al., 2016), where A has both low-valued λ and χ, B has λ ≫ χ, and C has χ ≫ λ, i.e., “chimera-like” state. Figure 10 shows the dynamical behaviors of p for the three points. The snapshots of the system state in the bottom confirm the corresponding behaviors, where only the point C shows the feature of chimera state, i.e., the coexistence of coherent and incoherent domains.

Moreover, Majhi et al. studied the chimera state in uncoupled HR neurons induced by a multilayer network (Majhi et al., 2016), where the neurons in the upper layer is unconnected but can share information through the neurons in the lower layer. This topology is related to the remote synchronization (Bergner et al., 2012) and thus may help us to understand brain functions. Bera et al. reported a new type of non-stationary chimera pattern in coupled HR neurons (Bera et al., 2019), called spike chimera.

In summary, these studies show different ways to induce chimera states, including both cross-coupling and single-variable coupling. At the same time, these studies also show that chimera states are available in different neuronal models, indicating the robustness to neuronal models. The drawbacks or limitations are that all the considered network structures are artificial but not empirical brain networks. And the models are for individual neurons but not the average behaviors of an ensemble of neurons, which are the only available time series in experiments such as EEG data.

Except the above extensive studies of chimera states in artificial neural systems, recently, some attention has been paid to the networks of human cerebral cortex measured by DTI. The former helps us to understand the mechanism of chimera state such as how chimera state is induced and what is the condition for chimera state to show up. While the latter highlights a way to cure or control brain diseases such as schizophrenia, Alzheimer's disease and brain tumors. It is well-known that in the fields of nonlinear science and complex network, brain functions can be represented by their corresponding dynamical patterns, i.e., a variety of patterns of partial synchronization. These patterns have a close relationship with chimera states and can be considered as a natural link between coherent and incoherent dynamics. Thus, the studies of chimera states on empirical brain networks is very helpful for exploring the mechanism of brain functions such as cognition and memory. We here make a brief summary for those results on the empirical brain networks, i.e., from the brain networks with smaller size to the middle and then to the larger ones.

Firstly, we introduce the study on an empirical brain network with smaller size. In this case, Bansal et al. considered an empirical brain network consisting of 76 brain regions (or nodes) and paid attention to how brain structure influences the dynamical patterns produced by stimulation (Bansal et al., 2019). They divided this network into nine cognitive systems by using personalized brain network models, named attention, visual, cingulo-opercular, subcortical, medial default mode, somatosensory and motor, frontoparietal, ventral temporal association, and auditory systems. Each of the nine cognitive systems is consist of the coactivated regions for supporting a generalized class of cognitive functions. Then, they presented a chimera-based, cognitively informed framework to study how large-scale brain structure influences brain dynamics and functions, called cognitive chimera states. In their study, the dynamics of each node was modeled by the Wilson-Cowan oscillators (Wilson and Cowan, 1972), a biologically motivated neural mass model, represented as follows (20)τdEidt=-Ei(t)+(SEm-Ei(t))SE(c1Ei(t)-      c2Ii(t)+c5∑jAijEj(t-τdij)+Pi(t)),τdIidt=-Ii(t)+(SIm-Ii(t))SI(c3Ei(t)-      c4Ii(t)+c6∑jAijIj(t-τdij)), where (21)SE,I(x)=11+e-aE,I(x-θE,I)-11+eaE,IθE,I and A_ij is the weighted coupling matrix. c₅ and c₆ represent the excitatory and inhibitory coupling strength, respectively, with c₆ = c₅/4. P_i(t) is the external stimulation. τdij=dij/td is the time-delay, where d_ij is the spatial distance between nodes i and j and t_d = 10m/s is the signal transmission speed. Other parameters are biologically taken as c₁ = 16, c₂ = 12, c₃ = 15, c₄ = 3, τ = 8, θ_I = 3.7, θ_E = 4, a_I = 2, and a_E = 1.3.

In numerical simulations, the dynamical behaviors of the nine cognitive systems were investigated by stimulating a brain region with P_i = 1.15 (Bansal et al., 2019). The stimulation gradually spreads to other parts of the brain network by the links of the stimulated node and form a dynamical state. By this way, different dynamical states are observed when different brain regions are stimulated. Then, Bansal et al. calculated a cognitive system-based order parameter ρ_si,sj from Equation (1). This parameter ρ_si,sj measures the degree of synchrony among all oscillators within the two cognitive systems s_i and s_j. By this way, a cognitive system-based 9 × 9 matrix can be obtained. Further, two cognitive systems s_i and s_j are considered as synchronized if ρ_si,sj exceeds a threshold value ρ_Th = 0.8. The whole brain network is a coherent state when all nine cognitive systems are synchronized, a cognitive chimera state when some cognitive systems are synchronized while the other systems are desynchronized, and a metastable state when no stable synchrony between cognitive systems is formed (Bansal et al., 2019).

Figure 11 shows the relative contribution of the nine cognitive systems for the three states after nodes stimulation. We see that the contributions to the three states are different from one cognitive system to another, implying a close relationship between cognitive systems and brain functions. Figure 11A shows that coherent states are mainly from the nodes stimulation within subcortical and medial default mode systems. Figure 11C shows that the frontoparietal and cingulo-opercular systems are the main contributions for metastable states, while the ventral temporal association and auditory systems also contribute substantially to metastable states. Figure 11B shows that all nine systems produce chimera states, implying that chimera states have higher possibility to be observed than either coherent or metastable states (Bansal et al., 2019). Bansal et al. also pointed that the regions of coherent states are distributed more closely to the midline of the brain, and the regions of metastable states are distributed farther from the midline, i.e., along the edges of the hemispheres, while the regions of chimera states are relatively uniformly distributed within the brain space. The metastable states enable segregated neural processing, while coherent states enable integrated neural processing. As the brain system must integrate information across spatially distributed, segregated regions to implement cognitive tasks, a balance between integration and segregation is required for adaptive cognition. This balance is automatically satisfied in chimera states and thus enables segregation and integration in brain dynamics, which guarantees the diverse processing requirements.

Secondly, Chouzouris et al. considered a slightly larger empirical brain network consisting of 90 nodes (Chouzouris et al., 2018), motivated by studies of epileptic seizures. They let the neurons be the FitzHugh-Nagumo oscillators. As this empirical brain network has a topology of complex network, Equation (11) has to be modified into (22)ϵdukdt=uk-uk33-vk+σ∑j=1NGkj[buu(uj-uk)      +buv(vj-vk)],dvkdt=uk+a+σ∑j=1NGkj[bvu(uj-uk)      +bvv(vj-vk)], where ϵ = 0.05 and G is the adjacency matrix. The rotational coupling matrix B = (b_uu) is taken as the same as in Equation (12).

Chouzouris et al. found that when a is in the range a ∈ (0, 0.8), the network exhibits chimera states for small coupling strength σ (Chouzouris et al., 2018). When the parameters are taken as a = 0.5, σ = 0.2, and N = 90, the network shows a stationary moving chimera. That is, the order parameter r from Equation (1) has a strong fluctuation or changes in time. Figure 12A shows such an example, where the parameters are a = 0.5 and σ = 0.6. The average of order parameter is 〈r〉 ≈ 0.5. Except the strong fluctuation, another feature is that the highest value of r appears right before its drop. Both effects were discovered in the synchronization of epileptic seizures (Jiruska et al., 2013). Chouzouris et al. further pointed out that the high coherence events can be controlled (Chouzouris et al., 2018). Larger σ increases the probability for chimera states to occur. Figure 12B shows the result that the variation of coupling results in a switching between the chimera state and synchronization, which controls the epileptic seizures.

This empirical brain network of 90 nodes was further studied by Ramlow et al. (2019), where a dynamical asymmetry between the hemispheres was addressed by considering natural structural asymmetry.

Finally, we discuss the situation of empirical brain network with larger size. This network is much larger than the above two and is constructed by the data of Hagmann et al. (2008) where the cerebral cortex was divided into 989 nodes and 17, 865 connections. For this network, Huo et al. first considered the case of adaptive coupling (Huo et al., 2019), based on the fact that in empirical brain network, both the coupling strength and neural activities influence each other and thus change with time. In Huo's model, the coupling matrix is adaptively evolved with the dynamics of neurons. They found that the adaptive coupling finally reaches a self-organized state and induces chimera states. This kind of self-organization may support the high flexibility of brain functions. In details, they let the nodes be the FitzHugh-Nagumo oscillators and the dynamics be represented as follows (23)ϵu.i=ui-ui33-vi(t)      +1N∑j=1Nλij[buu(uj-ui)+buv(vj-vi)],v.i=ui+a+1N∑j=1Nλij[bvu(uj-ui)+bvv(vj-vi)], where ϵ = 0.05, a = 0.5, and the rotational coupling matrix B is defined as in Equation (12). The adaptive coupling λ_ij is set as (24)λ.ij=-γ[sin(uj-ui+β)+λij], where γ is a small constant. Equation (24) does not influence the topology of network but only change the value of λ_ij.

Huo et al. chose γ = 0.01 and found that it is possible for chimera states to appear in this realistic network (Huo et al., 2019). Figure 13 represents the resulted chimera state for β = −0.5π + 1.2 and ϕ = −π + 4.45, where (Figure 13A) is the initial matrix of λ_ij, (Figure 13B) the stabilized matrix of λ_ij, (Figure 13C) the evolutionary pattern of u_i(t), and (Figure 13D) a snapshot of the fast variable u_i(t). Comparing Figure 13A with (Figure 13B), we see that the stabilized matrix λ_ij in (Figure 13B) is significantly different from the initial matrix λ_ij in (Figure 13A). This may help us to understand how the brain network is self-organized into chimera states.

To reflect the dependence of chimera states on the parameters ϕ and β, Huo et al. calculated the measure g₀ of chimera state by Equation (10). Figure 14 shows its phase diagram on the parameter plane of ϕ and β. The stabilized behaviors consist of disorder, coherent, and chimera states but the fraction of chimera states is the largest one on the phase diagram.

Kang et al. further considered the empirical brain network of 989 nodes as a two-layered network where the left and right hemispheres of cerebral cortex are considered as different layers, respectively (Kang et al., 2019). In their model, the intra- and inter-coupling strengths are considered to be different. Very interesting, they found that the model can reproduce the phenomenon of unihemispheric sleep where one hemisphere is completely synchronized while the other is completely desynchronized. This finding provides an explanation for the first-night effect in human sleep (Tamaki et al., 2016). Their model for the first layer-A is represented as follows (25)ϵu.ia=uia-(uia)33-via+λinkin,ia∑j=1NAij[duu(uja-uia)+duv(vja-via)]+λoutkout,ia∑j=1N(AB)ij[duu(ujb-uia)+duv(vjb-via)],v.ia=uia+a+λinkin,ia∑j=1NAij[dvu(uja-uia)+dvv(vja-via)]+λoutkout,ia∑j=1N(AB)ij[dvu(ujb-uia)+dvv(vjb-via)], where ϵ = 0.05 and a = 0.5. kout,ia and kin,ia are the inter- and intra-degrees of node i, respectively. (AB)_ij and A_ij denote the inter- and intra-coupling matrices, respectively. The other quantities are the same as in Equations (11) and (12).

The similar dynamics equations can be written for the network B.

Figure 15 shows the results for four typical cases where ω_i is calculated by Equation (13), the up and down panels are for the network-A and network-B, respectively, and the insets are their corresponding dynamics of u_i. The panels of Figure 15 represent four typical cases where (Figures 15A,E) are for the case of disorder with λ_in = 0.1 and λ_out = 0.3; (Figures 15B,F) for the case of chimera state with λ_in = 0.1 and λ_out = 1.8; (Figures 15C,G) for the case of disordered network-A and synchronized network-B, with λ_in = 0.4 and λ_out = 3.5; and (Figures 15D,H) for the case of synchronization with λ_in = 4.0 and λ_out = 3.5. These four cases represent different states. The first case of Figures 15A,E and the last case of Figures 15D,H denote the two extreme states of desynchronized and synchronized states, respectively. The second case of Figures 15B,F represents a chimera state where there is a plateau of ω_i in both the up and down panels and their insets show a coexistence of synchronized and unsynchronized u_i(t). The most interesting is the third case of Figures 15C,G where the network-A is disordered but the network-B is synchronized, marking the unihemispheric sleep. Kang et al. further showed that the parameter region for the state of unihemispheric sleep is much smaller than that of chimera state, implying that it is usually difficult to observe the phenomenon of unihemispheric sleep. This is consistent with the first-night effect (Tamaki et al., 2016), which can be observed only in the first-night sleep when a person is located in an unfamiliar place.

In sum, these studies showed the recent progress of chimera states in empirical brain networks, but did not pay much attention to the aspects of characteristic features of brain networks, such as the heterogeneous communities and hub nodes of rich-club, and deeper connection to concrete brain functions, such as cognitive and memory etc.

Chimera state is in fact one of the three kinds of partial synchronization. The other two of them are the cluster synchronization and remote synchronization. Cluster synchronization represents the case where the oscillators of network are automatically evolved into different synchronized clusters but the oscillators in different clusters are not synchronized each other (Schaub et al., 2016; Cao et al., 2018). The relationship between the synchronized cluster and network symmetry is discussed recently (Pecora et al., 2014; Sorrentino et al., 2016), including the case where the synchronized cluster is not directly from the symmetry but due to the same total amounts of inputs received from their neighboring nodes (Siddique et al., 2018). While remote synchronization represents the synchrony among the leaf nodes of a hub but not synchronized with the hub, i.e., the synchronized nodes are not directly connected (Bergner et al., 2012). However, these three partial synchronization are not completely independent of each other but may sometimes represent the same phenomenon. For example, chimera state can appear simultaneously with cluster synchronization in some systems (Hart et al., 2016; Cho et al., 2017; Bansal et al., 2019). It is also possible for remote synchronization to be related to cluster synchronization (Bergner et al., 2012; Kang et al., 2020; Wang and Liu, 2020). For example, in a star network with remote synchronization, if the leaf nodes of the hub is considered as a cluster, their synchronization is in fact the cluster synchronization. Based on these results, it is an open but promising direction to highlight the mechanisms of brain functions such as cognition, memory, and signal spreading etc, from these three aspects of partial synchronization.

One purpose of studying chimera states in neural systems is for its potential applications. On one hand, some attention has been paid to the phenomenon of unihemispheric sleep (Rattenborg et al., 2000; Ma et al., 2010). On the other hand, the study of chimera state may help us to understand neuronal diseases such as epileptic seizures, Parkinson's disease, schizophrenia, Alzheimer's disease and brain tumors (Uhlhaas and Singer, 2006). For example, a therapy for Parkinson's disease is external electric stimulation at high frequencies, called deep brain stimulation (Benabid et al., 1991). These studies are still very primary. More deeper studies are expected, such as the connections of chimera states to the mechanisms of cognitive and memory and the control of various brain diseases in clinical medicine etc.

ZL and ZW conceived the research project and wrote the paper. All authors contributed to the article and approved the submitted version.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.