Macroscopic complexity from an autonomous network of networks of theta neurons

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

$\begin{array}{cc}\dot{\theta }=(1-\cos \theta )+(1+\cos \theta )\eta ,& \left(1\right)\end{array}$

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:

$\begin{array}{cc}\stackrel{.}{{\theta }_j}=\left(1-\cos {\theta }_j\right)+\left(1+\cos {\theta }_j\right)\left[{\eta }_j+I_{syn}\right],& \left(2\right)\end{array}$

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

$\begin{array}{cc}{I}_{syn}=\frac{k}{N}{\displaystyle \sum _{i=1}^N{P}_n({\theta }_i)},& \left(3\right)\end{array}$

where P_n(θ) = a_n (1 − cosθ)ⁿ, n ∈ ℕ, and a_n is a normalization constant¹ such that

${\displaystyle {\int }_0^{2\pi }{P}_n(\theta )d\theta }=2\pi .$

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,

$\begin{array}{cc}g(\eta )=\frac{1}{\pi }\frac{\Delta }{{(\eta -{\eta }_0)}^2+{\Delta }^2},& \left(4\right)\end{array}$

where η₀ 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

$\begin{array}{cc}\tilde{z}(t)={\displaystyle \sum _{j=1}^N{e}^{i{\theta }_j}},& \left(5\right)\end{array}$

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

$\begin{array}{cc}\dot{z}=-i\frac{{(z-1)}^2}{2}+\frac{{(z+1)}^2}{2}\left\{-\Delta +i\left[{\eta }_0+k{H}_n(z)\right]\right\},& \left(6\right)\end{array}$

where

$\begin{array}{cc}{H}_n(z)=I_{syn}/k=a_n\left(A_0+{\displaystyle \sum _{q=1}^n{A}_q(z^q+z^{*q})}\right),& \left(7\right)\end{array}$

$\begin{array}{cc}{A}_q={\displaystyle \sum _{j,m=0}^n{\delta }_{j-2m,q}{Q}_{jm}},& \left(8\right)\end{array}$

and

$\begin{array}{cc}{Q}_{jm}=\frac{{(-1)}^{j-2m}n!}{2^{j}m!(n-j)!(j-m)!}.& \left(9\right)\end{array}$

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, η₀, 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:

