# Emergent Dynamics and Spatio Temporal Patterns on Multiplex Neuronal Networks

- Department of Physics, Indian Institute of Science Education and Research Tirupati, Tirupati, India

We present a study on the emergence of a variety of spatio temporal patterns among neurons that are connected in a multiplex framework, with neurons on two layers with different functional couplings. With the Hindmarsh-Rose model for the dynamics of single neurons, we analyze the possible patterns of dynamics in each layer separately and report emergent patterns of activity like in-phase synchronized oscillations and amplitude death (AD) for excitatory coupling and anti-phase mixed-mode oscillations (MMO) in multi-clusters with phase regularities when the connections are inhibitory. When they are multiplexed, with neurons of one layer coupled with excitatory synaptic coupling and neurons of the other layer coupled with inhibitory synaptic coupling, we observe the transfer or selection of interesting patterns of collective behavior between the layers. While the revival of oscillations occurs in the layer with excitatory coupling, the transition from anti-phase to in-phase and vice versa is observed in the other layer with inhibitory synaptic coupling. We also discuss how the selection of these spatio temporal patterns can be controlled by tuning the intralayer or interlayer coupling strengths or increasing the range of non-local coupling. With one layer having electrical coupling while the other synaptic coupling of excitatory(inhibitory)type, we find in-phase(anti-phase) synchronized patterns of activity among neurons in both layers.

## 1. Introduction

The complexity underlying the patterns of dynamical behavior in the brain is a fascinating and challenging research area in recent times (Sporns, 2013). The complexity arises not only from a large number of neurons involved but also from the variety and plasticity of connections or interactions among them during any type of neuronal or cognitive activity (Pereda, 2014; Ashwin et al., 2016). The interactions can be electrical via gap junction and excitatory or inhibitory interaction via chemical synapses. The collective behavior or synchronization among a large number of neurons is essential for various neurobiological processes, which mostly appear due to the inter neuronal synaptic interactions (Pikovsky et al., 2001). Also, various brain disorders, such as Alzheimer's disease, schizophrenia, Parkinson's disease, and epilepsy, have been linked to the abnormal patterns of synchronization among the neurons (Uhlhaas and Singer, 2006; Jalili et al., 2007; Knyazeva et al., 2010). The nature of the collective dynamics can have different forms of oscillatory patterns that include in-phase oscillations, anti-phase oscillations, multi-cluster oscillations, etc. (Jalan and Singh, 2016; Pournaki et al., 2019). In addition, coupled neurons also show quiescent states due to suppression of activity or amplitude death (AD) (Saxena et al., 2021).

We find the multiplex framework is ideal for describing the collective dynamics of neurons since an assembly of neurons can have excitatory or inhibitory types of electrical or chemical synaptic interactions (Boccaletti et al., 2014; Verma and Ambika, 2021). Then, analysis can be done with the same set of neurons distributed in different layers, based on the nature of interactions among them. In the present study, we consider the framework of multiplex networks to study the activity patterns that can emerge or get selected when neurons in one layer interact with each other through excitatory synaptic couplings and neurons in the other layer interact with each other through inhibitory synaptic couplings. Equally interesting and realistic is the case where one layer of neurons interact electrically while in the other layer, the interaction is synaptic or chemical of excitatory or inhibitory type. We begin by studying the patterns of collective dynamics in each layer separately and observe how excitatory synaptic coupling induces completely synchronized oscillations and AD, while inhibitory synaptic coupling induces anti-phase synchronized oscillations for local connections and multi-cluster oscillations with relative phase ordering for non-local connections. In this context, we note that anti-phase synchronization is observed in neuronal networks in human and animal brains (de la Iglesia et al., 2000; Ueda et al., 2002; Ohta et al., 2005), climactic networks (Hinnov et al., 2002; Saenko et al., 2002), food web (Vandermeer, 2004), and lasers (Wiesenfeld et al., 1990). We note in multiplex neuronal networks with attractive and repulsive interactions, anti-phase synchronization is reported recently (Chowdhury et al., 2021) and chimera states are found to occur in multilayer networks of neurons (Majhi et al., 2016, 2017, 2019).

