Neurons in the brain, process information through diverse neural dynamics emergent from interactions among neurons via synapses. How neural networks formed by the interactions can generate functional dynamics, such as oscillatory or synchronized activity and more specifically, cross-frequency coupling (CFC), largely remains to be explored.

A variety of oscillatory phenomena in the brain have been studied often through the measurement of the local field potential or an electroencephalogram, both of which include macroscopic information of neural networks (cell assemblies), different from single-unit recordings. The recorded macroscopic neural activity can provide useful indices of distinctive brain functions, where such indices show a possibility that the oscillatory phenomena correlate with the brain functions (Buzsáki and Draguhn, 2004). Further, the recorded neural activity represents various types of oscillatory waveforms and can be categorized according to the frequency, whose bands, for example, theta and gamma bands, can temporally coexist in the same or different brain regions (Steriade, 2001; Csicsvari et al., 2003). Among the sorted frequency bands, neighboring bands recorded in the same brain region may differ with the brain functions (Klimesch, 1999; Kopell et al., 2000; Engel et al., 2001; Csicsvari et al., 2003). The oscillatory frequency is nearly inversely proportional to the power in general as observed from the power spectrum (Freeman et al., 2000). This inversely proportional property of the power suggests that spatially widespread slow oscillations can modulate the local neural activity (Steriade, 2001; Csicsvari et al., 2003; Sirota et al., 2003). Specifically, when the slow and fast oscillatory components interact with each other, CFC phenomena emerge.

The mechanism underlying the oscillatory phenomena is thought to be based on the following three elements: neurons, synapses, and connectivity, all of which are necessary to constitute neural networks and to generate diverse macroscopic oscillations such that sensory input reflecting external environmental input can be correctly encoded in the oscillations.

First, a neuron generates spike trains stochastically and irregularly; that is, the neuron itself may not possess the reliability to generate clear rhythms, although evidence that a single neuron can produce different rhythms has been obtained on the other hand. There exist, at least, two types of neurons in the cerebral cortex, namely, excitatory and inhibitory neurons, which exhibit different response properties (Lux and Pollen, 1966; Connors et al., 1982; Connors and Gutnick, 1990; Kawaguchi and Kubota, 1997; Nowak et al., 2003; Tateno et al., 2004). In particular, an excitatory neuron typically shows a low firing rate whereas an inhibitory neuron shows a high firing rate (Tateno et al., 2004).

The second fundamental element for the oscillatory phenomena is a synapse, intermediating between neurons with its transmission efficacy. It had been considered that the synaptic efficacy is nearly constant over time or changing very slowly; more specifically, peaks of the synaptic current, induced by releases of neurotransmitters from presynaptic vesicles, had been thought to be independent of the timing of neural firings. During the last few decades, many reports have shown, however, that synapses involve fast plasticity mechanisms, enabling neurons with such synapses to generate more flexible dynamics compared to those with ‘static’ synapses. In particular, synapses that exhibit rapid changes in the coupling strength between neurons with a short-term plasticity mechanism are called dynamic synapses (Markram and Tsodyks, 1996; Markram et al., 1998; Thomson, 2000; Wang et al., 2006). Two types of dynamic synapses exist, namely, short-term depression and facilitation synapses, which are characterized by the transiently decreasing ratio of releasable neurotransmitters and by the transiently increasing calcium concentration in presynaptic terminals, respectively. The synaptic transmission efficacy decreases or increases depending on the ratio of the two time constants associated with the recovery from the aforementioned transient decrease or increase (Markram et al., 1998; Thomson, 2000). The distribution of the dynamic synapses in the brain differs among brain regions. For example, many depression synapses appear in the parietal lobe, while many facilitation synapses appear in the prefrontal lobe (Wang et al., 2006).

The third element for generating the oscillatory phenomena is connectivity, which enables neurons to interact with each other and to form a neural network via synapses. Although each neuron fires stochastically and the resulting spike trains do not show clear rhythmicity, the neural network can macroscopically generate rhythmic oscillations. When excitatory neurons aggregate to be structured as a cell assembly receiving input from inhibitory neurons, rhythmic dynamics on networks macroscopically appears and possesses reliability (Yoshimura et al., 2005). In particular, interactions between two types of networks, i.e., excitatory and inhibitory networks, would underlie the emergence of CFC (Aru et al., 2015; Hyafil et al., 2015). Notably, a network composed of the excitatory and inhibitory subnetworks can generate either slow or fast oscillations. Thus, the interaction between this network and another, which is also composed of the excitatory and inhibitory subnetworks, can be one origin of CFC; this coupling structure has been called as bidirectional coupling (Hyafil et al., 2015).

