# Macroscopic complexity from an autonomous network of networks of theta neurons

- 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

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:

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

where *P _{n}*(θ) =

*a*(1 − cosθ)

_{n}^{n},

*n*∈ ℕ, and

*a*is a normalization constant

_{n}^{1}such that

The parameter *n* defines the sharpness of the pulse-like synapse in that *P _{n}*(θ) 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,

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

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 |$\tilde{{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

where

and

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 *H _{n}*(

*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:

where θ_{1,j} and θ_{2,j} denote the *j*th 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 $\tilde{{z}}$_{1}(*t*), $\tilde{{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:

with *H*_{n}(*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 *N*_{1}, *N*_{2} → ∞.

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 *H _{n}* to ease notation.

### 3.2. The Driver-Response System

To put our network in the driver-response form, we set *k*_{12} = 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 *k*_{21} ≠ 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. The driver-response configuration**. *k*_{11} and *k*_{22} are the intra-population coupling strengths for populations 1 and 2, respectively, and *k*_{21} is the uni-directional coupling strength between the driver population (1) and the response population (2).

Writing the governing equation of population 2 as

with

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 *k*_{21}, and the state of the driver *z*_{1}.

Note that η_{eff} depends linearly on both η_{2} and *k*_{21} and non-linearly on the driver's state *z*_{1} through *H*(*z*_{1}). Additionally, η_{eff} may be time-dependent if population 1 exhibits a CPW state, since in that case *z*_{1} 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 *k*_{22} > 0, and the “inhibitorily coupled” response system has *k*_{22} < 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. (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, *k*_{11} = −2, and *k*_{22} = 9. **(B)** The one-dimensional bifurcation diagram showing the asymptotic values of *y*_{2} = Im(*z*_{2}) vs. *k*_{21}. 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. (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, *k*_{11} = −2, and *k*_{22} = −9. **(B)** The one-dimensional bifurcation diagram showing the asymptotic value of *x*_{2} = Re(*z*_{2}) vs. *k*_{21}. 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. (A)** The non-linear behavior of η_{eff} as a function of *k*_{11} 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 *k*_{21} = 2.0. **(B)** The one-dimensional bifurcation diagram showing the asymptotic value of *y*_{2} = Im(*z*_{2}) vs. *k*_{11}. Solid and dashed curves indicate stable and unstable equilibria, respectively. The parameters are as in **(A)**, with Δ_{2} = 0.5 and *k*_{22} = 9.

**Figure 5. (A)** The non-linear behavior of η_{eff} as a function of *k*_{11} 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 *x*_{2} = Im(*z*_{2}) vs. *k*_{11}. 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 *k*_{22} = −9. The inter-population coupling is fixed at *k*_{21} = 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 *k*_{11} = −2, which corresponds to a PSR state. Thus, *z*_{1} remains fixed at a constant value. We examine the behavior of the two response system configurations as we vary the inter-population coupling parameter, *k*_{21}. From Equation (13), η_{eff} varies linearly with respect to *k*_{21}.

#### 4.1.1. Excitatorily-coupled response system

We set the response system's internal coupling to *k*_{22} = 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 *k*_{21} = 0, η_{eff} = η_{2}, and the response system is situated at the solid black point marked in Figure 2A. As *k*_{21} 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 [*y*_{2} = Im(*z*_{2})] changes with respect to *k*_{21}, 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 *k*_{21} 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 *k*_{21} 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 *k*_{22} = −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 [*x*_{2} = Re(*z*_{2})] changes with respect to *k*_{21}. As before, η_{eff} increases linearly as *k*_{21} 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 *k*_{21} increases.

It is interesting to note that in both cases described above, the same bifurcation structure would be encountered if, instead of varying *k*_{21} with a fixed value η_{2}, we varied η_{2} with a fixed value of *k*_{21}. 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 *k*_{21}. We now turn our attention to the effects incurred by altering the value of the driver influence function *H*(*z*_{1}) in Equation (13). We do this by varying the driver's internal coupling strength *k*_{11}, thus causing the driver's asymptotic macroscopic mean field *z*_{1} to change. This manipulation has the effect of changing η_{eff} *non-linearly* with respect to *k*_{11}.

For simplicity, we only consider a range of *k*_{11} such that the driver always remains on a macroscopic equilibrium state, and we fix the inter-population coupling at *k*_{21} = 2.

We begin with the case of the excitatorily-coupled response system considered above, with η_{2} = −10, Δ_{2} = 0.5, and *k*_{22} = 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 *k*_{11} is varied. Even though we are considering *k*_{11} 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 *k*_{11} changes, η_{eff} moves back and forth along the dotted line non-linearly. In particular, Figure 4A shows that for very negative values of *k*_{11}, η_{eff} is near −5, which corresponds to a point in Figure 2A to the right of the SN curves. As *k*_{11} 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 [*y*_{2} = Im(*z*_{2})] as a function of *k*_{11}. This shows that as *k*_{11} 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 *k*_{11} < 0 reversed. Note that the two regions are not symmetrical. This is due to the non-symmetric behavior of η_{eff} as *k*_{11} 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 *k*_{22} = −9. Figure 5A shows how η_{eff} changes as *k*_{11} 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 *k*_{11} 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 *k*_{11} < 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 *k*_{11} is increased. Note also the presence of an attracting limit cycle CPW state in intervals of both positive and negative *k*_{11}.

### 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, *k*_{11} = −9, and Δ_{1} = 0.5, which results in a CPW driver state for which *H*(*z*_{1}) oscillates periodically in time. In particular, we have *H*(*z*_{1}) > 0 for all time. Thus, according to Equation (13), η_{eff} also oscillates periodically for *k*_{21} ≠ 0, and both the centroid and the amplitude of the η_{eff} oscillation increase as *k*_{21} 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 *k*_{21}.

#### 4.2.1. Periodic behavior in the response system

We begin by considering the excitatorily coupled response system, with Δ_{2} = 0.5 and *k*_{22} = 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 *k*_{21} 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 *x*_{2} = Re(*z*_{2}) vs. *k*_{21}. As *k*_{21} 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. (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 *k*_{21}. The curves are local maxima and minima of *x*_{2} = Re(*z*_{2}). The driver parameters are η_{1} = 10.75, Δ_{1} = 0.5, and *k*_{11} = −9, and the response parameters are η_{2} = −20, Δ_{2} = 0.5, and *k*_{22} = 9. **(B)** Slightly more complicated periodic behavior obtained at the same parameters, except with η_{2} = −5. The curves are local maxima and minima of *y*_{2} = Im(*z*_{2}).

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 *k*_{21} 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 *y*_{2} = Im(*z*_{2}) vs. *k*_{21}. Figure 7 shows *y*_{2} vs. *x*_{2} plots of the periodic orbits at *k*_{21} = 6 (upper panels) and *k*_{21} = 10 (lower panels). As *k*_{21} increases from zero, a periodic orbit with winding number two emerges (similar to that shown in Figure 7A) and grows in amplitude, peaking near *k*_{21} ≈ 2.5. The amplitude subsequently decreases to a minimum near *k*_{21} ≈ 7.2, and then slowly increases again. Note that the four curves in Figure 6B for *k*_{21} ∈ [0, 7.2] correspond to two pairs of alternating local maxima and minima in the time series of *y*_{2}, as shown in Figure 7B.

**Figure 7. The response system's behavior at parameters corresponding to Figure 6B at k_{21} = 6 (A,B) and k_{21} = 10 (C,D), with z_{2} = x_{2} + iy_{2}**.

Interestingly, near *k*_{21} ≈ 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 *k*_{21} ⪆ 7.2, and the two additional local maxima and minima in the time series of *y*_{2} in Figure 7D.

#### 4.2.2. Multistability in the response system

Continuing with the excitatorily coupled response system (with *k*_{22} = 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 *k*_{21} 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 *y*_{2} vs. *k*_{21} for this case. The first feature to emerge as *k*_{21} increases from zero is a simple periodic orbit whose amplitude increases, similar to the example in Figure 6A. At *k*_{21} ≈ 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 *y*_{2} vs. *x*_{2} plots of these two limit cycles at *k*_{21} = 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. Multistability in the response system driven by a CPW state of the driver. (A)** Local maxima and minima of *y*_{2} = Re(*z*_{2}) vs. the inter-population coupling *k*_{21}. **(B)** *y*_{2} vs. *x*_{2} plots showing two coexisting limit cycles of the response system at *k*_{21} = 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 *k*_{21} = 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, *k*_{11} = −9; η_{2} = −10, Δ_{2} = 0.5, *k*_{22} = 9.

Figure 8C shows, in black, the asymptotic states of *y*_{2} vs. η_{eff} for *fixed* values of η_{eff}, with *k*_{21} = 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 *k*_{22} = −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 *k*_{21} 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 *k*_{21} ≈ 5.2, η_{eff} begins sweeping across the homoclinic and the Andronov-Hopf bifurcation curves. Eventually, for sufficiently large *k*_{21}, η_{eff} sweeps across all four bifurcation curves (SN, AH, HC, and SN).

Figure 9A shows the local maxima and minima of *x*_{2} = Re(*z*_{2}) vs. *k*_{21}. We initially see the emergence of a simple periodic orbit that grows slowly in amplitude. However, at *k*_{21} ≈ 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 *k*_{21} for which there is a positive Lyapunov exponent, indicating the presence of macroscopic chaos.

**Figure 9. Emergence of macroscopic chaos in the response system driven by a CPW state of the driver. (A)** Local minima and maxima of *x*_{2} = Re(*z*_{2}) vs. the inter-population coupling *k*_{21}. **(B)** Magnification of the chaotic region (top), with a plot of the largest two Lyapunov exponents. Parameters are: Δ_{1} = 0.5, *k*_{11} = −9, η_{1} = 10.75; Δ_{2} = 0.5, *k*_{22} = −9, η_{2} = 5.

As *k*_{21} increases, the first chaotic band, beginning at *k*_{21} ≈ 5.28, coexists with the simple periodic loop that was present for smaller *k*_{21} (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 *k*_{21} = 5.48. This second band terminates at approximately *k*_{21} = 5.65, after which a series of reverse period-doubling cascades are seen.

The *y*_{2} vs. *x*_{2} plot of the chaotic attractor present at *k*_{21} = 5.296, for which the largest Lyapunov exponent is approximately 0.2118, is shown in Figure 10A.

**Figure 10. (A)** Chaotic (*y*_{2} vs. *x*_{2}) and **(B)** Quasiperiodic (*y*_{2} vs. *x*_{1} vs. *y*_{1}) attractors in the response system driven by a CPW state of the driver. Parameters are: *k*_{11} = k_{22} = −9, with **(A)** η_{2} = 5, Δ_{1} = Δ_{2} = 0.5 and *k*_{21} = 5.296, and **(B)** η_{2} = 10.75, Δ_{1} = 0.5, Δ_{2} = 0.3, and *k*_{21} = 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 *k*_{21} is increased, various phase-locked and quasiperiodic states are seen. An example of quasiperiodic behavior in the response system for *k*_{21} = 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. ^${{a}}_{{n}}{=}{2}{\pi}{/}{\displaystyle {{\int}}_{{-}{\pi}}^{{\pi}}{(}{1}{-}{\mathrm{cos}}{(}}{x}{)}{{)}}^{{n}}{=}{n}{!}{/}{(}{2}{n}{-}{1}{)}{!}{!}$

## 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

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

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

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

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

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.

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

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

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

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

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.

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

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, USACopyright © 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