When both layers are multiplexed, we find transfer or selection of activity patterns across the layers, with the revival of oscillations from AD state in the first layer and a transition from anti-phase to in-phase in the second layer. Depending on the strength of intralayer coupling, activity patterns corresponding to the stronger interaction get selected and stabilized across the neurons in both layers. When one layer has electrical coupling and the other layer with synaptic coupling, in-phase or anti-phase oscillations are induced depending on whether synaptic coupling is excitatory or inhibitory. These activity patterns have rhythmic dynamics with mixed-mode oscillations (MMO), which are complex periodic forms of activity. We note such MMOs are experimentally observed and analyzed in neurophysiological studies (Del Negro et al., 2002; Desroches et al., 2013; Ghosh et al., 2020). We study the transitions between such patterns of activity and how the relevant parameters can be tuned for a specific pattern to get selected across the layers.

## 2. Multiplex Neuronal Networks

We consider a multiplex network of neurons with two layers, each of them consisting of an ensemble of *N* Hindmarsh-Rose (HR) neurons coupled on a regular ring network. We take the neurons in the first layer (L1) to be interacting with each other with excitatory synaptic coupling and those in the second layer (L2) interacting through an inhibitory synaptic coupling. The neurons in L1 interact with neurons in L2 with multiplex like *i* to *i* coupling via feedback. The dynamics of the multiplex network of neurons is thus modeled as shown in Equation (1)

where ${B}_{i,j}=a{x}_{i,j}^{2}-{x}_{i,j}^{3}-{y}_{i,j}-{z}_{i,j}$, *i* = 1, 2, …, *N* and *j* = 1, 2 (Majhi et al., 2016). The variable *x*_{i,j} represents the action potential, and the variables *y*_{i,j} and *z*_{i,j} represent the transport of ions across the membrane through fast and slow channels, respectively. The function Γ(*x*_{i,j}) = 1/{1 + exp[−β(*x*_{i,j} − ϕ_{s})]} is the sigmoidal chemical synaptic coupling function with *V*_{s} as reversal potential. Here, we take the reversal potential *V*_{s} = 2 such that *V*_{s} > *x*_{i}(*t*) can be satisfied. We choose the synaptic threshold ϕ_{s} = −0.25 and β = 10 in the sigmoidal function. Also, *p*_{1} and *p*_{2} take care of the range of interactions, whether it is local or non-local, with *p*_{1,2} = 1 being local. The other system parameters are *a* = 2.8, α = 1.6, *b* = 9, *c* = 0.001, and *e* = 5 such that the individual HR neurons show regular square-wave bursting dynamics. In the present work, the emergent dynamics of Hindmarsh-Rose (HR)neurons are studied by solving Equation (1), using fourth-order Runge-Kutta method, with initial conditions are chosen randomly between −1 and 1, for various cases as presented below.

### 2.1. Spatio Temporal Patterns on a Single Layer

We begin by considering the emergent dynamics or patterns of activity that can develop in each layer in the absence of multiplexing with ϵ = 0 and number of neurons *N* = 50 in Equation (1).

In layer L1 with excitatory synaptic coupling among neurons, we observe that, for sufficient strength of synaptic coupling, they settle to a completely synchronized oscillatory state, which is shown in Figure 1A at λ_{1} = 1.5. However, the nature of oscillations is changed from intrinsic bursts to varied forms like bursts of decreasing amplitudes and broad spikes as λ_{1} is increased. With stronger coupling, at λ_{1} = 2.9, these spikes are suppressed, and the layer goes to AD. We note AD phenomenon has been reported earlier in globally coupled HR neurons (Prasad et al., 2010). Here, we find that AD can occur for all values of *p*_{1}, local, non-local, and global, with sufficient strength of coupling. To detect the transition to AD, we compute the average amplitude of the spikes of all the neurons using Equation (2) (Verma et al., 2017).

This is plotted in Figure 1B for *p*_{1} = 1 with increasing λ_{1}. We find that the average amplitude increases with λ_{1} initially, reaches a maximum, and then decreases. At λ_{1} = 2.9, there is a sudden transition to AD. The nature of the burst patterns in these regions differs as spikes of decreasing amplitude in each burst that change to square bursts before reaching AD. We repeat the study by increasing N to 100 and 500 and find qualitatively similar results.

**Figure 1. (A)** Time series showing synchronized bursts of action potentials from all the neurons on layer L1, coupled with excitatory synaptic coupling for coupling strength λ_{1} = 1.5, **(B)** Average amplitude of oscillations for increasing coupling strength λ_{1}. The transition from oscillatory state to amplitude death (AD) state occurs at λ_{1} = 2.9. Here, *p*_{1} = 1, ϵ = 0, and *N* = 50.