Here, stochastic neural network models have been utilized to elucidate the relationship between irregular spike trains and rhythmic macroscopic oscillations. The process of the macroscopic oscillations generated from such irregular dynamics can be explained by a network receiving the following two types of input: strong external noise and strong recurrent inhibition (Brunel and Hakim, 1999; Brunel and Wang, 2003). The fast repetition of such noisy input causes short-term synaptic depression, and then, the network induces the destabilization of attractors (Cortes et al., 2006) and chaotic oscillations (Marro et al., 2007).

Additionally, the stochastic neural network models with dynamic synapses have been intensively investigated. Synaptic depression triggers a large fluctuation in sustained periods between the up and down states (Mejias et al., 2010), and the level of the synaptic depression changes the property of the sustained activities of these two different states (Benita et al., 2012). Moreover, the synaptic depression contributes to the destabilization of network activity, the generation of an oscillatory state, and the spontaneous state transitions among multiple patterns in an associative memory network (Katori et al., 2013). Additionally, the synaptic depression can be a suitable mechanism to explain critical avalanches in self-organized neural networks (Bonachela et al., 2010). On the other hand, synaptic facilitation enhances the working memory function (Mongillo et al., 2008). Both synaptic depression and facilitation are related to the storage capacity of attractor neural networks (Bibitchkov et al., 2002; Torres et al., 2002; Matsumoto et al., 2007; Mejias and Torres, 2009; Otsubo et al., 2011; Mejias et al., 2012). Bibitchkov et al. have shown that depression synapses may reduce the storage capacity (Bibitchkov et al., 2002). Torres et al. investigated this negative effect of depression synapses on the storage capacity in general (Torres et al., 2002). However, Matsumoto et al. argued that the storage capacity is not influenced by synaptic depression, when noise is not considered in the network (Matsumoto et al., 2007). Otsubo et al. reported that a network with both depression synapses and the noise would reduce the storage capacity (Otsubo et al., 2011). Then, Mejias and Torres found that the combination of depression and facilitation synapses can enhance the storage capacity (Mejias and Torres, 2009) and furthermore, Mejias et al. generated phase diagrams, indicating that synaptic facilitation enlarges the memory phase region (Mejias et al., 2012). Further, dynamic synapses play a role in stochastic resonance, where a weak input signal to a network can be detected in an output signal under certain conditions (Pantic et al., 2003; Mejias and Torres, 2008, 2011; Torres et al., 2011; Pinamonti et al., 2012; Torres and Marro, 2015). Pantic et al. have shown that a neuron with depression synapses is capable of detecting noisy input signals with a wider frequency range, compared to one with static synapses, under a certain firing threshold (Pantic et al., 2003). Mejias and Torres found that the inclusion of facilitation dynamics in depression synapses would enhance the detection performance (Mejias and Torres, 2008). Moreover, they have shown that such combination of depression and facilitation synapses can generate two suitable noise levels to detect input signals; where it has been argued that this bimodal resonance is caused by the interplay between the adaptively varying firing threshold and the dynamic synapses (Mejias and Torres, 2011). Torres et al. demonstrated that a model with this interplay can predict experimental data of stochastic resonance (Torres et al., 2011). Furthermore, Pinamonti et al. demonstrated that stochastic resonance is well enhanced near phase transitions among patterns in an associative memory network (Pinamonti et al., 2012). Then, Torres and Marro generated a detailed phase diagram embedding many patterns associated with stochastic resonance, such that multiple noise levels are well responsible for optimizing input signals (Torres and Marro, 2015). Additionally, it has been reported that the combination of depression and facilitation synapses contributes to flexible information representation, from the viewpoint of both data analysis and mathematical modeling (Katori et al., 2011).

