Chimera states in uncoupled neurons induced by a multilayer structure

Spatial coexistence of coherent and incoherent dynamics in network of coupled oscillators is called a chimera state. We study such chimera states in a network of neurons without any direct interactions but connected through another medium of neurons, forming a multilayer structure. The upper layer is thus made up of uncoupled neurons and the lower layer plays the role of a medium through which the neurons in the upper layer share information among each other. Hindmarsh-Rose neurons with square wave bursting dynamics are considered as nodes in both layers. In addition, we also discuss the existence of chimera states in presence of inter layer heterogeneity. The neurons in the bottom layer are globally connected through electrical synapses, while across the two layers chemical synapses are formed. According to our research, the competing effects of these two types of synapses can lead to chimera states in the upper layer of uncoupled neurons. Remarkably, we find a density-dependent threshold for the emergence of chimera states in uncoupled neurons, similar to the quorum sensing transition to a synchronized state. Finally, we examine the impact of both homogeneous and heterogeneous inter-layer information transmission delays on the observed chimera states over a wide parameter space.

Scientific RepoRts | 6:39033 | DOI: 10.1038/srep39033 subject to a common noisy field 34 and also under neuron's membrane potential stimulation 35 . Experimentally it was observed that synchronization in different neurons may appear in the same region of the brain and even in different regions in the brain 36 . On the other hand, most of the earlier works on chimera states assume that the oscillators are on a single layer whatever be the network topology is. Again, there are many physical systems that do not interact directly but exchange the information through a common medium, for instance, in the Huygen's experiment the two pendulum clocks were interacting through the common wooden beam from which they were hanging. Also this type of indirect interaction is particularly important in biological systems, e.g., populations of cells in which oscillatory reactions are taking place 37 through the interaction via chemicals that diffuse in the surrounding medium.
In this paper, we investigate an architecture where the neuronal oscillators in the upper layer (the layer of our interest) have no connection among them, while they interact with each other via another layer of similar oscillators, thus forming a multilayer structure. To the best of our knowledge, this is the first observation of chimera states in uncoupled neurons in the form of multilayer network. Recent research and reviews attest to the fact that multilayer networks are the next frontier in network science [38][39][40][41][42][43][44][45][46][47][48][49][50][51][52] . In particular, not only are the interactions between the constituents of a complex system limited and thus best described by networks, it is also a fact that the processes happening in one network may vitally affect another network, and moreover, that a node in one network is likely part of another network. From the world economy and transportation systems to social media and biological systems, it is clear that such interdependencies exist, thus making networks of networks or multilayer networks a much better and realistic description of such systems.
This type of multilayer structure is also quite evident in neuronal networks. So it would be of obvious importance to study a network having an architecture of the this type taking neuronal models as the nodes of both the layers. We take the neuron dynamics in terms of Hindmarsh-Rose neuronal model which exhibits various types of bursting dynamics such as spiking, plateau bursting, square-wave bursting (periodic and chaotic) and mixed mode bursting etc. depending on the system parameters. Two types of synapses i.e. electrical and chemical synapses exist through which neurons communicate with each others. Moreover, neurons may not be connected with each other with the same type of synapses everywhere. In fact, a recent work on chimera-like states has been done on neural network inspired by the neuronal connection of the C. elegans soil worm, organized into six interconnected communities, assuming that the neurons are connected with electrical synapses within the communities and with chemical synapses across them 12 . In this study, we consider that the neurons are connected with electrical synapses within the bottom layer (the medium) and with chemical synapses across the layers. Again in neuronal networks, delays arise due to finite propagation times along the axons or to reaction times at chemical synapses. Also depending on the physiological properties of axons and synapses, these delays in signal transmission between different cells of the network may differ. Thus time delay in inter-layer information transmission process is indisputable and it can be heterogeneous too. Hence our aim in this article is to examine how these two types of synaptic connection affect the dynamics of upper layer and how chimera pattern emerges due to that competing effects of two synapses and finally to study the influence of inter-layer synaptic delay (both homogeneous and heterogeneous) on the upper layer dynamics.
We consider N identical isolated neurons which are connected through a common medium. We assume that isolated neurons are situated on same layer (upper layer) and medium as globally connected neurons (lower layer). Each isolated neuron interact directly with one neuron (its replica) in the medium. As we are considering global coupling in the multi-layering layer (common medium) and the uncoupled neurons are only interacting with its replica in the common medium, so the spatial order of the neurons in the upper layer is same as the spatial order of the neurons in the common medium. The schematic diagram of the network is shown in Fig. 1. Therefore, in terms of recently developed multi-layer networks, our proposed network is a multilayer network with two layers. We consider each neuron as Hindmarsh-Rose neuron model. The neurons in the medium are connected through electrical synapses as they allow direct and passive flow of electron via gap junctions. These gap junctions permit for mutual instantaneous transmission of electron between the neurons which are spatially very close to each other. The main goal for considering electrical synapses in the lower layer (medium) is to synchronize (or coherent motion) electrical activity among the neurons. We assume electrical synapses among the globally coupled neurons in the common medium which is homogeneous distribution of the medium. This assumption is more realistic in biological and chemical systems, as homogeneous medium is observed either in biological systems by the fast diffusion of the small molecules or in chemical systems by stirring the solution. On the other hand, the isolated neurons are connected with the common medium through chemical synapses. The chemical synapses typically function in a longer range compared to electrical synapses, it would thus be more likely to connect across the layers. Thus the simultaneous effect of electrical and chemical synaptic coupling is best represented by multilayer structure where the neurons in the medium are connected through electrical synapses while across the layers through chemical synapses.