$\begin{array}{cc}\begin{array}{l}{\dot{\theta }}_{1,j}=1+{\eta }_{1,j}-(1-{\eta }_{1,j})\cos {\theta }_{1,j}+a_n(1+\cos {\theta }_{1,j})\\ \left[\frac{k_{11}}{N_1}{\displaystyle \sum _{p=1}^{N_1}{(1-\cos {\theta }_{1,p})}^n}+\frac{k_{12}}{N_2}{\displaystyle \sum _{q=1}^{N_2}{(1-\cos {\theta }_{2,q})}^n}\right]\text{}\text{},\\ {\dot{\theta }}_{2,j}=1+{\eta }_{2,j}-(1-{\eta }_{2,j})\cos {\theta }_{2,j}+a_n(1+\cos {\theta }_{2,j})\\ \left[\frac{k_{21}}{N_1}{\displaystyle \sum _{p=1}^{N_1}{(1-\cos {\theta }_{1,p})}^n}+\frac{k_{22}}{N_2}{\displaystyle \sum _{q=1}^{N_2}{(1-\cos {\theta }_{2,q})}^n}\right]\text{}\text{},\end{array}& \left(10\right)\end{array}$

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 η₁, η₂ and widths Δ₁, Δ₂, 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}$₁(t), $\tilde{z}$₂(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:

$\begin{array}{cc}\begin{array}{l}{\dot{z}}_1=-i\frac{{(z_1-1)}^2}{2}+\frac{{(z_1+1)}^2}{2}\\ \left\{-{\Delta }_1+i\left[{\eta }_1+k_{11}{H}_n(z_1)+k_{12}{H}_n(z_2)\right]\right\},\\ {\dot{z}}_2=-i\frac{{(z_2-1)}^2}{2}+\frac{{(z_2+1)}^2}{2}\\ \left\{-{\Delta }_2+i\left[{\eta }_2+k_{21}{H}_n(z_1)+k_{22}{H}_n(z_2)\right]\right\}.\end{array}& \left(11\right)\end{array}$

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₁, N₂ → ∞.

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₁₂ = 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₂₁ ≠ 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

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

Writing the governing equation of population 2 as

$\begin{array}{cc}{\dot{z}}_2=-i\frac{{(z_2-1)}^2}{2}+\frac{{(z_2+1)}^2}{2}\left\{-{\Delta }_2+i\left[{\eta }_{eff}+k_{22}H(z_2)\right]\right\}& \left(12\right)\end{array}$

with

$\begin{array}{cc}{\eta }_{eff}\equiv {\eta }_2+k_{21}H(z_1),& \left(13\right)\end{array}$

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 η₂, the inter-population coupling k₂₁, and the state of the driver z₁.

Note that η_eff depends linearly on both η₂ and k₂₁ and non-linearly on the driver's state z₁ through H(z₁). Additionally, η_eff may be time-dependent if population 1 exhibits a CPW state, since in that case z₁ 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₂₂ > 0, and the “inhibitorily coupled” response system has k₂₂ < 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

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: η₁ = −0.2, Δ₁ = 0.1, k₁₁ = −2, and k₂₂ = 9. (B) The one-dimensional bifurcation diagram showing the asymptotic values of y₂ = Im(z₂) vs. k₂₁. 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 η₂ = −10 and Δ₂ = 0.5.

FIGURE 3

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: η₁ = −0.2, Δ₁ = 0.1, k₁₁ = −2, and k₂₂ = −9. (B) The one-dimensional bifurcation diagram showing the asymptotic value of x₂ = Re(z₂) vs. k₂₁. 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 η₂ = 5 and Δ₂ = 0.5.

FIGURE 4

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

FIGURE 5

Figure 5. (A) The non-linear behavior of η_eff as a function of k₁₁ 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₂ = Im(z₂) vs. k₁₁. 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: η₁ = −0.05, Δ₁ = 0.05, η₂ = 5, Δ₂ = 0.5, and k₂₂ = −9. The inter-population coupling is fixed at k₂₁ = 3.5.

4.1. Driver on a Macroscopic Equilibrium

We begin by fixing the driving population's parameters at η₁ = −0.2, Δ₁ = 0.1, and k₁₁ = −2, which corresponds to a PSR state. Thus, z₁ 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₂₁. From Equation (13), η_eff varies linearly with respect to k₂₁.

4.1.1. Excitatorily-coupled response system

We set the response system's internal coupling to k₂₂ = 9, and show in Figure 2A the two-parameter bifurcation diagram of the response system with respect to Δ₂ 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 η₂ = −10 and Δ₂ = 0.5. Thus, for k₂₁ = 0, η_eff = η₂, and the response system is situated at the solid black point marked in Figure 2A. As k₂₁ 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₂ = Im(z₂)] changes with respect to k₂₁, 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₂₁ 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₂₁ 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₂₂ = −9, i.e., inhibitory, and η₂ = 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 Δ₂ 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₂ = Re(z₂)] changes with respect to k₂₁. As before, η_eff increases linearly as k₂₁ 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₂₁ increases.

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

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

We begin with the case of the excitatorily-coupled response system considered above, with η₂ = −10, Δ₂ = 0.5, and k₂₂ = 9, and choose the remaining driver parameters to be η₁ = −0.05 and Δ₁ = 0.05. Figure 4A shows the non-linear behavior of η_eff as k₁₁ is varied. Even though we are considering k₁₁ 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₁₁ changes, η_eff moves back and forth along the dotted line non-linearly. In particular, Figure 4A shows that for very negative values of k₁₁, η_eff is near −5, which corresponds to a point in Figure 2A to the right of the SN curves. As k₁₁ 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 Δ₂ = 0.5 in Figure 2A crosses the SN curves).

Figure 4B shows the behavior of the asymptotic state of the response system [y₂ = Im(z₂)] as a function of k₁₁. This shows that as k₁₁ 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₁₁ < 0 reversed. Note that the two regions are not symmetrical. This is due to the non-symmetric behavior of η_eff as k₁₁ 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 η₂ = 5 and k₂₂ = −9. Figure 5A shows how η_eff changes as k₁₁ 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₁₁ 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₁₁ < 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₁₁ is increased. Note also the presence of an attracting limit cycle CPW state in intervals of both positive and negative k₁₁.

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

4.2.1. Periodic behavior in the response system