In particular, the stochastic neural network models and the corresponding mean field models have been effectively used for analyzing the properties of neural networks, including those with dynamic synapses (Pantic et al., 2002; Torres et al., 2007; Igarashi et al., 2010; Katori et al., 2012). While Pantic et al. and Torres et al. have derived mean field models by taking the population average, with respect to stochastic variables, based on the assumption of ergodicity (Pantic et al., 2002; Torres et al., 2007), Igarashi et al. and Katori et al. have recently introduced dynamic mean field theory in which two types of mean field models, i.e., microscopic and macroscopic mean field models, are in turn derived (Igarashi et al., 2010; Katori et al., 2012).

However, one of the key oscillatory phenomena, CFC described above, has not been well understood, although this coupling phenomenon could contribute to complex information processing in the brain, thanks to the presence of oscillations with two different timescales. On the other hand, much attention to phase-amplitude CFC has been attracted, because it has been suggested that this kind of CFC plays a crucial role in adjusting neural communications among distant brain regions (Canolty et al., 2006; Jansen and Colgin, 2007; Canolty and Knight, 2010). To elucidate possible mechanisms generating the CFC, the continuous-time neuron models have been investigated so far (Malerba and Kopell, 2012; Fontolan et al., 2013), but little is known about discrete-time neuron models, which have been shown to be suitable for simulating, reproducing, and predicting neural phenomena in the brain (Rulkov, 2002; Ibarz et al., 2011). Thus, in this study, we extend the previously proposed discrete-time network model (Katori et al., 2012), which only includes excitatory neurons, to that including also inhibitory neurons, and clarify the bifurcation structure underlying the phase-amplitude CFC.

In this study, we hypothesize that the CFC is generated by the bidirectional coupling between the two subnetworks as described above, where one subnetwork includes an excitatory population while the other includes an inhibitory population, both of which receive input from another excitatory or inhibitory population and generate slow or fast oscillations (Figure 1A) (White et al., 2000; Kramer et al., 2008; Roopun et al., 2008; Hyafil et al., 2015). Each subnetwork can be regarded as a pure excitatory or a pure inhibitory network because the unidirectional input may be equivalent to change of the neural firing threshold. Therefore, hereinafter, we call the above two subnetworks as an excitatory network/subnetwork and an inhibitory network/subnetwork, respectively, for the sake of simplicity (Figure 1A). Accordingly, we focus on a stochastic network composed of excitatory and inhibitory neurons with dynamic synapses. Furthermore, the proposed model considers the decay process of the synaptic current. We analyze a macroscopic mean field model reproducing the overall network dynamics associated with the stochastic model. Upon the adjustment of parameters specifying the properties of the synaptic current and dynamic synapses, rich bifurcation structures of the network dynamics are expected to be found. In the following sections, first, we describe a network model with stochastic neurons connected via dynamic synapses. Then, we derive a macroscopic mean field model capturing the macroscopic dynamics of the network model. Subsequently, we analyze the bifurcation structures of the present model and illustrate various solutions included in the dynamical systems of not only an excitatory or an inhibitory network, but also a network composed of both the excitatory and the inhibitory subnetworks. Finally, we discuss the dynamical properties revealed from the standpoint of neuroscience, and consider possible future directions of this research.

This section consists of three parts. The first part describes the mechanism of signal transmission on synapses with short-term plasticity. The second part explains the stochastic model, composed of excitatory and inhibitory binary neurons and dynamic synapses. The third part introduces the mean field theory in which the stochastic model is converted into the corresponding microscopic and macroscopic mean field models. The microscopic mean field approximation is used for extracting the average neural activity on realization of stochastic variables. In contrast, the macroscopic mean field approximation is used for extracting the average neural activity over population dynamics, and thereby, a low-dimensional discrete-time dynamical system is derived, such that we can identify bifurcation structures.

In this section, we analyze a variety of bifurcation structures arising from the proposed model by considering τaξ, J0ξξ, J0ξη, I^ξ, and τRξ/τFξ to be bifurcation parameters and T^ξ and Useξ to be constants; we set T^ξ = 0.8 and Useξ=0.1, respectively, because the activation function gβξ[·] used here is an idealized sigmoid form different from the experimentally known frequency-current relationship (Tateno et al., 2004) and Useξ merely determines the steady value of the calcium concentration. The bifurcation analysis focuses on two types of dynamic synapses, namely short-term depression (τRξ/τFξ≫1) and facilitation (τRξ/τFξ≪1) synapses, where τRξ has been fixed at 70 and τFξ ranges between 1 and 140 for the sake of simplicity. In this study, first, we individually investigate an excitatory and an inhibitory networks, by setting J0EI=0 and J0EE≥0 for the excitatory network and J0IE=0 and J0II<0 for the inhibitory network. Next, we analyze a network composed of both the excitatory and the inhibitory subnetworks. Throughout the study, when the stochastic model is simulated to be compared with the corresponding macroscopic mean field model, we use Nξ=104. In the following, we omit the superscripts attached to the variables and the parameters for the sake of simplicity when examining the aforementioned two networks independently.