Results
In this work, we study the emergence of symmetry breaking pattern in the upper layer as a result of the co-action of the two types of synapses. Here we are considering multilayer network and the number of neurons in both the layers are same. We investigate two different cases based on inter-layer coupling delays. In first case, we consider the instantaneous inter-layer chemical synaptic coupling between the layers and later the effect of delays (homogeneous and heterogeneous) present in the inter-layer chemical synaptic coupling.
Instantaneous inter-layer chemical synaptic interaction. In this section, we mainly investigate the dynamics of the isolated neurons in absence of inter-layer chemical synaptic coupling delay. For K ch = 0 (in Eq. (4) of Method section), the bottom layer gets synchronized quite rapidly due to the global (all-to-all type) interaction between the neurons with electrical synapses (particularly, K el = 1.0) and the neurons in the upper-layer behave Scientific RepoRts | 6:39033 | DOI: 10.1038/srep39033 according to their individual rhythms (i.e., the rhythm of square wave bursting). So, our objective is to explore the dynamics of the uncoupled neurons in the upper layer while activating the inter-layer chemical synaptic coupling strength K ch and keeping the bottom layer neurons synchronized. In this case the uncoupled neurons are connected with the common medium. Now switching on K ch , we initially see the incoherent (desynchronized) behaviors of the upper layer neurons and they remain incoherent for 0 ≤ K ch < 1.075. But as we increase K ch , the upper layer network spontaneously splits into two coexisting domains, one of which is coherent and the other one is incoherent which portrays the feature of chimera states. If we increase K ch further, we observe that all the neurons get synchronized for K ch > 1.230. Figure 2(a,b,c) show the snapshots of membrane potentials of all the uncoupled neurons in the upper layer exhibiting incoherent, chimera and coherent states at K ch = 1.0, 1.13 and 1.30 respectively. At these points, the coherent behavior of all the neurons in the common medium are illustrated in insets of Fig. 2(a,b,c). Middle panel of Fig. 2 shows the behaviors of particular three neurons x 20,1 , x 50,1 and x 90,1 illustrated by green, blue and red colors respectively. At incoherent states, all the upper layer neurons follow chaotic square-wave bursting dynamics ( Fig. 2(d)). At higher value of inter-layer chemical synaptic coupling K ch = 1.13, a typical pattern of chimera state with one group of coherent neurons and another group of incoherent neurons coexist. In this case, a neuron in the coherent group and a neuron in the incoherent group have the same time series form i.e. chaotic square-wave bursting in nature shown in Fig. 2(e). At higher value of k ch = 1.3, all the uncoupled neurons in the upper layer are found to be in coherent state and the neurons are in plateau bursting states ( Fig. 2(f)). Additionally, we calculate mean angular frequency 53 of the i-th neuron as, is the geometric phase for the fast variables x i,1 and y i,1 of the i-th neuron, which is a good approximation as long as c is small (≪ 1) and 〈 … 〉 t denotes long term time average. The mean angular frequencies corresponding to incoherent, chimera and coherent states are shown in Fig. 2(g,h,i) respectively. To calculate mean angular frequencies ω i,1 , the time interval is taken over 5 × 10 5 time units after an initial transient of 3 × 10 5 units, throughout the paper. The mean angular frequency corresponding to the neurons in incoherent group are randomly scattered whereas for coherent group of neurons they are same. These mean angular frequency profiles clearly distinguished coherent and incoherent groups in the chimera state. Figure 3 depicts the variation of strength of incoherence (SI) (refer to the Method section) with respect to inter-layer chemical synaptic coupling strength K ch . As can be seen, in the region I = {K ch :1.0 ≤ K ch < 1.075}, the value of SI remains unity characterizing the incoherent (disordered) neurons but as we increase K ch beyond K ch = 1.075, we observe chimera state characterized by the values 0 < SI < 1 in the region II = {K ch :1.075 ≤ K ch ≤ 1.230}. Although 0 < SI < 1 may represent other dynamical states like cluster state 54 , splay state etc., rigorous verification of the snapshots and time series of the neurons in the parameter range II (also for all the parameter region plots in the following sections) have confirmed the existence of chimera patterns only. With further increase, the value of K ch leads to the coherent state as the values of SI becomes zero in the region III = {K ch :K ch > 1.230}.
Quorum sensing mechanism for chimera states. Next, we find the density-dependent threshold for the emergence of chimera states in the upper-layer. This density-dependent threshold is a similar entity like that in quorum-sensing transition to synchronization 55,56 . This mechanism plays a key role in bacterial infection, biofilm formation and bioluminescence 57 . In the context of neuronal network, a similar quorum sensing mechanism involve local field potential 58,59 which may exist through a different level in cortical hierarchy and play an important role in the synchronization of group of neurons. This mechanism also exists in the synchronization of   Fig. 2(a,b and c) respectively. chemical oscillators 60 and cold atoms 61 . In fact, many natural synchronization phenomena where the individual oscillators are not directly coupled, but coupled rather through a common medium experience different synchronization regimes as a function of the number of uncoupled nodes or their density. In our work, by increasing the number of uncoupled neurons in the upper layer (as well as the number of neurons in the lower layer) which are interacting through lower layer (medium), an emergence of chimera states is observed in the upper layer. For small number of uncoupled neurons, say N = 22, we observe chimera states for very small range of K ch (1.12 ≤ K ch ≤ 1.13) and at higher value of it gives coherent state. Chimera states are not identified for N < 22 and in this case all the neurons are either in disordered or coherent state depending on the value of K ch . The snapshot of the membrane potential x i,1 for N = 22 and K ch = 1.125 is shown in Fig. 4(a). If we increase the number of uncoupled neurons in the upper layer, the range of K ch increases for the existence of chimera states. For N = 40, chimera emerges for the inter-layer chemical synaptic coupling strength K ch in (1.13 ≤ K ch ≤ 1.175) and snapshot of x i,1 at K ch = 1.15 are shown in Fig. 4(b). Figure 4(c) shows the chimera states for N = 60 and K ch = 1.16. Figure 4(d,e,f) show the mean angular frequencies ω i,1 corresponding to the chimera states in Fig. 4(a,b,c) respectively. To calculate mean angular frequencies ω i,1 , the time interval is taken over 5 × 10 5 time units after an initial transient of 3 × 10 5 . It is observed that the chimera states persist for long time range in the case of small number of uncoupled neurons as well. To explore the complete dynamics of the uncoupled neurons, we plot the phase diagram in the N − K ch parameter space (Fig. 4(g)) for the range of N ∈ [10, 100] and K ch ∈ [1.0, 1.5]. From this figure, it is seen that the region of chimera states in the inter-layer synaptic coupling strength K ch increases with increasing density N of the neurons.
Homogeneous inter-layer chemical synaptic coupling delay. As mentioned earlier that the inter-layer information transmission may not be instantaneous, in general, and so the time-delay in this process should not be neglected. Previous notable works include enhancement of synchrony in a network of Hindmarsh-Rose neuronal oscillators with time-delayed coupling 62 . Influence of time delay in the context of control of synchronization is studied in coupled excitable neurons in ref. 63. Impact of information transmission delay on the synchronization transitions of modular networks in presence of both electrical and chemical synapses has also been studied in ref. 64. Here our aim is to analyze how the dynamical state in the upper layer varies in presence of time-delay in the inter-layer chemical synapses. At first we take the homogeneous inter-layer delay, i.e. the delay in the information transmission from upper to lower layer and from lower to upper layer are same i.e. τ 1 = τ 2 = τ. Keeping τ fixed at a certain value and varying K ch we observe chimera patterns again as a link between incoherence and coherence as before. In Fig. 5, we fix the inter-layer chemical synaptic coupling delay as τ = 0.4 and vary the inter-layer chemical synaptic coupling strength K ch . Here again for small values of K ch , the upper layer remains disordered but as it increases, chimera pattern arises and sustains for a longer range of K ch compared to the previous case of instantaneous inter-layer coupling. Further increase in K ch gives rise to coherent dynamics in the layer. In the left panels of Fig. 5, the snapshot of amplitudes (membrane potential) depicting disordered, chimera and coherent states for K ch = 0.43, 0.73 and 1.10 are shown in Fig. 5(a-c) respectively. Corresponding mean angular frequencies for disordered, chimera and coherent states are given in Fig. 5(d-f).
As a characterization of the chimera states, we plotted the strength of incoherence SI, for different values of K ch . As in Fig. 6, the upper layer remains disordered in the region I = {K ch :0.4 ≤ K ch < 0.57} with SI = 1, chimera pattern emerges in the region II = {K ch :0.57 ≤ K ch ≤ 0.92} where 0 < SI < 1 and the layer gets ordered if we further increase the value of K ch as the value of SI turns into zero.   Fig. 5(a,b and c)  To get the complete understanding of the simultaneous effect of K ch and τ, we rigorously plot the states of the upper layer network in the K ch − τ parameter space for the range K ch ∈ [0, 1.5] and τ ∈ [0, 1.0] in Fig. 7. The strength of incoherence (SI) is used to distinguish between different dynamical states by changing K ch and τ simultaneously. Blue, red and yellow colors stand for the region of coherent, chimera, and incoherent states respectively. As K ch increases, the successive scenario of incoherent state, chimera state followed by coherent states remains unaltered for almost all the values of τ in the parameter space. The widening of the region reflecting chimera pattern due to the introduction of time delay τ is observed in the parameter space as well. These proves that the chimera patterns persist as a natural link between incoherence and coherence even in the presence of homogeneous delay in the inter-layer interaction, irrespective of the amount of time-delay. Also we must note that for the case of instantaneous inter-layer chemical synaptic coupling, Fig. 3 shows chimera state starts occurring only when K ch ≥ 1.075. But the introduction of the information transmission delay τ in the inter-layer coupling brings about chimera pattern even when K ch is much smaller than 1.075 which is clear from Fig. 7.
Heterogeneous inter-layer chemical synaptic coupling delays. Next we consider the most general form of interaction between the two layers of oscillators by taking heterogeneous delays in the inter-layer synaptic coupling. Formerly impact of heterogeneous time delay is investigated in the context of different types of cluster synchronization 65,66 , synchrony 67 and amplitude death 68 of coupled neural networks. Information transmission delays from upper to lower layer i.e., τ 1 and from lower to upper layer i.e., τ 2 are different in this case. These heterogeneous time delays can also be reduced to homogeneous time delay τ = τ τ + 2 1 2 , by a time shift transformation, recently proposed by Lucken et al. 69 . But our main concern will be with the emergence of chimera patterns due to the co-action of these different time-delays τ 1 and τ 2 .
In order to do that, first we keep K ch fixed at 0.8, where chimera was observed for a long range of τ in the homogeneous delay case (Fig. 7). We plot the τ 1 − τ 2 parameter space for the range τ 1 ∈ [0, 1.4] and τ 2 ∈ [0, 1.4] in Fig. 8(a). Initially, as τ 1 and τ 2 increase simultaneously, from the state of incoherence the upper layer network may achieve chimera state when τ 1 and τ 2 passes a certain value (satisfying  τ τ + . 0 26 1 2 ). In fact, then we found a region (in black color) where incoherent and chimera pattern coexist in the parameter space. For further increment in τ 1 and τ 2 , coexistence of chimera and coherence is observed (green region) and even higher values of τ 1 and τ 2 leads the upper layer to the state of coherence (blue region). The more the value of τ 1 , less the value of τ 2 is needed (and vice versa ) for the network to attain coherent state as in Fig. 8(a). This is how in this case, chimera may be found in a wide range of the τ 1 − τ 2 parameter space.
Next, taking K ch from the regime where synchrony is observed for a long range of τ (the homogeneous delay), we again discover a convincing enough chimera region in the τ 1 − τ 2 parameter space followed by a synchronous region. For fixed K ch = 1.0, in the τ 1 − τ 2 parameter space with τ 1 ∈ [0, 0.5] and τ 2 ∈ [0, 0.5], disordered state (yellow region) appears only for very small values of τ 1 , τ 2 as shown in Fig. 8(b). But as both the delays increase, the region of coexistence of incoherence and chimera (in black color) and a region of coexistence of coherence and chimera pattern (in green color) develops as in the previous case followed by the coherent state (blue region).
Finally, we fix K ch = 0.35 at such a value where no chimera was observed for almost any value of the homogeneous time-delay τ (Fig. 7). Taking K ch = 0.35, for τ 1 and τ 2 satisfying  τ τ + . 0 66 1 2 , the chimera state can emerge in the uncoupled neurons as in Fig. 8(c). The region (in black color) of coexistence of disordered and chimera states is identified. But coherence has not been observed for almost any value of τ 1 and τ 2 in this case. In fact, even if we fix K ch at other smaller or larger values than the above values taken, we will have the similar type of τ 1 − τ 2 space except only the range of τ 1 and τ 2 at which disordered, chimera and coherent states show up. This is how the information transmission delays can play a crucial role as far as the formation of chimera pattern is  Fig. 5(a,b,c) respectively.
Scientific RepoRts | 6:39033 | DOI: 10.1038/srep39033 concerned. By introducing heterogeneous delay, we tried to make our model as general as possible, however, the effect of both τ 1 and τ 2 are similar as seen in Fig. 8.

