Impact Factor 1.821

The Frontiers in Neuroscience journal series is the 1st most cited in Neurosciences

Original Research ARTICLE

Front. Comput. Neurosci., 18 November 2014 | https://doi.org/10.3389/fncom.2014.00145

Macroscopic complexity from an autonomous network of networks of theta neurons

Tanushree B. Luke, Ernest Barreto and Paul So*
  • School of Physics, Astronomy, and Computational Sciences and The Krasnow Institute for Advanced Study, George Mason University, Fairfax, VA, USA

We examine the emergence of collective dynamical structures and complexity in a network of interacting populations of neuronal oscillators. Each population consists of a heterogeneous collection of globally-coupled theta neurons, which are a canonical representation of Type-1 neurons. For simplicity, the populations are arranged in a fully autonomous driver-response configuration, and we obtain a full description of the asymptotic macroscopic dynamics of this network. We find that the collective macroscopic behavior of the response population can exhibit equilibrium and limit cycle states, multistability, quasiperiodicity, and chaos, and we obtain detailed bifurcation diagrams that clarify the transitions between these macrostates. Furthermore, we show that despite the complexity that emerges, it is possible to understand the complicated dynamical structure of this system by building on the understanding of the collective behavior of a single population of theta neurons. This work is a first step in the construction of a mathematically-tractable network-of-networks representation of neuronal network dynamics.

1. Introduction

The brain is a complex hierarchical network of networks (Zhou et al., 2006; Bullmore and Sporns, 2009; Meunier et al., 2010). Neurons are organized into different neuronal assemblies, and these neuronal assemblies interact with each other, forming larger assemblies (Sherrington, 1906; Hebb, 1949; Harris, 2005). But while there is a wealth of knowledge on the microscopic scale regarding the dynamics of individual neurons, the macroscopic behavior of such interacting populations of neurons is not well understood. Indeed, the functional and information-processing activity of the brain, from perception to consciousness, is thought to result from the emergent collective behavior of these assemblies.

In recent years, the mathematical study of networks of this kind, based on globally-coupled populations of simple phase oscillators, has advanced significantly. This is in large part due to new analytical techniques (Ott and Antonsen, 2008, 2009; Marvel et al., 2009; Ott et al., 2011; Pikovsky and Rosenblum, 2011). These techniques enable the derivation of low-dimensional dynamical systems that reveal the collective emergent behavior of the full discrete population (in the limit of an infinite number of interacting elements). In the context of computational neuroscience, these methods were applied to autonomous globally-coupled networks of canonical Type-I neurons (i.e., theta neurons) by Luke et al. (2013), and to non-autonomous theta neuron networks by So et al. (2014). More recently, Laing (2014) extended these results to include space-dependent coupling. A similar approach, based on phase-response curves, was pursued by Pazó and Montbrió (2014).

Of course, such networks lack the intricate connectivity found in real biological networks. Nevertheless, they are ideal building blocks for the construction of a more realistic, yet mathematically tractable, network-of-networks representation of the brain. In the current study, we consider the simplest hierarchical structure as a first step in this process. Using two globally-coupled networks of theta neurons, we arrange for the activity of one population to drive the second population. Thus, the overall network has an autonomous driver-response configuration. We demonstrate that even in this simplest network-of-networks, the collective behavior of the response network can exhibit a full range of complex behavior, from simple collective rhythms to temporally chaotic dynamics. Most importantly, we provide a complete non-linear dynamical analysis of this system, including predictive bifurcation diagrams for the behavior of the response population in terms of the driver's dynamics and the network characteristics.

2. Recap of Single Population Results

2.1. The Theta Neuron

Neurons are typically classified into two types, based on the nature of the onset of spiking as a constant injected current exceeds an effective threshold (Hodgkin, 1948; Ermentrout, 1996; Izhikevich, 2007). Type-I neurons begin to spike at an arbitrarily low rate, whereas Type-II neurons spike at a non-zero rate as soon as the threshold is exceeded. Neurophysiologically, excitatory pyramidal neurons are often of Type-I, and fast-spiking inhibitory interneurons are often of Type-II (Nowak et al., 2003; Tateno et al., 2004). Near the onset of spiking, Type-I neurons can be represented by a canonical phase model that features a saddle-node bifurcation on an invariant cycle, or SNIC bifurcation (Ermentrout and Kopell, 1986; Ermentrout, 1996). This model has come to be known as the theta neuron, and is given by

θ˙=(1cosθ)+(1+cosθ)η,    (1)

where θ is a phase variable on the unit circle and η is a bifurcation parameter related to the injected current. For η < 0, the neuron is attracted to a stable equilibrium which represents the resting state. An unstable equilibrium is also present, representing the threshold. If an external stimulus pushes the neuron's phase across the unstable equilibrium, θ will move around the circle and approach the resting equilibrium from the other side. When θ crosses θ = π, the neuron is said to have spiked. Thus, for η < 0, the neuron is excitable. As the parameter η increases, these equilibria approach each other and merge via the SNIC bifurcation at η = 0. At this point, the equilibria disappear, leaving a limit cycle. The neuron spikes regularly for η > 0. In the following, we call η the “excitability parameter.”

2.2. A Network of Theta Neurons

We formulate a single population of N theta neurons as follows:

θj.=(1cosθj)+(1+cosθj) [ηj+Isyn],    (2)

where j = 1, …, N is the index for the j-th neuron. The neurons are coupled via a pulse-like synaptic current

Isyn=kNi=1NPn(θi),    (3)

where Pn(θ) = an (1 − cosθ)n, n ∈ ℕ, and an is a normalization constant1 such that

02πPn(θ)dθ=2π.