To help the readers to understand the bifurcation analysis results below, we summarize the Result section briefly here referring to Figures 3–12. First, we analyze the excitatory and inhibitory networks independently, so that a parameter space generating the oscillatory state becomes clear (Figures 3, 5). The oscillation on the inhibitory network is faster than that on the excitatory network (Figure 6), and the frequency of both oscillations changes depending on parameters I, τ_a, and τ_R/τ_F (Figures 3, 5). Next, we analyze a network composed of the above two subnetworks, so that a parameter space generating two types of phase-amplitude CFCs becomes clear (Figure 7). The two CFC modes differ in the underlying attractors (Figures 8, 10, 11) and in the modulation properties (Figure 12), depending on parameters J0EI and J0IE (Figure 7), but this coupling phenomenon disappears if inhibitory input is large enough (Figures 7, 9). While these analyses are applied to the macroscopic mean field model, the stochastic model also yields the consistent results (Figure 4).

The excitatory and the inhibitory networks, each of which reflects the property of depression synapses with τ_R/τ_F = 11.7, exhibit distinctive oscillatory states, as shown in the (J₀, I) phase diagram (Figure 3A). The oscillatory state on the excitatory network (OSE), which was already found in the previous model of Katori et al. (2012), changes into the steady state (the fixed point) via the Neimark-Sacker (NS) bifurcation (Kuznetsov, 2004), when J₀ decreases/increases or I increases from the region of the OSE state. On the other hand, when J₀ decreases and I increases from the region of the steady state, the oscillatory state on the inhibitory network (OSI) is likewise generated by the NS bifurcation; this OSI state has been newly observed here. It is clear from the bifurcation diagrams that these oscillatory states emerge via the NS bifurcation from the steady state (see Figures 3B,C). The bifurcation diagram, with respect to J₀ on (I, τ_a) = (−1, 2.5) (Figure 3B), shows the OSE state emergent from the steady state. The mean neural activity, m₀, in the steady state, increases with an increase in J₀, whereas the steady state destabilizes via the first NS bifurcation at J₀ = 1.63 and then, the OSE state appears. Because of the second NS bifurcation at J₀ = 3.48, the OSE state disappears, and accordingly, the steady state reappears. In contrast, the bifurcation diagram, with respect to J₀ on (I, τ_a) = (1, 2.5) (Figure 3C), shows the OSI state emergent from the steady state. The mean neural activity, m₀, in the steady state decreases with a decrease in J₀, whereas the steady state destabilizes via the NS bifurcation at J₀ = −4.73 and accordingly, the OSI state appears. The OSE and the OSI states, generated from the steady state via the NS bifurcation, likewise appear on the network with facilitation synapses. The aforementioned two bifurcation diagrams, with respect to J₀, show good agreement with those generated from the stochastic model (Figures 4A,B).

The model has been proposed as a discrete-time system, such that time is represented by an arbitrary unit to flexibly describe real data. To characterize oscillatory time-scale from this kind of time, first, we converted time courses generated from the model into the power spectrum, where 4096 time steps were used in calculation for a time course. Second, the frequency of the first typical peak in the spectrum was translated into the corresponding time steps; we call this time step value as “period”. Based on the period, we can say that the oscillation on the inhibitory network tends to be faster than that on the excitatory network (Figure 3A). Specifically, the period of the OSI state ranges from 4.99 to 6.00 time steps, whereas that of the OSE state ranges from 33.9 to 78.8 time steps.