For the dynamics on the second layer L2 with inhibitory synaptic coupling among neurons, we first study the case when *p*_{2} = 1, i.e., the system has only local interactions. We find that the emergent dynamics in this case shows anti-phase synchronized oscillations, which is clear from the time series, and spatio-temporal plots shown for coupling strength, λ_{2} = 1, in Figures 2a,b. First, we note that the nature of dynamics is changed from intrinsic bursts, in this case also, revealing MMO. Moreover, we find the neurons in one cluster, say at all even number sites, are all synchronized completely but are in anti-phase with those in the other cluster, at odd number sites. This is made more explicit by plotting the time series of all odd number of neurons and even number of neurons separately in Figures 2c,d, that display the pattern of anti-phase synchronized oscillations among the adjacent neurons.We also show the time series of the other two variables *y*_{i} and *z*_{i} in Figures 2e,f, respectively.

**Figure 2. (a)** Mixed-mode oscillations (MMO) of action potentials and **(b)** spatio temporal plot of neurons coupled with inhibitory synaptic coupling on layer L2 at coupling strength λ_{2} = 1. **(c)** Time series of all the odd number nodes *x*_{2i−1}, showing complete synchronized oscillations, for *i* = 1;2 : : : *n*/2 and **(d)** time series of all the even number nodes *x*_{2i}, where *i* = 1;2 : : : *n*/2 that indicate completely synchronized oscillations among them. Time series of the other variables *y*_{i} and *z*_{i} are plotted in **(e,f)**, respectively. Here, other parameters are kept as *p*_{2} = 1, ϵ = 0, and *N* = 50. The anti-phase nature of the oscillations in adjacent nodes and in-phase nature in alternate nodes separates the network into two clusters. The color bar in **(b)** (also in spatio-temporal plots in later figures) indicates the values of action potential *x*_{i}.

For a detailed characterization of the observed phase order in temporal dynamics, we calculate the phase of each neuron from its time series, *x*_{i}. We note the time ${T}_{k}^{i}$, (*k* = 1, 2, …) at which *x*_{i} crosses the chosen threshold value, and then, we calculate the phase of the *i*^{th} neuron using the following equation (Pikovsky et al., 2001):

where *i* = 1, 2, …, *N*. In Figure 3A, the phase of each neuron calculated relative to that of the first neuron is plotted. It is clear that every odd neuron is in anti-phase with every even neuron. The snapshot of *x*_{i} at a given time τ is shown in Figure 3B, which further confirms the anti-phase pattern of the mixed mode oscillations. This is induced by the range (nearest neighbor) and the nature (inhibitory) of the coupling chosen in this context. Thus, the neurons in effect form two clusters such that synchronized oscillations in one cluster are anti-phase with that in the other cluster. Further, we calculate the spike frequency of the large amplitude oscillations of *i*^{th} neuron as shown in Equation (4) (Mozumdar and Ambika, 2019):

where *K*_{i} refers to the number of spikes for the *i*^{th} neuron in each burst and ${t}_{k}^{i}$ corresponds to time of the maximum of the *k*^{th} spike. Then, the average frequency obtained from this, is plotted in Figure 3C with increasing coupling strength λ_{2}. Here, we can see that the average frequency increases with increasing λ_{2}. We also show how the average amplitude < *A* > of coupled neurons increases with λ_{2}, for the range considered as shown in Figure 3D. Both the frequency and amplitude calculated here relate to the large amplitude spikes of the mixed-mode oscillatory states of the neurons. We note such activity patterns of synchronized oscillations with amplification are reported in multiplex networks in a different context (Njougouo et al., 2020).

**Figure 3. (A)** Phase of each neuron in layer L2, coupled with inhibitory synaptic coupling calculated relative to that of its first neuron for coupling strength λ_{2} = 1 and *p*_{2} = 1. **(B)** Snapshot of *x*_{i} at a given time τ, shows that every odd neuron is in anti-phase with every even neuron. **(C)** Average frequency and **(D)** average amplitude of the large amplitude oscillations of the MMO for increasing coupling strength λ_{2}.

As the range of coupling increases or the coupling becomes non-local, we observe traveling wave-like patterns. In Figure 4A, we plot the time series of the action potential from node 1 and node 2 and in Figure 4B from node 1 and node 4. It is clear that nodes 1 and 4 are almost synchronized but with a small phase shift. We find this shift in the phase depends on the coupling strength and the range of coupling. The spatio-temporal plot, in this case, shows traveling wave like patterns, as shown in Figures 4C,D, for *p*_{2} = 2 and *p*_{2} = 5, and λ_{2} = 3. For larger sizes of networks also, we find qualitatively similar emergent dynamics.