Discussion
In summary, we have inspected how multi-layering can bring about chimera states in a network of uncoupled neurons where the multi-layering layer (lower layer) of neurons are globally coupled. The neurons in the multi-layering layer has been assumed to be connected with electrical synapses whereas inter layer connection has been supposed to be of chemical synaptic type. The coaction of these two types of synapses leads the uncoupled layer (upper layer) to a chimera state. We discussed the existence of density dependent threshold for the emergence of chimera states in uncoupled neurons. It is identified that delay in the inter layer coupling may enlarge the range of inter layer coupling strength for which the chimera pattern appears, compared to instantaneous inter layer coupling. We also obtained chimera states when we took other types of coupling in the common medium i.e. the multi-layering layer (see Supplementary Information). Our results therefore seem to be relevant for brain dynamics where coexistence of coherent and incoherent behaviors of the neurons appear.

Methods
In the present work our object is to study the behavior of the uncoupled neurons which are not directly coupled rather they are communicated with each others via a common medium. The neuronal dynamics can be controlled by the coaction of two synapses, namely of electric and chemical synapses.
Hindmarsh-Rose neuronal model. We consider each neuron in the multilayer network with Hindmarsh-Rose neuronal model dynamics. The Hindmarsh-Rose neuronal model is a popular for its chaotic bursting behavior and the original form is as follows: where x-variable represents the membrane potential, y and z represent the transport of ions across the membrane through the fast and slow channels respectively. The variable z corresponds the controls of speed of variation of the slow current and this speed is control by the small parameter c. Here the parameter I denotes an external current that enters the neuron and x 0 controls delaying and advancing the activation of the slow current in the modeled neuron. For the sake of simplicity, after parameter redefinition or linear transformation 70 x → x, y → 1 − y, z → 1 + I + z, d → a + α, e → − 1 − I − bx 0 , the above Eq. (1) can be rewritten as The transformed model (2) is a phenomenological model that gives all the common dynamical features in a number of biophysical modeling studies of various bursting.
We consider Hindmarsh-Rose models as the nodes of the network (schematic diagram is shown in Fig. 1) where both types of synapses (electrical and chemical) are present, the equations governing the dynamics of upper layer becomes and for the lower layer as follows where (x i,1 , y i,1 , z i,1 ) and (x i,2 , y i,2 , z i,2 ) represent the state vectors for the neurons in the upper and lower layers respectively, = ⋅ ⋅ ⋅ i N 1,2, , ; N being the number of neurons in each of the layers of the network. Here τ 1 and τ 2 are the time-delays required to propagate the information from upper to lower layer and lower to upper layer respectively. The variables x i,k represent the membrane potentials, and the variables y i,k and z i,k are the transport of ions across the membrane through the fast and slow channels, respectively for upper and lower layers for k = 1, 2. We consider c a small positive parameter so that z i,k varies much slower than x i,k and y i,k (k = 1, 2). The regular square-wave bursting is observed for the set of parameter values: a = 2.8, α = 1.6, c = 0.001, b = 9, and e = 5 (time series shown in Fig. 2(a) in the Supplementary information). This system is monostable, that is, the coexistence of a stable equilibrium point and a limit cycle has not been observed for this set of parameter values. The synapses are excitatory or inhibitory for the reversal potential v s greater or less than x i,k (t) for all x i,k (t) and all times t. If i-th and j-th neurons are connected through electrical synapses then E(x j,2 , x i,2 ) = x j,2 − x i,2 , i, j = 1, 2, … , N. From physicist's perspective, at electrical synapse, gap junction between two neurons allows electron to move from one to another neuron via intercellular channels. This synapse is bidirectional and of a local character and occurring between those neurons which are spatially very close. By the mutual interaction through these synapses, neurons exhibit coherence or phase synchronization very easily and resulting into a group of synchronized neurons. The coupling strength associated with these synapses is K el . Whereas, the chemical synaptic coupling function Γ (x) is modeled by the sigmoidal nonlinear input-output function as x ( ) s with λ determining the slope of the function and Θ s is the synaptic threshold. There is no such intercellular continuity at chemical synapses and no direct flow of electron from one neuron to another. The space gap between presynaptic and postsynaptic neurons is substantially greater at chemical synapses than electric synapses. Synaptic current flows from presynaptic neuron to postsynaptic neuron only in response to the secretion of neurotransmitters (e.g. acetylchosine, glutamate etc.). This synapse is either excitatory or inhibitory that depends on the neurotransmitters. We assume the chemical synapses are in excitatory for v s = 2.0 as it has a important function in information processing within the brain and throughout the peripheral nervous system. We choose the threshold Θ s = − 0.25 so as to make every spike in the isolated neuron burst to reach the threshold and we fix the value λ = 10. Here K ch is the coupling strength associated with the chemical synapses. We assume, for simplicity, both the synapses, namely electrical and chemical synapses transmit the electron bidirectionally from one to another neuron. The systems (3) and (4) are integrated using fifth-order Runge-Kutta-Fehlberg integration scheme with integration time step 0.01 for non-delayed cases i.e. τ 1 = τ 2 = 0.0. In presence of synaptic coupling delay, we integrate systems (3) and (4) using modified Heun method with integration time step 0.01.
Characterization of chimera state: strength of incoherence. To characterize the disordered, chimera and coherent states, we use the statistical measures using the time series of the network 54 . In order to do that we measure the strength of incoherence (SI) using a local standard deviation analysis. To calculate SI, we firstly define the transformations w i,1 = x i,1 − x i+1,1 , i = 1, 2, ···, N. We divided the total number of neurons in upper layer into M (even) bins of equal length n = N/M and σ 1 (m) is the local standard deviation in each of these bins as follows where Θ (·) is the Heaviside step function and δ is a predefined threshold. Consequently, the values SI = 1, SI = 0 and 0 < SI < 1 correspond to the incoherent, coherent and chimera states respectively.