Figures 3D,E show the bifurcation structures of the excitatory and the inhibitory networks, respectively, and represent the quasi-periodic oscillations on the invariant closed curve for both the OSE and the OSI states. The (τ_a, I) phase diagram in Figure 3D shows the distribution of the OSE state, where the period of the OSE state ranges from 44.0 to 75.9 time steps and increases as τ_a and I decrease. In contrast, the (τ_a, I) phase diagram in Figure 3E shows the distribution of the OSI state, where the period of the OSI state ranges from 3.94 to 8.00 time steps and tends to increase as τ_a and I increase.

The property of the dynamic synapses similarly affects the period on both the excitatory and the inhibitory networks (Figure 5). As the influence of short-term depression on the excitatory network increases (as τ_R/τ_F increases), the period for the OSE state increases, and the network can oscillate with less inhibitory input (Figure 5A). In contrast, τ_R/τ_F does not have much effect on the period for the OSI state (Figure 5B).

The typical OSE and OSI states represent the relatively slow and the fast oscillations, respectively, and are observed in both the stochastic and the macroscopic mean field models (Figure 6). The population averages of stochastic variables, s0(t)[=(1/N)∑iNsi(t)], a0(t)[=(1/N)∑iNai(t)], x0(t)[=(1/N)∑iNxi(t)], and u0(t)[=(1/N)∑iNui(t)], correspond to macroscopic variables (Figures 6A,D). Most excitatory (inhibitory) neurons fire together at each time step at a low (high) frequency, such that the macroscopic variables exhibit oscillations with large amplitude. The attractor, composed of the variables, is the closed curve for both the excitatory and the inhibitory networks (Figures 6C,F), but fundamental frequency components differ between these networks (Figures 6B,E). The dynamics of such macroscopic variables shows good agreement with that of the stochastic variables, in terms of: the time course (Figures 6A,D); the power spectrum (Figures 6B,E); and the trajectory in the state space (Figures 6C,F).

Figure 7 shows the dynamical structure of a network composed of the abovementioned excitatory and the inhibitory subnetworks, where the following parameter set, corresponding to the typical OSE and OSI states described above, has been used: (J0EE,IE,τaE)=(2,-1,2.5) for the excitatory subnetwork, (J0II,II,τaI)=(-10,1,12.5) for the inhibitory subnetwork, and τRξ/τFξ=11.7 for the depression synapses. This network with (J0EI,J0IE)=(0,0) is a direct product system, composed of the excitatory and the inhibitory subnetworks, so that a part of the system is equivalent to the previous model of Katori et al. (2012). However, the network with (J0EI,J0IE)≠(0,0) shows the following four types of distinctive neural dynamics: the steady state (SS); the oscillatory state with a single frequency component on a closed curve (OS1C); that with two frequency components on a two-dimensional torus (OS2T); and that with two frequency components on a closed curve (OS2C), where the OS2T and OS2C states did not appear in the previous model. It is clear from the number of zero-Lyapunov exponents whether the attractor underlying the network is a closed curve or a two-dimensional torus (see the (J0EI,J0IE) phase diagram in Figure 7A). The network dynamics involves only one zero-exponent for the emergence of the OS1C or the OS2C state, whereas the dynamics involves two zero-exponents for the emergence of the OS2T state. The dynamics is not affected by the property of the dynamic synapses characterized by parameter τRξ/τFξ, which merely changes the period of oscillations on subnetwork ξ.

Here we use the term MTd to investigate the bifurcation structure distinguishing the above four states clearly. Then, the MT1 arising from the model corresponds to the OS1C or the OS2C state, while the MT2 corresponds to the OS2T state. While the OS1C and the OS2C states are generated from the SS state via the NS bifurcation, these states change into the OS2T state via the two types of bifurcation, namely, the NS bifurcation of MT1 (MT1NS) and the saddle-node cycle bifurcation of MT1 (MT1SNC) (Kamiyama et al., 2014; Komuro et al., 2016). Because the MT1NS bifurcation here is subcritical, a saddle-node bifurcation of MT2 (MT2SN) intermediates between the OS1C and the OS2T states (Komuro et al., 2016), where there exists a specific hysteresis region between these states. We show below the qualitative difference between the states, before and after the bifurcation of the MT1 and the MT2, by observing the bifurcation diagrams generated from quasi-periodic points collected in the slice.