The parameter n defines the sharpness of the pulse-like synapse in that Pn(θ) becomes more and more sharply peaked as n increases. We assume that the synaptic strength k is the same for all neurons.

Note that the connectivity described by Equations (2) and (3) includes self-coupling terms. These have negligible effect on the collective network dynamics (data not shown), which is to be expected since they represent only one out of N inputs to any given neuron. Nevertheless, we note that these self-connections have real-world analogs in “autapses,” which have been found in several regions of the brain (e.g., Bacci et al., 2003; Bekkers, 2003).

Neurons in real biological networks exhibit a range of different intrinsic dynamics. We model this by taking the excitability parameter ηj of each neuron to be different, with each ηj being drawn randomly from a distribution g(η). In the following analysis, we assume a Lorentzian distribution,

g(η)=1πΔ(ηη0)2+Δ2,    (4)

where η0 is the center of the distribution, and Δ, the half-width at half-maximum, describes the degree of heterogeneity in the population.

2.3. Reduction and Asymptotic States of the Single Population

The macroscopic behavior of our network can be quantified by the “macroscopic mean field,” or order parameter, defined as

z˜(t)=j=1Neiθj,    (5)

where the tilde indicates that the sum is over a finite population of N oscillators. (Below we will drop the tilde in the case of an infinite network.) The magnitude of the order parameter |z˜(t)| ∈ [0, 1] quantifies the degree of synchronization present at time t.

In Luke et al. (2013), we used the Ott-Antonsen method (Ott and Antonsen, 2008, 2009; Ott et al., 2011) to derive a low-dimensional dynamical system whose asymptotic dynamics can be shown to coincide with that of the order parameter of the single-population network defined above (Equations 2–4), in the limit N → ∞. This reduced dynamical system is

z˙=i(z1)22+(z+1)22{Δ+i[η0+kHn(z)]},    (6)

where

Hn(z)=Isyn/k=an(A0+q=1nAq(zq+z*q)),    (7)
     Aq=j,m=0nδj2m,qQjm,    (8)

and

Qjm=(1)j2mn!2jm!(nj)!(jm)!.    (9)

In these equations, z* denotes the complex conjugate of z, and δi,j is the Kronecker delta function on the indices (i,j). Note that Hn(z) = H*n(z) is a real-valued function.

The analysis of Equations (6–9) reported in Luke et al. (2013) showed that the theta neuron network can exhibit three types of asymptotic states. These correspond to a node, a focus, and a limit cycle in the order parameter. A complete bifurcation analysis describing how these states change as the parameters k, η0, and Δ change was also reported. For our purposes in the current work, we now briefly describe the three possible collective macroscopic states.

We called the node, focus, and limit cycle solutions the “Partially Synchronous Rest” (PSR), “Partially Synchronous Spiking” (PSS), and “Collective Periodic Wave” (CPW) states, respectively. In the PSR state, most neurons remain at rest, while in the PSS state, most neurons spike continuously. Nevertheless, in both these states, the macroscopic mean field (or order parameter) sits at an equilibrium. In contrast, the CPW state corresponds to periodic oscillations of the complex order parameter, and typically, both |z(t)| and arg (z) oscillate in time indicating that the individual neurons clump together and spread apart in a periodic fashion. We refer the interested reader to Luke et al. (2013) for further details, including movies that illustrate both the microscopic and macroscopic behaviors of these collective states.

3. Formulation of the Driver-Response Network

In this work, we are interested in the dynamics exhibited by a network of two coupled populations of theta neurons. We formulate the general case, but restrict analysis to the simplest such configuration: a driver-response network.

3.1. General Two-Population Model

Extending the model described above, a general formulation of a pair of interacting populations of theta neurons can be expressed as follows:

θ˙1,j=1+η1,j(1η1,j)cosθ1,j+an(1+cosθ1,j)           [k11N1p=1N1(1cosθ1,p)n+k12N2q=1N2(1cosθ2,q)n],θ˙2,j=1+η2,j(1η2,j)cosθ2,j+an(1+cosθ2,j)           [k21N1p=1N1(1cosθ1,p)n+k22N2q=1N2(1cosθ2,q)n],    (10)

where θ1,j and θ2,j denote the jth neuron in the first and second populations, respectively, and the extension to any number of interacting populations is straightforward. The excitability parameters η1,j and η2,j are randomly drawn from two independent Lorentzian distributions as in Equation (4), with medians η1, η2 and widths Δ1, Δ2, respectively. We take the sharpness parameter of the pulse-like synaptic interaction, n, to be the same for both populations. Macroscopic mean field parameters z˜1(t), z˜2(t) can be defined for each population by analogy with Equation (5).

Adapting the procedures described in Luke et al. (2013), we derived the Ott-Antonsen reduction of the coupled networks of Equation (10). This resulted in the following dynamical system:

z˙1=i(z11)22+(z1+1)22         {Δ1+i[η1+k11Hn(z1)+k12Hn(z2)]},z˙2=i(z21)22+(z2+1)22         {Δ2+i[η2+k21Hn(z1)+k22Hn(z2)]}.    (11)

with Hn(z) defined as in Equations (7–9). As before, the asymptotic dynamics of Equation (11) can be shown to coincide with that of the order parameters of the populations in the network of Equation (10), in the limit N1, N2 → ∞.

We showed in Luke et al. (2013) that the dynamical structure of the single population depends rather weakly on the synaptic sharpness parameter n. Furthermore, we argued that a modest sharpness is more biophysically plausible than the δ-function coupling obtained in the limit n → ∞. Thus, from here on, we fix n = 2 and drop the subscript on Hn to ease notation.

3.2. The Driver-Response System