**Figure 4**. Time series of action potentials from **(A)** nodes 1 and 2 and **(B)** nodes 1 and 4 in layer L2 coupled with inhibitory synaptic coupling for *p*_{2} = 2. We find nodes *i* and *i* + 3 show inphase synchronized oscillations with small phase shift. The spatio-temporal plot of neurons showing traveling wave like patterns **(C)** for *p*_{2} = 2 **(D)** for *p*_{2} = 5. Here λ = 3 and *N* = 50.

### 2.2. Dynamics of the Multiplex Network of Neurons With Excitatory and Inhibitory Synaptic Couplings

With the two layers of neurons multiplexed, we study how different emergent activity patterns of dynamics get selected across the layers as parameters are varied. We first consider the case where neurons of layer L1 are uncoupled, while those of layer L2 are coupled with inhibitory synaptic coupling and both layers are coupled to each other via *i* to *i* connections with feedback coupling of strength ϵ as given in Equation (1). In this case, with *p*_{2} = 1, λ_{2} = 6 for L2, ϵ = 1, the patterns of synchronized oscillations that are anti-phase for adjacent nodes on second layer L2, get selected as such in the first layer L1 also. This is clear from Figures 5A1–B2, where the time series and spatio-temporal plots of both layers are given. Also for *p*_{2} = 2 and λ_{2} = 10, both layers show traveling wave-like oscillations (Figures 5C1–D2). Thus, the emergent dynamics and the corresponding activity patterns get transferred from one layer to other layer when the layers are multiplexed.

**Figure 5**. Transfer of dynamical patterns from L2 to L1, in the 2-layer multiplex network with neurons on L2 coupled with inhibitory synaptic coupling and neurons in L1 uncoupled. Time Series of action potential and spatio-temporal plot are shown for different values of *p*_{2} and λ_{2}: **(A1,B1)** first layer and **(A2,B2)** second layer with *p*_{2} = 1 and λ_{2} = 6. **(C1,D1)** first layer and **(C2,D2)** second layer, with *p*_{2} = 2 and λ_{2} = 10. Here λ_{1} = 0, ϵ = 2, and *N* = 50.

Next, we consider the neurons of first layer L1 coupled with excitatory synaptic coupling, with neurons of L2 still coupled with inhibitory synaptic coupling, both with local couplings as *p*_{1} = 1, and *p*_{2} = 1. With the interlayer coupling strength at ϵ = 1, for λ_{1} = 0.1, and λ_{2} = 4, we observe that both layers L1 and L2 exhibit anti-phase synchronized oscillations, with phase ordering which is shown in Figures 6A1,B1, respectively. When we set λ_{1} = 3.0 and λ_{2} = 0.1, we observe in-phase synchronized oscillations in both layers, which is shown in Figures 6A2,B2, respectively. Thus, we see that for strong inhibitory synaptic coupling strength, both layers show anti-phase synchronized oscillations in adjacent nodes, while for strong excitatory synaptic coupling, both layers show in-phase synchronized oscillations.

**Figure 6**. Time series of the action potentials of multiplex HR neurons for both layers L1 (left panel) and L2 (right panel) for different values of λ_{1} and λ_{2}: **(A1,B1)** λ_{1} = 0.1 and λ_{2} = 1.0, show anti-phase synchronized oscillations in two clusters, and **(A2,B2)** λ_{1} = 3 and λ_{2} = 0.1, show in-phase synchronized oscillations. Here ϵ = 1, *p*_{1} = 1, *p*_{2} = 1, and *N* = 50. The pattern of the dynamics on the layer of larger intralayer coupling strength is selected across both layers.

Also, as couplings become non-local, with *p*_{2} = 2 and 3, both layers show phase shifted oscillations and spatio-temporal dynamics that are transferred from L2 to L1 for larger λ_{2}. We also observe that these states are selected by layer 1 for all values of *p*_{1} up to *p*_{1} = 10, λ_{1} = 0.1. The spatio-temporal plots for *p*_{2} = 2, and λ_{2} = 6, shown in Figures 7A1,B1, and for *p*_{2} = 3 and λ_{2} = 10, in Figures 7A2,B2, indicate the transfer of dynamical patterns across the layers. However, the nature of spikes and bursts in layers L1 and L2 differs due to difference in the parameter chosen. So, the selection of the specific activity patterns on both layers depends on the relative intra-layer coupling strengths and follows the spatio-temporal dynamics of the layer with larger intra-layer coupling strength. This is further illustrated for other types of emergent dynamics below.

