Macroscopic complexity from an autonomous network of networks of theta neurons

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.


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 informationprocessing 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 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.

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."

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 n (1 − cos θ ) n , n ∈ N, and a n is a normalization constant 1 such that 2π 0 P n (θ)dθ = 2π.
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, 1 a n = 2π/ π −π (1 − cos (x)) n = n!/(2n − 1)!! where η 0 is the center of the distribution, and , the half-width at half-maximum, describes the degree of heterogeneity in the population.

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 |z(t)| ∈ [0, 1] quantifies the degree of synchronization present at time t.
In Luke et al. (2013), we used the Ott-Antonsen method Antonsen, 2008, 2009;Ott et al., 2011) to derive a lowdimensional 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 iṡ 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.

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.

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 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:ż 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.

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.
Writing the governing equation of population 2 aṡ 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 nonlinearly on the driver's state z 1 through H(z 1 ). Additionally, η eff 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).

Frontiers in Computational Neuroscience
www.frontiersin.org November 2014 | Volume 8 | Article 145 | 3 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.

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.

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

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

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.

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 = Frontiers in Computational Neuroscience www.frontiersin.org November 2014 | Volume 8 | Article 145 | 5 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 saddlenode 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.

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 twodimensional 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 .

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 .

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 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 ).
as that of the driver throughout this range of interpopulation coupling.
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.
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.

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

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

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.

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., 2009Barreto 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.