To put our network in the driver-response form, we set k12 = 0, so that population 1 receives no input from population 2. Therefore, the macrostates and bifurcations of population 1 are identical to those explored in Luke et al. (2013), described above. However, we allow k21 ≠ 0. Our goal is to examine the consequences of the influence of population 1 on population 2. We call population 1 the “driver” and population 2 the “response” system. See Figure 1.

FIGURE 1
www.frontiersin.org

Figure 1. The driver-response configuration. k11 and k22 are the intra-population coupling strengths for populations 1 and 2, respectively, and k21 is the uni-directional coupling strength between the driver population (1) and the response population (2).

Writing the governing equation of population 2 as

z˙2=i(z21)22+(z2+1)22{Δ2+i[ηeff+k22H(z2)]}    (12)

with

ηeffη2+k21H(z1),    (13)

and comparing to Equation (6), we see that the behavior of population 2 is the same as that of a single population of theta neurons with an effective median excitability parameter ηeff. This effective parameter depends on the median excitability parameter intrinsic to population 2 η2, the inter-population coupling k21, and the state of the driver z1.

Note that ηeff depends linearly on both η2 and k21 and non-linearly on the driver's state z1 through H(z1). Additionally, ηeff may be time-dependent if population 1 exhibits a CPW state, since in that case z1 oscillates periodically. In the following, we will examine all these cases.

4. Results

We will examine the behavior of population 2 as various parameters are varied. We organize the presentation of our results by first considering the case in which the driver population exhibits an equilibrium state. Later, we consider the case in which the driver population exhibits periodic behavior.

We will mainly consider two configurations of the response system. The “excitatorily coupled” response system has k22 > 0, and the “inhibitorily coupled” response system has k22 < 0. Other parameters are as noted below.

The bifurcation diagrams that appear below in Figures 2, 3, 4B, 5B, 8C were obtained using XPPAUT (Ermentrout, 2002). Data for all other figures were generated using custom-designed code.

FIGURE 2
www.frontiersin.org

Figure 2. (A) A two-dimensional bifurcation diagram of the excitatorily-coupled response system. The heavy black lines are saddle-node (SN) bifurcation curves, and the solid dot denotes the parameters of the response system when decoupled from the driver. In the cases considered in the main text, the driver causes ηeff to vary along the horizontal dotted line. The parameters are: η1 = −0.2, Δ1 = 0.1, k11 = −2, and k22 = 9. (B) The one-dimensional bifurcation diagram showing the asymptotic values of y2 = Im(z2) vs. k21. Solid and dashed curves indicate stable and unstable equilibria, respectively, corresponding to partially synchronous spiking (PSS) and partially synchronous resting (PSR) states. The parameters are as in (A), with η2 = −10 and Δ2 = 0.5.

FIGURE 3
www.frontiersin.org

Figure 3. (A) The two-dimensional bifurcation diagram of the inhibitorily-coupled response system. The heavy black lines are saddle-node (SN) bifurcation curves, green is a homoclinic (HC) bifurcation curve, and red is an Andronov-Hopf (AH) bifurcation curve. The latter two curves emerge from a Bogdanov-Takens (BT) point. The solid dot denotes the parameters of the response system when decoupled from the driver. In the cases considered in the main text, the driver causes ηeff to vary along the horizontal dotted line. The parameters are: η1 = −0.2, Δ1 = 0.1, k11 = −2, and k22 = −9. (B) The one-dimensional bifurcation diagram showing the asymptotic value of x2 = Re(z2) vs. k21. Solid curves denote stable equilibria; dashed black curves are unstable equilibria. Green represents the maxima and minima of a collective periodic wave (CPW) limit cycle. The parameters are as in (A), with η2 = 5 and Δ2 = 0.5.

FIGURE 4
www.frontiersin.org

Figure 4. (A) The non-linear behavior of ηeff as a function of k11 for the excitatorily-coupled response system. ηeff is plotted horizontally to facilitate comparison with Figure 2A. The parameters are: η1 = −0.05, Δ1 = 0.05, η2 = −10, with the inter-population coupling fixed at k21 = 2.0. (B) The one-dimensional bifurcation diagram showing the asymptotic value of y2 = Im(z2) vs. k11. Solid and dashed curves indicate stable and unstable equilibria, respectively. The parameters are as in (A), with Δ2 = 0.5 and k22 = 9.

FIGURE 5
www.frontiersin.org

Figure 5. (A) The non-linear behavior of ηeff as a function of k11 for the inhibitorily-coupled response system. ηeff is plotted horizontally to facilitate comparison with Figure 3A. (B) The one-dimensional bifurcation diagram showing the asymptotic value of x2 = Im(z2) vs. k11. Solid and dashed black curves indicate stable and unstable equilibria, respectively, and green represents the maxima and minima of a CPW limit cycle state. The parameters are: η1 = −0.05, Δ1 = 0.05, η2 = 5, Δ2 = 0.5, and k22 = −9. The inter-population coupling is fixed at k21 = 3.5.

4.1. Driver on a Macroscopic Equilibrium

We begin by fixing the driving population's parameters at η1 = −0.2, Δ1 = 0.1, and k11 = −2, which corresponds to a PSR state. Thus, z1 remains fixed at a constant value. We examine the behavior of the two response system configurations as we vary the inter-population coupling parameter, k21. From Equation (13), ηeff varies linearly with respect to k21.

4.1.1. Excitatorily-coupled response system

We set the response system's internal coupling to k22 = 9, and show in Figure 2A the two-parameter bifurcation diagram of the response system with respect to Δ2 and ηeff. Two saddle-node bifurcation curves which meet at a cusp are seen. To the left of these curves, the response network exhibits a PSR state, and to the right, a PSS state. These states coexist inside the approximately triangular region.