**Figure 7**. Spatio-temporal plots of multiplex HR neurons for layer L1 (left panel) and L2 (right panel) for different values of *p*_{2}, andλ_{2}: **(A1,B1)** *p*_{2} = 2 and λ_{2} = 6, **(A2,B2)** *p*_{2} = 3, and λ_{2} = 10. Here, λ_{1} = 0.1, ϵ = 1, *p*_{1} = 10, and *N* = 50. The spatio-temporal dynamics of the layer with larger coupling strength gets selected in both the layers. However, the nature of spikes and bursts are different in L1 and L2 due to difference in intralayer parameters.

As reported earlier, when ϵ = 0, both L1 and L2 function as independent layers, and for higher synaptic coupling strength λ_{1} = 3, layer L1 goes to AD and at λ_{2} = 0.3, layer L2 shows anti-phase synchronized oscillations in two clusters (Figures 8A1,B1). But when both layers are multiplexed with ϵ = 1, we observe a revival of oscillations from death state on layer L1 and in-phase oscillation on layer L2, as shown in Figures 8A2,B2, respectively. Further, we observe that the activity pattern in L2 undergoes a transition from in-phase to anti-phase as λ_{2} is tuned. This transition from in-phase to anti-phase with an increase in λ_{2} in layer L2 is shown in Figure 9A, where the average phase difference is calculated as $<\varphi >=\frac{1}{N}\sum _{i=1}^{N}({\varphi}_{i}-{\varphi}_{i+1})$, with ϕ_{i} obtained for each neuron from Equation (3). The inhibitory synaptic coupling in one layer can revive the oscillations from the suppressed state on the other layer. The variety of interesting activity patterns of spatio-temporal dynamics and their selection across layers happens at low to moderate values of interlayer coupling strengths. When the interlayer coupling strength ϵ is increased, to say ϵ = 10, both layers settle to AD states (Figures 9B,C), and the time series near the transition point is as shown in Figure 9D. Thus, the selection of activity patterns in both layers due to multiplexing depends on the nature and strengths of intralayer and interlayer couplings, and therefore, the coupling strengths and range of couplings can be tuned to select any desired pattern of activity.

**Figure 8**. Revival of activity in layer L1 due to multiplexing with L2 and transition from in phase to anti-phase in L2 as parameters are tuned. Time series of action potentials for multiplex HR neurons in layer L1 (left panel) and L2 (right panel) are plotted for different values of λ_{1} and λ_{2}: **(A1)** at λ_{1} = 3, ϵ = 0 neurons in layer L1 exhibit AD, **(B1)** at λ_{2} = 0.3, and ϵ = 0: neurons in layer L2 show anti-phase oscillations in two clusters: **(A2,B2)** λ_{2} = 0.3 and ϵ = 1: observed revival of oscillations in L1 and transition from anti-phase to in-phase in L2. Here, λ_{1} = 3, *p*_{1} = 1, *p*_{2} = 1, and *N* = 50.

**Figure 9. (A)** Transition from in-phase to anti-phase for synchronized activity in L2. Here, average phase of neurons on the second layer for varying coupling strength λ_{2}, is shown with neurons in L1 coupled locally with *p*_{1} = 1 and λ_{1} = 3. Suppression of activity on increasing the interlayer coupling strength. Average amplitudes of neurons on **(B)** L1 and **(C)** L2 are shown with varying interlayer coupling strength ϵ. The nature of spikes and bursts near the transition for ϵ = 8.1 is shown in **(D)**. Here, λ_{1} = 1, λ_{2} = 1.0, *p*_{1} = 1, *p*_{2} = 1, and *N* = 50.

### 2.3. Dynamics of the Multiplex Network of Neurons With Electrical and Synaptic Coupling

Now, we consider the case where neurons in the first layer (L1) interact with each other with electrical coupling and those in the second layer (L2) interact through synaptic coupling. The dynamics of the multiplex network of neurons thus modeled is given as follows,

Here, we define a parameter E whose sign decides the nature of synaptic coupling, for *E* = 1 neurons in L2 are coupled with excitatory synaptic coupling, and for *E* = −1, second layer are coupled with inhibitory synaptic coupling.