First, we consider the bifurcations intermediating between the OS1C and the OS2T states. The bifurcation diagram of trajectories in the slice, with respect to J0EI on J0IE=2, shows a transition between the OS1C and the OS2T states, where section A0I(t)=1.7 and width ϵ = 0.001 have been used for the slice (Figure 7B). As J0EI increases from J0EI=-0.145, the ST0, corresponding to the OS1C state, destabilizes and immediately an oscillation starts along with a jump to the ST1, corresponding to the OS2T state; this abrupt change is attributed to the subcritical bifurcation. The bifurcation property becomes clearer by a simultaneous observation of both trajectories in the slice, applied to a state just before the bifurcation, and that just after the bifurcation. This simultaneous analysis clarifies that the ST1 appears sufficiently far from the ST0 when the bifurcation occurs, where J0EI=-0.142 for the OS1C state, and J0EI=-0.139 for the OS2T state (Figure 7F). In contrast, as J0EI decreases from J0EI=-0.135, the stable ST1 and the saddle ST1 collide, and the ST0 appears. By combining the aforementioned slice analysis with the Lyapunov exponent analysis, we can identify the bifurcation types more clearly (Figure 7D). The negative exponents with multiplicity of two approach zero simultaneously as J0EI increases from J0EI=-0.145; this is a property of the NS bifurcation. However, only one exponent starts to decrease as J0EI decreases from J0EI=-0.135; this is a property of the saddle-node bifurcation. Taken together, the bifurcation from the OS1C to OS2T states can be identified as the subcritical NS bifurcation of ST0 (ST0NS), while that from the OS2T to OS1C states as the saddle-node bifurcation of ST1 (ST1SN), within the slice (Komuro et al., 2016). Thus, these two bifurcations can be interpreted as the MT1NS and the MT2SN bifurcations, respectively, outside the slice (Komuro et al., 2016).

Next, we consider the bifurcation intermediating between the OS2C and the OS2T states. The bifurcation diagram of trajectories in the slice, with respect to J0IE on J0EI=-0.05, shows a transition between the OS2C and the OS2T states, where section A0E(t)=0.8 and width ϵ = 0.001 have been used for the slice (Figure 7C). As J0IE decreases from J0IE=9, the ST0, corresponding to the OS2C state, changes into the ST1 at a bifurcation point. It is clear that there does not exist a hysteresis region between the OS2C and the OS2T states, because the ST1 changes into the ST0 at the same bifurcation point as J0IE increases from J0IE=0. By observing both a state just before the bifurcation and that just after the bifurcation simultaneously, we can find that the ST1 rightly covers the ST0, where J0IE=3.12 for the OS2C state and J0IE=3 for the OS2T state (Figure 7G). The bifurcation involving this feature has been recently identified as the saddle-node cycle bifurcation of ST0 (ST0SNC), within the slice (Figure 8A) (Komuro et al., 2016), or as the MT1SNC bifurcation, outside the slice (Figure 8B) (Komuro et al., 2016). Before the ST0SNC (MT1SNC) bifurcation occurs, there exists a pair of a stable ST0 (MT1) and a saddle ST0 (MT1) and therefore, an unstable set of the saddle ST0 (MT1) generates a one-dimensional (two-dimensional) torus. The ST0SNC (MT1SNC) bifurcation, where the stable ST0 (MT1) and the saddle ST0 (MT1) collide, stabilizes the unstable set and accordingly, an ST1 (MT2) appears. The mechanism of this MT1SNC bifurcation can be likewise verified by the Lyapunov exponent analysis (Figure 7E); only one negative exponent approaches zero as J0IE decreases from J0IE=9.

Figures 9–11 show typical dynamics emergent from the OS1C, the OS2T, and the OS2C states, respectively, where each is observed in both the stochastic and the macroscopic mean field models. The population averages of stochastic variables, s0ξ(t)[=(1/Nξ)∑iNξsiξ(t)], a0ξ(t)[=(1/Nξ)∑iNξaiξ(t)], x0ξ(t)[=(1/Nξ)∑iNξxiξ(t)], and u0ξ(t)[=(1/Nξ)∑iNξuiξ(t)], correspond to macroscopic variables. The dynamics of the macroscopic variables shows good agreement with that of the stochastic variables, in terms of: the time course (Figures 9A,B, 10A,B, 11A,B); the power spectrum (Figures 9C, 10C, 11C); the trajectory in the state space (Figures 9D, 10D, 11D); and its distribution in the slice (Figures 9E, 10E, 11E).