We set the remaining parameters of the response system to η2 = −10 and Δ2 = 0.5. Thus, for k21 = 0, ηeff = η2, and the response system is situated at the solid black point marked in Figure 2A. As k21 increases from zero, ηeff increases linearly along the dotted line in Figure 2A, starting from the black point. In so doing, it traverses the SN bifurcation curves. Figure 2B shows how the imaginary part of the response's asymptotic macroscopic mean field [y2 = Im(z2)] changes with respect to k21, illustrating the coexistence of the stable PSR and PSS states, along with an unstable PSR state (uPSR).

The point marked “SN/NF” in Figure 2B indicates that as k21 increases, a saddle node bifurcation is encountered, corresponding to the left SN curve in Figure 2A. This creates a stable and an unstable PSS state. However, the unstable PSS state converts into an unstable PSR state at a value of k21 very slightly beyond the SN bifurcation. That is, the node corresponding to the unstable PSS state becomes a unstable PSR focus, a transition we called a Node-Focus (NF) transition in Luke et al. (2013). The distinction between these events is indistinguishable in the figure.

4.1.2. Inhibitorily-coupled response system

We performed a similar analysis for the case in which the response system's internal coupling is k22 = −9, i.e., inhibitory, and η2 = 5. The remaining parameters were unchanged. The results are shown in Figure 3. In this case, the two-dimensional bifurcation diagram of the response system with respect to Δ2 and ηeff (Figure 3A) shows a similar (but mirror-image) cusp of saddle-node curves. A new feature is the occurrence of a codimension-2 Bogdanov-Takens (BT) point on the left SN curve, and the emergence of homoclinic (HC; green) and Andronov-Hopf (AH; red) bifurcation curves from the BT point.

Figure 3B shows how the real part of the response's asymptotic macroscopic mean field [x2 = Re(z2)] changes with respect to k21. As before, ηeff increases linearly as k21 increases, starting from the black solid point in Figure 3A and moving toward the right, traversing the various bifurcation curves along the dotted line. Note the presence of the attracting limit cycle CPW state in Figure 3B, which emerges at the HC bifurcation and terminates at the AH bifurcation as k21 increases.

It is interesting to note that in both cases described above, the same bifurcation structure would be encountered if, instead of varying k21 with a fixed value η2, we varied η2 with a fixed value of k21. While this is obvious from Equation (13) since H(z_1) is constant in these cases, this leads to the non-obvious conclusion that by modifying either the inter-population coupling or the intrinsic median excitability of the response population—two rather different system characteristics—one obtains identical transitions in the response network.

4.1.3. Variation of the driver's macroscopic equilibrium

In the cases we considered previously, ηeff changed linearly with respect to the inter-population coupling k21. We now turn our attention to the effects incurred by altering the value of the driver influence function H(z1) in Equation (13). We do this by varying the driver's internal coupling strength k11, thus causing the driver's asymptotic macroscopic mean field z1 to change. This manipulation has the effect of changing ηeff non-linearly with respect to k11.

For simplicity, we only consider a range of k11 such that the driver always remains on a macroscopic equilibrium state, and we fix the inter-population coupling at k21 = 2.

We begin with the case of the excitatorily-coupled response system considered above, with η2 = −10, Δ2 = 0.5, and k22 = 9, and choose the remaining driver parameters to be η1 = −0.05 and Δ1 = 0.05. Figure 4A shows the non-linear behavior of ηeff as k11 is varied. Even though we are considering k11 to be the independent parameter, we plot ηeff horizontally so that it may be easily compared to Figure 2A; recall that this shows the two-dimensional bifurcation diagram of the response system. Now, as k11 changes, ηeff moves back and forth along the dotted line non-linearly. In particular, Figure 4A shows that for very negative values of k11, ηeff is near −5, which corresponds to a point in Figure 2A to the right of the SN curves. As k11 increases, ηeff decreases to approximately −10, thus crossing both SN curves in Figure 2A from right to left in the process. ηeff subsequently increases, and goes back across the SN curves from left to right. Note that Figure 4A includes vertical lines marking the position of the SN bifurcations (i.e., the values of ηeff at which the horizontal line at Δ2 = 0.5 in Figure 2A crosses the SN curves).

Figure 4B shows the behavior of the asymptotic state of the response system [y2 = Im(z2)] as a function of k11. This shows that as k11 increases, the response system passes through two separate regions of bistability, corresponding to the two traversals of the triangular bistable region in Figure 2A. Thus, Figure 4B is qualitatively similar to two copies of Figure 2B, with the structure for k11 < 0 reversed. Note that the two regions are not symmetrical. This is due to the non-symmetric behavior of ηeff as k11 changes.

Next, we examine how the same manipulation of the driver system affects the inhibitorily-coupled response system. The parameters are as above, but with η2 = 5 and k22 = −9. Figure 5A shows how ηeff changes as k11 is varied, again plotted with ηeff on the horizontal axis for ease of comparison with Figure 3A. Note the vertical lines in Figure 5A marking the SN, HC, and AH bifurcations.

The one-dimensional bifurcation diagram depicting the asymptotic state of the response system as a function of k11 is shown in Figure 5B. A situation similar to the previous case is observed. Two distorted versions of the structure of Figure 3B, with the features for k11 < 0 being reversed, are seen. Again, this is due to the non-linear and asymmetric behavior of ηeff as it traverses the bifurcations in Figure 3A twice: first right to left, and then left to right, as k11 is increased. Note also the presence of an attracting limit cycle CPW state in intervals of both positive and negative k11.

4.2. Driver on a Macroscopic Limit Cycle