With *E* = 1, and the excitatory coupling strength at λ_{2} = 5, we observe that the coupled system shows in-phase synchronized oscillations, in both layers L1 and L2, as shown in Figures 10A1,B1, respectively. Next, with *E* = −1, the coupled system shows anti-phase synchronized oscillations in both layers L1 and L2 (Figures 10A2,B2). Further, we also observe that along with the transfer of the emergent phenomena from one layer to another, the node of layer L1 shows in-phase synchronization with the same node of layer L2. To indicate this, we show the time series of the 5^{th} node of both layers L1 and L2, where layer L1 coupled with electrical coupling and L2 coupled with excitatory coupling in Figure 11A and L2 coupled inhibitory synaptic coupling in Figure 11B, respectively. Also, for strong electrical coupling strength (λ_{1}) and weak synaptic coupling (λ_{2}) (inhibitory or excitatory), we observe traveling wave patterns on both layers.

**Figure 10**. Time series of action potentials for multiplex HR neurons from layer L1 (left panel) to L2 (right panel): **(A1,B1)** in-phase patterns for neurons of L1 coupled with electrical coupling and that of L2 coupled with excitatory synaptic coupling. **(A2,B2)** Anti-phase activity for neurons on L1 coupled with electrical coupling and that on L2 coupled with inhibitory synaptic coupling. The dynamics of L2 are transferred to L1 in both cases. λ_{1} = 0.5, λ_{2} = 5, *p*_{1} = 1, *p*_{2} = 1, and *N* = 50.

**Figure 11**. In-phase synchronized dynamics between similar nodes on layers L1 and L2. Time series of node 5 from both layers L1 and L2 are plotted with **(A)** neurons of L1 coupled with electrical coupling and that of L2 coupled with excitatory synaptic coupling; **(B)** neurons of L1 coupled with electrical coupling and that of L2 coupled with inhibitory synaptic coupling. The other parameters are λ_{1} = 0.5, λ_{2} = 5, *p*_{1} = 1, *p*_{2} = 1, and *N* = 50.

## 3. Conclusion

In this study, we report the selection of various activity patterns as the emergent spatio-temporal dynamics on a multiplex neuronal network of HR neurons where the nature of interaction in each layer can be different. This framework can thus model the plasticity and variability of connections among neurons which can exist as synaptic or electrical in nature with excitatory or inhibitory connections. By tuning the strengths of connections in each layer and across layers, the network can select various activity patterns and induce the pattern from one layer to the other.

We first present the pattern of dynamics on the first layer L1, where neurons are coupled through excitatory synaptic couplings. By tuning the synaptic coupling strength, the coupled neurons can be in completely synchronized oscillations, while for strong synaptic coupling strength, the oscillations are suppressed to the state of AD. The phenomenon of AD is observed for all values of *p*_{1}, corresponding to the local, non-local, and global types of couplings. The second layer of neurons, coupled with inhibitory synaptic coupling, shows anti-phase synchronized oscillations with amplification when the neurons are locally coupled, i.e., *p*_{2} = 1. The anti-phase synchronized oscillations are interesting in two aspects. First, the nature of oscillations are MMO with enhanced frequency and amplitude with large amplitude spikes, and second, the phase relationship among them occurs in an orderly way, with alternate neurons being in phase and neighboring ones being in anti-phase. Thus, the whole network splits into two clusters, every odd node belonging to one cluster and every even node to the other cluster. For *p*_{2} = 2 and 3, we get traveling wave type of oscillations over the network.

When the two layers are multiplexed, for sufficient inhibitory coupling strength, we observe mixed-mode synchronized oscillations that are phase-shifted get selected on both layers. In general, the selection of the specific pattern of activity on both layers can be controlled by tuning the relative intra-layer coupling strengths.

Also, multiplexing can revive the oscillations from the AD state on the first layer by changing the inhibitory coupling strength on the second layer. We also report the transition from anti-phase to the in-phase type of MMO, and vice versa that get selected as the excitatory and inhibitory coupling strengths are tuned to specific values. We repeat the study by increasing the size of the networks in both layers to 100 and 500 and find qualitatively similar results.

With the nature of coupling among neurons in one layer L1 electrical, while the other layer L2 has neurons with synaptic connections, we observe in-phase synchronized activity in both layers when L2 has excitatory connections and anti-phase activity when it has inhibitory connections. We also find neurons at similar nodes in both layers are synchronous with in-phase oscillations.

We note the variety of activity patterns presented here that occur for a collection of neurons forming a multiplex network, corresponding to experimentally observed patterns of activity reported recently (Crofts et al., 2016). Also, modulation of neuronal oscillation frequency is reported to occur during sensory information processing (Lee et al., 2018). Thus, the study provides a better understanding of the mechanism underlying such patterns known to occur in brain networks that incorporate multiplex network architecture naturally (Frolov et al., 2020).

## Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

## Author Contributions

GA and UK conceived and designed the study. UK performed the numerical analysis and drafted the manuscript. All authors edited and approved the manuscript.

## Conflict of Interest

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.

## Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

## References

Ashwin, P., Coombes, S., and Nicks, R. (2016). Mathematical frameworks for oscillatory network dynamics in neuroscience. *J. Math. Neurosci.* 6:2. doi: 10.1186/s13408-015-0033-6

Boccaletti, S., Bianconi, G., Criado, R., Del Genio, C. I., Gomez-Gardenes, J., Romance, M., et al. (2014). Explosive transitions in complex networks' structure and dynamics: percolation and synchronization. *Phys. Rep.* 544, 1–122. doi: 10.1016/j.physrep.2014.07.001

Chowdhury, S. N., Rakshit, S., Buldu, J. M., Ghosh, D., and Hens, C. (2021). Antiphase synchronization in multiplex networks with attractive and repulsive interactions. *Phys. Rev. E* 103:032310. doi: 10.1103/PhysRevE.103.032310

Crofts, J. J., Forrester, M., and O'Dea, R. D. (2016). Structure-function clustering in multiplex brain networks. *EPL* 116:18003. doi: 10.1209/0295-5075/116/18003

de la Iglesia, H. O., Meyer, J., Carpino, A. Jr., and Schwartz, W. J. (2000). Antiphase oscillation of the left and right suprachiasmatic nuclei. *Science* 290, 799–801. doi: 10.1126/science.290.5492.799

Del Negro, C. A., Wilson, C. G., Butera, R. J., Rigatto, H., and Smith, J. C. (2002). Periodicity, mixed-mode oscillations, and quasiperiodicityin a rhythm-generating neural network. *Biophys. J.* 82:206. doi: 10.1016/S0006-3495(02)75387-3

Desroches, M., Kaper, T. J., and Krupa, M. (2013). Mixed-mode bursting oscillations: dynamics created by a slow passage through spike-adding canard explosion in a square-wave burster. *Chaos* 23:046106. doi: 10.1063/1.4827026

Frolov, N., Maksimenko, V., and Hramov, A. (2020). Revealing a multiplex brain network through the analysis of recurrences. *Chaos* 30:121108. doi: 10.1063/5.0028053

Ghosh, S., Mondal, A., Ji, P., Mishra, A., Dana, S. K., Antonopoulos, C. G., et al. (2020). Emergence of mixed mode oscillations in random networks of diverse excitable neurons: the role of neighbors and electrical coupling. *Front. Comput. Neurosci.* 14:49. doi: 10.3389/fncom.2020.00049

Hinnov, L. A., Schulzb, M., and Yiouc, P. (2002). Interhemispheric space-time attributes of the Dansgaard-Oeschger events between 100 and 0 ka. *Quat. Sci. Rev.* 21:1228. doi: 10.1016/S0277-3791(01)00140-8

Jalan, S., and Singh, A. (2016). Cluster synchronization in multiplex networks. *EPL* 113:30002. doi: 10.1209/0295-5075/113/30002

Jalili, M., Lavoie, S., Deppen, P., Meuli, R., Do, K. Q., Cuenod, M., et al. (2007). Resiliency of EEG-based brain functional networks. *PLoS ONE* 2:e1059. doi: 10.1371/journal.pone.0001059

Knyazeva, M. G., Jalili, M., Brioschi, A., Bourquin, I., Fornari, E., Hasler, M., et al. (2010). Topography of EEG multivariate phase synchronization in early Alzheimer's disease. *Neurobiol. Aging*. 31:1132. doi: 10.1016/j.neurobiolaging.2008.07.019

Lee, B., Shin, D., Gross, S. P., and Cho, K. H. (2018). Combined positive and negative feedback allows modulation of neuronal oscillation frequency during sensory processing, *Cell Rep.* 25, 1548–1560. doi: 10.1016/j.celrep.2018.10.029

Majhi, S., Bera, B. K., Ghosh, D., and Perc, M. (2019). Chimera states in neuronal networks: a review. *Phys. Life Rev.* 28, 100–121. doi: 10.1016/j.plrev.2018.09.003

Majhi, S., Perc, M., and Ghosh, D. (2016). Chimera states in uncoupled neurons induced by a multilayer structure. *Sci. Rep.* 6:39033. doi: 10.1038/srep39033