We begin by considering the excitatorily coupled response system, with Δ₂ = 0.5 and k₂₂ = 9, but with η₂ = −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₂₁ is increased from zero to eight, η_eff oscillates back and forth along the horizontal line in Figure 2A at Δ₂ = 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₂ = Re(z₂) vs. k₂₁. As k₂₁ 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

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₂₁. The curves are local maxima and minima of x₂ = Re(z₂). The driver parameters are η₁ = 10.75, Δ₁ = 0.5, and k₁₁ = −9, and the response parameters are η₂ = −20, Δ₂ = 0.5, and k₂₂ = 9. (B) Slightly more complicated periodic behavior obtained at the same parameters, except with η₂ = −5. The curves are local maxima and minima of y₂ = Im(z₂).

We now change the response system such that η₂ = −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₂₁ increases, η_eff oscillates back and forth along the Δ₂ = 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₂ = Im(z₂) vs. k₂₁. Figure 7 shows y₂ vs. x₂ plots of the periodic orbits at k₂₁ = 6 (upper panels) and k₂₁ = 10 (lower panels). As k₂₁ 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₂₁ ≈ 2.5. The amplitude subsequently decreases to a minimum near k₂₁ ≈ 7.2, and then slowly increases again. Note that the four curves in Figure 6B for k₂₁ ∈ [0, 7.2] correspond to two pairs of alternating local maxima and minima in the time series of y₂, as shown in Figure 7B.

FIGURE 7

Figure 7. The response system's behavior at parameters corresponding to Figure 6B at k₂₁ = 6 (A,B) and k₂₁ = 10 (C,D), with z₂ = x₂ + iy₂.

Interestingly, near k₂₁ ≈ 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₂₁ ⪆ 7.2, and the two additional local maxima and minima in the time series of y₂ in Figure 7D.

4.2.2. Multistability in the response system

Continuing with the excitatorily coupled response system (with k₂₂ = 9 > 0), we set η₂ = −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₂₁ increases, η_eff again sweeps back and forth along the horizontal line at Δ₁ = 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₂ vs. k₂₁ for this case. The first feature to emerge as k₂₁ increases from zero is a simple periodic orbit whose amplitude increases, similar to the example in Figure 6A. At k₂₁ ≈ 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₂ vs. x₂ plots of these two limit cycles at k₂₁ = 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

Figure 8. Multistability in the response system driven by a CPW state of the driver. (A) Local maxima and minima of y₂ = Re(z₂) vs. the inter-population coupling k₂₁. (B) y₂ vs. x₂ plots showing two coexisting limit cycles of the response system at k₂₁ = 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₂₁ = 1.5. Green curves are coexisting limit cycles of the response system when coupled to the driver. Parameters are: η₁ = 10.75, Δ₁ = 0.5, k₁₁ = −9; η₂ = −10, Δ₂ = 0.5, k₂₂ = 9.

Figure 8C shows, in black, the asymptotic states of y₂ vs. η_eff for fixed values of η_eff, with k₂₁ = 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 η₂ = 5, Δ₂ = 0.5, and k₂₂ = −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₂₁ increases, η_eff sweeps across the same horizontal line at Δ₂ = 0.5 with increasing amplitude and centroid, initially crossing just the left SN bifurcation curve. At k₂₁ ≈ 5.2, η_eff begins sweeping across the homoclinic and the Andronov-Hopf bifurcation curves. Eventually, for sufficiently large k₂₁, η_eff sweeps across all four bifurcation curves (SN, AH, HC, and SN).

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

FIGURE 9

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₂ = Re(z₂) vs. the inter-population coupling k₂₁. (B) Magnification of the chaotic region (top), with a plot of the largest two Lyapunov exponents. Parameters are: Δ₁ = 0.5, k₁₁ = −9, η₁ = 10.75; Δ₂ = 0.5, k₂₂ = −9, η₂ = 5.

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

The y₂ vs. x₂ plot of the chaotic attractor present at k₂₁ = 5.296, for which the largest Lyapunov exponent is approximately 0.2118, is shown in Figure 10A.

FIGURE 10

Figure 10. (A) Chaotic (y₂ vs. x₂) and (B) Quasiperiodic (y₂ vs. x₁ vs. y₁) attractors in the response system driven by a CPW state of the driver. Parameters are: k₁₁ = k₂₂ = −9, with (A) η₂ = 5, Δ₁ = Δ₂ = 0.5 and k₂₁ = 5.296, and (B) η₂ = 10.75, Δ₁ = 0.5, Δ₂ = 0.3, and k₂₁ = 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 Δ₂ = 0.3. As the inter-population coupling strength k₂₁ is increased, various phase-locked and quasiperiodic states are seen. An example of quasiperiodic behavior in the response system for k₂₁ = 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-\cos (}x){)}^n=n!/(2n-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

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