We now focus on the behavior of the response population when the driver is on a CPW state, which is a limit cycle of the driver's macroscopic mean field (or order parameter). Throughout this section, we fix the driver parameters at η1 = 10.75, k11 = −9, and Δ1 = 0.5, which results in a CPW driver state for which H(z1) oscillates periodically in time. In particular, we have H(z1) > 0 for all time. Thus, according to Equation (13), ηeff also oscillates periodically for k21 ≠ 0, and both the centroid and the amplitude of the ηeff oscillation increase as k21 increases.

We show below that in this configuration, the response population can exhibit periodic, multistable, chaotic, and/or quasiperiodic behavior, depending on the response system's parameters and the interpopulation coupling strength k21.

4.2.1. Periodic behavior in the response system

We begin by considering the excitatorily coupled response system, with Δ2 = 0.5 and k22 = 9, but with η2 = −20. When decoupled from the driver, this places the response system at a point well to the left in the parameter space of Figure 2A. Thus, the response system in isolation asymptotes to a PSR state. As k21 is increased from zero to eight, ηeff oscillates back and forth along the horizontal line in Figure 2A at Δ2 = 0.5, but always stays to the left of the SN curves shown in that figure. Thus, the driver simply pushes the response system's PSR state back and forth, avoiding any bifurcations. The result is simple periodic behavior in the driven response system. Figure 6A shows a plot of the maximum and minimum of x2 = Re(z2) vs. k21. As k21 increases, the amplitude of this simple periodic behavior increases. We observe that the frequency of the response system's oscillation is the same as that of the driver throughout this range of interpopulation coupling.

FIGURE 6
www.frontiersin.org

Figure 6. (A) Simple periodic behavior in the response system driven by a CPW state of the driver as a function of the inter-population coupling strength k21. The curves are local maxima and minima of x2 = Re(z2). The driver parameters are η1 = 10.75, Δ1 = 0.5, and k11 = −9, and the response parameters are η2 = −20, Δ2 = 0.5, and k22 = 9. (B) Slightly more complicated periodic behavior obtained at the same parameters, except with η2 = −5. The curves are local maxima and minima of y2 = Im(z2).

We now change the response system such that η2 = −5, and leave all other parameters the same as above. This change places the response system at a point to the right of the SN curves in Figure 2A, and for these parameters, the uncoupled response system asymptotes to a PSS state. Once again, as k21 increases, ηeff oscillates back and forth along the Δ2 = 0.5 line in Figure 2A, but this time it does so always staying to the right of the SN curves.

The result is multi-frequency periodic behavior in the response system that is more complicated than in the previous example. Figure 6B shows a plot of the local minima and maxima of y2 = Im(z2) vs. k21. Figure 7 shows y2 vs. x2 plots of the periodic orbits at k21 = 6 (upper panels) and k21 = 10 (lower panels). As k21 increases from zero, a periodic orbit with winding number two emerges (similar to that shown in Figure 7A) and grows in amplitude, peaking near k21 ≈ 2.5. The amplitude subsequently decreases to a minimum near k21 ≈ 7.2, and then slowly increases again. Note that the four curves in Figure 6B for k21 ∈ [0, 7.2] correspond to two pairs of alternating local maxima and minima in the time series of y2, as shown in Figure 7B.

FIGURE 7
www.frontiersin.org

Figure 7. The response system's behavior at parameters corresponding to Figure 6B at k21 = 6 (A,B) and k21 = 10 (C,D), with z2 = x2 + iy2.

Interestingly, near k21 ≈ 7.2, an additional loop appears in the orbit, as shown in Figure 7C. This is reflected in the additional inner curves in Figure 6B that appear for k21 ⪆ 7.2, and the two additional local maxima and minima in the time series of y2 in Figure 7D.

4.2.2. Multistability in the response system

Continuing with the excitatorily coupled response system (with k22 = 9 > 0), we set η2 = −10 and leave all other parameters unchanged. In this case the uncoupled response system is at a point just to the left of the left SN curve in Figure 2A, and as k21 increases, ηeff again sweeps back and forth along the horizontal line at Δ1 = 0.5. However, now this sweeping cuts across both SN curves. Thus, the response system sweeps back and forth across the approximately triangular multistable region bounded by the SN curves.

Figure 8A shows the maxima and minima of y2 vs. k21 for this case. The first feature to emerge as k21 increases from zero is a simple periodic orbit whose amplitude increases, similar to the example in Figure 6A. At k21 ≈ 0.5, a new and separate coexisting limit cycle appears, as indicated by the upper curves that emerge in Figure 8A. Figure 8B shows the y2 vs. x2 plots of these two limit cycles at k21 = 1.5, where the larger orbit corresponds to the upper two curves in Figure 8A. In this bistable region, the macroscopic dynamics of the response system approaches one or the other of these periodic states, depending on the initial conditions.

FIGURE 8
www.frontiersin.org

Figure 8. Multistability in the response system driven by a CPW state of the driver. (A) Local maxima and minima of y2 = Re(z2) vs. the inter-population coupling k21. (B) y2 vs. x2 plots showing two coexisting limit cycles of the response system at k21 = 1.5 (dotted vertical line in A). (C) The solid and dashed black curves show the asymptotic states of the response for fixed values of ηeff, with k21 = 1.5. Green curves are coexisting limit cycles of the response system when coupled to the driver. Parameters are: η1 = 10.75, Δ1 = 0.5, k11 = −9; η2 = −10, Δ2 = 0.5, k22 = 9.