Majhi, S., Perc, M., and Ghosh, D. (2017). Chimera states in a multilayer network of coupled and uncoupled neurons. *Chaos* 27:073109. doi: 10.1063/1.4993836

Mozumdar, K., and Ambika, G. (2019). Frequency locking and travelling burst sequences in community structured network of inhibitory neurons with differing time-scales. *Commun. Nonlinear Sci. Num. Simulat.* 69:320. doi: 10.1016/j.cnsns.2018.09.026

Njougouo, T., Camargo, V., Louodop, P., Ferreira, F. F., Talla, P. K., and Cerdeira, H. A. (2020). Dynamics of multilayer networks with amplification. *Chaos* 30:123136. doi: 10.1063/5.0025529

Ohta, H., Yamazaki, S., and McMahon, D. G. (2005). Constant light desynchronizes mammalian clock neurons. *Nat. Neurosci.* 8:267. doi: 10.1038/nn1395

Pereda, A. E. (2014). Electrical synapses and their functional interactions with chemical synapses. *Nat. Rev. Neurosci.* 15, 250–263. doi: 10.1038/nrn3708

Pikovsky, A. S., Rosenblum, M. G., and Kurths, J. (2001). *Synchronization: A Universal Concept in Nonlinear Sciences*. Cambridge: Cambridge University Press. doi: 10.1017/CBO9780511755743

Pournaki, A., Merfort, L., Ruiz, J., Kouvaris, N. E., Havel, P., and Hizanidis, J. (2019). Synchronization patterns in modular neuronal networks: a case study of *C. elegans*. *Front. Appl. Math. Stat.* 5:52. doi: 10.3389/fams.2019.00052

Prasad, A., Dhamala, M., Adhikari, B. M., and Ramaswamy, R. (2010). Amplitude death in nonlinear oscillators with nonlinear coupling. *Phys. Rev. E* 81:027201. doi: 10.1103/PhysRevE.81.027201

Saenko, A. O., Schmittner, A., and Weaver, A. J. (2002). The Atlantic–Pacific Seesaw. *J. Clim.* 17, 2033–2038. doi: 10.1175/1520-0442(2004)017<2033:TAS>2.0.CO;2

Saxena, G., Prasad, A., and Ramaswamy, R. (2012). Amplitude death: the emergence of stationarity in coupled nonlinear systems. *Phys. Rep.* 521:205. doi: 10.1016/j.physrep.2012.09.003

Sporns, O. (2013). Structure and function of complex brain networks. *Dialog. Clin. Neurosci.* 15, 247–262. doi: 10.31887/DCNS.2013.15.3/osporns

Ueda, H. R., Chen, W., Adachi, A., Wakamatsu, H., Hayashi, S., Takasugi, T., et al. (2002). A transcription factor response element for gene expression during circadian night. *Nature* 418:534. doi: 10.1038/nature00906

Uhlhaas, P. J., and Singer, W. (2006). Neural synchrony in brain disorders: relevance for cognitive dysfunctions and pathophysiology. *Neuron* 52, 155–168. doi: 10.1016/j.neuron.2006.09.020

Vandermeer, J. (2004). Coupled oscillations in food webs: balancing competition and mutualism in simple ecological models. *Am. Nat.* 163:857. doi: 10.1086/420776

Verma, U. K., and Ambika, G. (2021). Tipping induced by multiplexing on two-layer networks. *Eur. Phys. J. Spec. Top.* 230, 3299–3309. doi: 10.1140/epjs/s11734-021-00116-x

Verma, U. K., Sharma, A., Kamal, N. K., Kurths, J., and Shrimali, M. D. (2017). Explosive death induced by mean–field diffusion in identical oscillators. *Sci. Rep.* 7:7936. doi: 10.1038/s41598-017-07926-x

Keywords: multiplex network, neuronal network, synchronization, multi-cluster synchronization, mixed-mode oscillations

Citation: Kumar Verma U and Ambika G (2021) Emergent Dynamics and Spatio Temporal Patterns on Multiplex Neuronal Networks. *Front. Comput. Neurosci.* 15:774969. doi: 10.3389/fncom.2021.774969

Received: 13 September 2021; Accepted: 01 November 2021;

Published: 02 December 2021.

Edited by:

Dibakar Ghosh, Indian Statistical Institute, IndiaReviewed by:

Soumen Majhi, Bar-Ilan University, IsraelSajad Jafari, Amirkabir University of Technology, Iran

Copyright © 2021 Kumar Verma and Ambika. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: G. Ambika, g.ambika@labs.iisertirupati.ac.in