Figure 8C shows, in black, the asymptotic states of y2 vs. ηeff for fixed values of ηeff, with k21 = 1.5. These curves show that for a large interval of ηeff, a stable PSR coexists with a stable PSS and an unstable PSR state for the frozen (i.e., ηeff fixed) system. With the driver on the CPW state, ηeff sweeps from approximately −9.1 to −7.6 and back again–a range which is well within the bistable region. Superimposed in green in Figure 8C are projections of the two coexisting limit cycles onto this space, showing that the lower limit cycle is a simple periodic perturbation of the response system's underlying PSR state, and the upper limit cycle is a periodic perturbation of the underlying PSS state.

4.2.3. Chaos in the response system

We now switch to the inhibitorily coupled response system, with parameters η2 = 5, Δ2 = 0.5, and k22 = −9. The parameter space of this system is shown in Figure 3A, and the uncoupled response system resides at the solid black dot in that figure, to the left of all the bifurcations. As the interpopulation coupling strength k21 increases, ηeff sweeps across the same horizontal line at Δ2 = 0.5 with increasing amplitude and centroid, initially crossing just the left SN bifurcation curve. At k21 ≈ 5.2, ηeff begins sweeping across the homoclinic and the Andronov-Hopf bifurcation curves. Eventually, for sufficiently large k21, ηeff sweeps across all four bifurcation curves (SN, AH, HC, and SN).

Figure 9A shows the local maxima and minima of x2 = Re(z2) vs. k21. We initially see the emergence of a simple periodic orbit that grows slowly in amplitude. However, at k21 ≈ 5.2, chaos suddenly emerges through a crisis. Figure 9B shows a magnification of this region, with a plot of the two largest Lyapunov exponents. We see that there are significant intervals of k21 for which there is a positive Lyapunov exponent, indicating the presence of macroscopic chaos.

FIGURE 9
www.frontiersin.org

Figure 9. Emergence of macroscopic chaos in the response system driven by a CPW state of the driver. (A) Local minima and maxima of x2 = Re(z2) vs. the inter-population coupling k21. (B) Magnification of the chaotic region (top), with a plot of the largest two Lyapunov exponents. Parameters are: Δ1 = 0.5, k11 = −9, η1 = 10.75; Δ2 = 0.5, k22 = −9, η2 = 5.

As k21 increases, the first chaotic band, beginning at k21 ≈ 5.28, coexists with the simple periodic loop that was present for smaller k21 (this coexistence is not apparent in the figure). Outside of this band, there is a window dominated by periodic behavior of rather high period. A second chaotic band appears at approximately k21 = 5.48. This second band terminates at approximately k21 = 5.65, after which a series of reverse period-doubling cascades are seen.

The y2 vs. x2 plot of the chaotic attractor present at k21 = 5.296, for which the largest Lyapunov exponent is approximately 0.2118, is shown in Figure 10A.

FIGURE 10
www.frontiersin.org

Figure 10. (A) Chaotic (y2 vs. x2) and (B) Quasiperiodic (y2 vs. x1 vs. y1) attractors in the response system driven by a CPW state of the driver. Parameters are: k11 = k22 = −9, with (A) η2 = 5, Δ1 = Δ2 = 0.5 and k21 = 5.296, and (B) η2 = 10.75, Δ1 = 0.5, Δ2 = 0.3, and k21 = 0.1.

4.2.4. Quasiperiodicity in the response system

Finally, we consider the case in which the response system exhibits a CPW state when uncoupled from the driver, and ask what happens when this is driven by another CPW state in the driver. We use the same drive system parameters as above, and set the response system's parameters to be the same except for Δ2 = 0.3. As the inter-population coupling strength k21 is increased, various phase-locked and quasiperiodic states are seen. An example of quasiperiodic behavior in the response system for k21 = 0.1 is shown in Figure 10B.

5. Discussion

In this work, we have taken the first step toward designing a mathematically tractable modular network-of-networks representation of neuronal systems. Our approach is based on dynamical analysis techniques that enable a complete description of the emergent macroscopic behavior of large, heterogeneous discrete networks of globally-coupled phase oscillators. Building on previous results (Luke et al., 2013) in which we used these techniques to show that the collective dynamics of a single such population of theta neurons is relatively simple (exhibiting just equilibria and limit cycle states), we constructed the next simplest hierarchical structure: a driver-response configuration of theta neuron populations. Our results show that even in this simplest of configurations, the response system (and hence, the network as a whole) can exhibit a full range of dynamical behaviors and surprising complexity. A notable strength of our work is that despite the complexity that emerges from this arrangement, the behavior can be understood and explained in terms of what is known about a single population's dynamics and bifurcation structure.

With the driving system on a fixed equilibrium, we showed that the response system is equivalent to a single population with a simple shift in one parameter. Specifically, this parameter is the median of the distribution of excitability parameters in the response system, which indicates whether the response population is dominated by excitable or intrinsically-spiking neurons. Although this arrangement does not introduce any new dynamical features, we showed that the response system can nevertheless still exhibit an interesting bifurcation structure involving macroscopic equilibria, limit cycles, and multistability as the strength of the inter-population coupling varies. More interestingly, we found that the inter-population coupling strength is effectively equivalent to the response system's median intra-population excitability. By this we mean that changes in either of these rather different network parameters lead to identical bifurcation scenarios. This surprising result follows from the drive-response network configuration in particular.

The first level of additional complication arose when modestly altering an internal parameter of the drive system. This effectively led to a non-linear change in the response system's median excitability, causing a dramatic change in the response's bifurcation structure. Such bifurcation structures might be difficult to understand if encountered blindly, as might be the case when studying the dynamics of a network without knowledge of its internal structure. Experimental studies of neuronal networks often take a similar “black box” approach out of necessity, since detailed knowledge of connectivity (i.e., the “connectome”) is rarely available. In our case, however, we showed that knowledge of the non-linearity, along with knowledge of the bifurcation structure of a single network, leads to a natural explanation of the additional features that arise due to the network-of-networks structure. In our particular case studies, we observed multiple distorted and reversed copies of the bifurcation structure that is associated with a single population of theta neurons. We therefore speculate that in “black box” investigations, the observation of such repeated and/or distorted bifurcation structures might be indicative of driver-response-type connectivity in the network of study.

Finally, we investigated the consequences of placing the driver system on a collective rhythmic state (i.e., a macroscopic periodic orbit). Our results were consistent with previous results that studied non-autonomous phase oscillator (So and Barreto, 2011) and theta neuron systems (So et al., 2014). In those investigations, it was shown that networks of oscillators subjected to a sinusoidal variation of a network parameter led to complicated dynamics including quasiperiodicity and macroscopic chaos. Here, our driver-response arrangement of two separate interacting populations of theta neurons leads to an overall autonomous system, but with the response system being subjected to a periodic driving signal from the driver. Such arrangements might be found in real neuronal systems at the early stages of sensory input processing. For example, the lateral geniculate nucleus may be driven by a periodic visual signal delivered to the retina. Another candidate might be the trisynaptic circuit of the dentate gyrus and the CA3 and CA1 regions of the hippocampus (Kandel et al., 2000). More generally, the information-processing capabilities of the brain are thought to be regulated by collective rhythms, notably theta and gamma oscillations, which arise in various areas and periodically drive other areas (Buzsáki, 2006).

Our results may also have implications for populations of bursting neurons (So et al., 2014). Neuronal bursting in individual neurons is commonly understood to arise as the result of the interplay between a slowly oscillating neuronal parameter (or “slow variable”) and the neuron's fast spiking dynamics. Bursting arises if the slow parameter sweeps back and forth across bifurcations, and (Rinzel and Ermentrout, 1989) classified bursters as square, parabolic, or elliptic based on the bifurcations encountered in this process. It has also been demonstrated that slowly oscillating intra- and extra-cellular ion concentrations can lead to wide range of neuronal bursting behaviors (Cressman et al., 2009, 2011; Barreto and Cressman, 2011).

Finally, we note that our explorations in this work were limited to cases in which the driver population's parameters were either fixed or were varied only modestly. In the latter case, we changed the driver's median excitability parameter only to the extent that its collective equilibrium state was displaced but not altered. Significantly greater complexity in the response's dynamics would arise if the collective state of the driver were pushed across its own bifurcations, possibly resulting in topological changes and hysteretic effects in the driver's macroscopic state. As discussed above, such complexity would be difficult to understand if encountered in a “black box”-type investigation. Nevertheless, if it is known that the network of interest has a driver-response structure, it may be possible to comprehend the origin of such complexity in the manner that we have outlined here.

This study constitutes an initial attempt at building a mathematically tractable model to understand the collective behavior of a hierarchical “network-of-networks” arrangement of model neurons. In future work we plan to consider networks of networks that include feedback connections and additional populations in an effort to understand the emergence of macroscopic dynamical complexity in more realistic networks.

Author Contributions

Tanushree B. Luke, Ernest Barreto, and Paul So conceived and designed the investigation, analyzed the data, and wrote the paper. Tanushree B. Luke and Paul So performed the numerical computations.

Conflict of Interest Statement

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.

Acknowledgment

Publication of this article was funded by the George Mason University Libraries Open Access Publishing Fund.

Footnotes

1. ^an=2π/ππ(1cos(x))n=n!/(2n1)!!

References

Bacci, A., Huguenard, J. R., and Prince, D. A. (2003). Functional autaptic neurotransmission in fast-spiking interneurons: a novel form of feedback inhibition in the neocortex. J. Neurosci. 23, 859–866. Available online at: http://www.jneurosci.org/content/23/3/859.abstract

Barreto, E., and Cressman, J. (2011). Ion concentration dynamics as a mechanism for neuronal bursting. J. Biol. Phys. 37, 361–373. doi: 10.1007/s10867-010-9212-6

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Bekkers, J. M. (2003). Synaptic transmission: functional autapses in the cortex. Curr. Biol. 13, R433–R435. doi: 10.1016/S0960-9822(03)00363-4

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Bullmore, E., and Sporns, O. (2009). Complex brain network: graph theoretical analysis of structural and functional systems. Nat. Rev. Neurosci. 10, 186–198. doi: 10.1038/nrn2575

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Buzsáki, G. (2006). Rhythms of The Brain, 1 Edn. New York, NY: Oxford University Press. doi: 10.1093/acprof:oso/9780195301069.001.0001

CrossRef Full Text

Cressman, J., Ullah, G., Ziburkus, J., Schiff, S., and Barreto, E. (2009). The influence of sodium and potassium dynamics on excitability, seizures, and the stability of persistent states: 1. Single neuron dynamics. J. Comput. Neurosci. 26, 159–170. doi: 10.1007/s10827-008-0132-4

CrossRef Full Text | Google Scholar

Cressman, J., Ullah, G., Ziburkus, J., Schiff, S., and Barreto, E. (2011). Erratum to: the influence of sodium and potassium dynamics on excitability, seizures, and the stability of persistent states: I. single neuron dynamics. J. Comput. Neurosci. 30, 781. doi: 10.1007/s10827-011-0333-0

CrossRef Full Text | Google Scholar

Ermentrout, G. (1996). Type I membranes, phase resetting curves and synchrony. Neural Comput. 8, 979–1001. doi: 10.1162/neco.1996.8.5.979

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Ermentrout, G. (2002). Simulating, Analyzing, and Animating Dynamical Systems (A Guide to XPPAUT for Researchers and Students. Philadelphia, PA: SIAM. doi: 10.1137/1.9780898718195

CrossRef Full Text | Google Scholar

Ermentrout, G., and Kopell, N. (1986). Parabolic bursting in an excitable system coupled with a slow oscillation. SIAM J. Appl. Math. 46, 233–253. doi: 10.1137/0146017

CrossRef Full Text | Google Scholar

Harris, K. (2005). Neural signatures of cell assembly organization. Nat. Rev. Neurosci. 6, 399–407. doi: 10.1038/nrn1669

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Hebb, D. (1949). The Organization of Behavior: A Neuropsychological Theory, 1 Edn. New York, NY: Wiley.

Google Scholar

Hodgkin, A. (1948). The local electric charges associted with repetitive action in non-medulated axons. J. Physiol. 107, 165–181.

Pubmed Abstract | Pubmed Full Text | Google Scholar

Izhikevich, E. (2007). Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting. Cambridge, MA: MIT Press.

Google Scholar

Kandel, E., Schwartz, J., and Jessell, T. (2000). Principles of Neural Science, 4 Edn. New York, NY: McGraw-Hill.

Laing, C. (2014). Derivation of a neural field model from a network of theta neurons. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 90:010901. doi: 10.1103/PhysRevE.90.010901

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Luke, T., Barreto, E., and So, P. (2013). Complete classification of the macroscopic behavior of a heterogeneous network of theta neurons. Neural Comput. 25, 3207–3234. doi: 10.1162/NECO_a_00525

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Marvel, S. A., Mirollo, R. E., and Strogatz, S. H. (2009). Identical phase oscillators with global sinusoidal coupling evolve by mbius group action. Chaos 19:043104. doi: 10.1063/1.3247089

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Meunier, D., Lambiotte, R., and Bullmore, E. (2010). Modular and hierarchically modular organization of brain networks. Front. Neurosci. 4:200. doi: 10.3389/fnins.2010.00200

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Nowak, L., Azouz, R., Sanchez-Vives, M., Gray, C., and McCormick, D. (2003). Electrophysiological classes of cat primary visual cortical neurons in vivo as revealed by quantitative analyses. J. Neurophysiol. 89, 1541–1566. doi: 10.1152/jn.00580.2002

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Ott, E., and Antonsen, T. (2008). Low dimensional behavior of large systems of globally coupled oscillators. Chaos 18, 037113. doi: 10.1063/1.2930766

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Ott, E., and Antonsen, T. (2009). Long time evolution of phase oscillator systems. Chaos 19, 023117. doi: 10.1063/1.3136851

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Ott, E., Hunt, B. R., and Antonsen, T. M. (2011). Comment on long time evolution of phase oscillator systems [Chaos 19, 023117 (2009)]. Chaos 21, 025112. doi: 10.1063/1.3574931

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Pazó, D., and Montbrió, E. (2014). Low-dimensional dynamics of populations of pulse-coupled oscillators. Phys. Rev. X 4:011009. doi: 10.1103/PhysRevX.4.011009

CrossRef Full Text | Google Scholar

Pikovsky, A., and Rosenblum, M. (2011). Dynamics of heterogeneous oscillator ensembles in terms of collective variables. Phys. D Nonlinear Phenom. 240, 872–881. doi: 10.1016/j.physd.2011.01.002

CrossRef Full Text | Google Scholar

Rinzel, J., and Ermentrout, G. (1989). “Analysis of neural excitability and oscillations,” in Methods in Neuronal Modeling, eds C. Koch and I. Segev (Cambridge, MA: The MIT Press), 251–291.

Sherrington, C. (1906). The Integrative Action of The Nervous System, 1 Edn. New York, NY: Scribners.

Google Scholar

So, P., and Barreto, E. (2011). Generating macroscopic chaos in a network of globally coupled phase oscillators. Chaos 21, 033127. doi: 10.1063/1.3638441

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

So, P., Luke, T., and Barreto, E. (2014). Networks of theta neurons with time-varying excitability: macroscopic chaos, multistability, and final-state uncertainty. Phys. D 267, 16–26. doi: 10.1016/j.physd.2013.04.009

CrossRef Full Text | Google Scholar

Tateno, T., Harsch, A., and Robinson, H. (2004). Threshold firing frequency-current relationships of neurons in rat somatosensory cortex: type 1 and type 2 dynamics. J. Neurophysiol. 92, 2283–2294. doi: 10.1152/jn.00109.2004

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Zhou, C., Zemanova, L., Zamora, G., Hilgetag, C., and Kurths, J. (2006). Hierarchical organization unveiled by functional connectivity in complex brain networks. Phys. Rev. Lett. 97:238103. doi: 10.1103/PhysRevLett.97.238103

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Keywords: theta neuron, type-I neuron, hierarchical network, neural field, macroscopic behavior, coherence, synchrony, chaos

Citation: Luke TB, Barreto E and So P (2014) Macroscopic complexity from an autonomous network of networks of theta neurons. Front. Comput. Neurosci. 8:145. doi: 10.3389/fncom.2014.00145

Received: 31 August 2014; Accepted: 27 October 2014;
Published online: 18 November 2014.

Edited by:

Tobias Alecio Mattei, Ohio State University, USA

Reviewed by:

Jianbo Gao, Wright State University, USA
Yasuhiro Tsubo, Ritsumeikan University, Japan

Copyright © 2014 Luke, Barreto and So. 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) or licensor 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: Paul So, George Mason University, Mail Stop 2A1, Fairfax, VA 22030, USA e-mail: paso@gmu.edu