Home | About | Journals | Submit | Contact Us | Français |

**|**Front Comput Neurosci**|**v.11; 2017**|**PMC5371682

Formats

Article sections

- Abstract
- 1. Introduction
- 2. Materials and methods
- 3. Results
- 4. Discussion
- Author contributions
- References

Authors

Related links

Front Comput Neurosci. 2017; 11: 18.

Published online 2017 March 30. doi: 10.3389/fncom.2017.00018

PMCID: PMC5371682

Edited by: Carlo Laing, Massey University, New Zealand

Reviewed by: Joaquín J. Torres, University of Granada, Spain; Gennady Cymbalyuk, Georgia State University, USA

Received 2016 August 23; Accepted 2017 March 13.

Copyright © 2017 Sase, Katori, Komuro and Aihara.

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

We investigate a discrete-time network model composed of excitatory and inhibitory neurons and dynamic synapses with the aim at revealing dynamical properties behind oscillatory phenomena possibly related to brain functions. We use a stochastic neural network model to derive the corresponding macroscopic mean field dynamics, and subsequently analyze the dynamical properties of the network. In addition to slow and fast oscillations arising from excitatory and inhibitory networks, respectively, we show that the interaction between these two networks generates phase-amplitude cross-frequency coupling (CFC), in which multiple different frequency components coexist and the amplitude of the fast oscillation is modulated by the phase of the slow oscillation. Furthermore, we clarify the detailed properties of the oscillatory phenomena by applying the bifurcation analysis to the mean field model, and accordingly show that the intermittent and the continuous CFCs can be characterized by an aperiodic orbit on a closed curve and one on a torus, respectively. These two CFC modes switch depending on the coupling strength from the excitatory to inhibitory networks, via the saddle-node cycle bifurcation of a one-dimensional torus in map (MT1SNC), and may be associated with the function of multi-item representation. We believe that the present model might have potential for studying possible functional roles of phase-amplitude CFC in the cerebral cortex.

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 (Figure1A)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 (Figure1A).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.

Here, we describe the signal transmission mechanism on dynamic synapses with the short-term plasticity. This mechanism can be explained by the following three processes (Figure (Figure1B):1B): First, an action potential is generated on the presynaptic neuron and transmitted to the terminals of the synapses (①). Second, voltage-dependent calcium channels are opened by the action potential, and calcium ions flow into the synaptic terminals via the channels (②). Third, chemical reactions with these calcium ions cause the fusion of the presynaptic membrane and synaptic vesicles, including neurotransmitters, which are released into the synaptic cleft. The released neurotransmitters attach to the postsynaptic membrane, and then, a synaptic current is generated (③). The amplitude of the synaptic current decays with a certain time constant, dependent on the properties of the postsynaptic receptors.

Arrival of many action potentials within a short period of time causes a transient decrease or increase in the efficacy of the synaptic transmission. This is because of changes in the amount of neurotransmitters and in the calcium concentration in the presynaptic terminal. Finally, the synaptic vesicles are retrieved and the neurotransmitters are restored in the reusable synaptic vesicles (Markram and Tsodyks, 1996; Markram et al., 1998; Thomson, 2000; Wang et al., 2006).

We consider a discrete-time neural network model composed of the excitatory (E) and the inhibitory (I) subnetworks, which consist of *N*_{E} excitatory and *N*_{I} inhibitory neurons, respectively. The state of the *i*th neuron belonging to the network ξ (ξ{E, I}) at time *t*, denoted by ${s}_{i}^{\xi}(t)$, indicates either a quiescent state $\left[{s}_{i}^{\xi}(t)=0\right]$ or an active state $\left[{s}_{i}^{\xi}(t)=1\right]$. The state of the neuron is stochastically determined by total input for the neuron. This stochasticity is reflecting the biologically observed noisy neural activity. The evolution of the neuronal state is described as follows:

$$\begin{array}{l}\text{Prob}\left[{s}_{i}^{\xi}\left(t+1\right)=1\right]={g}_{{\beta}^{\xi}}\left[{h}_{i}^{\xi}\left(t\right)\right],\end{array}$$

(1)

$$\begin{array}{l}{g}_{{\beta}^{\xi}}\left[{h}_{i}^{\xi}(t)\right]=\frac{1}{2}\left\{1+tanh\left[{\beta}^{\xi}{h}_{i}^{\xi}(t)\right]\right\},\end{array}$$

(2)

where ${h}_{i}^{\xi}(t)=\sum _{j\ne i}^{{N}_{\xi}}\left[{J}_{ij}^{\xi \xi}{a}_{j}^{\xi}(t)\right]+\sum _{j=1}^{{N}_{\eta}}\left[{J}_{ij}^{\xi \eta}{a}_{j}^{\eta}(t)\right]+{I}^{\xi}$ with η{E, I|η≠ξ}. The variables ${h}_{i}^{\xi}(t)$ and ${a}_{j}^{\xi}(t)$ represent the total input for the *i*th neuron on network ξ and the synaptic activity with short-term plasticity, respectively; 1/β^{ξ} = *T*^{ξ} denotes the noise intensity; ${J}_{ij}^{\xi \xi}$, ${J}_{ij}^{\xi \eta}$, and *I*^{ξ} indicate the weight of the recurrent coupling from the *j*th neuron to the *i*th neuron on network ξ, the weight of the unidirectional coupling from the *j*th neuron on network η to the *i*th neuron on network ξ, and the external input, respectively.

The signal transmission mechanism in synapses with the short-term plasticity can be modeled by the synaptic activity ${a}_{i}^{\xi}(t)$, and two variables, ${x}_{i}^{\xi}(t)$ and ${u}_{i}^{\xi}(t)$, denoting the ratio of releasable neurotransmitters and the calcium concentration, respectively. The synaptic activity, ${a}_{i}^{\xi}(t)$, increases with the presynaptic neural activity. This increase is proportional to the product of ${x}_{i}^{\xi}(t)$ and ${u}_{i}^{\xi}(t)$, which represents the synaptic transmission efficacy (Markram et al., 1998; Mongillo et al., 2008; Mejias and Torres, 2009). If there is no synaptic activation, ${a}_{i}^{\xi}(t)$ converges to its steady state ${a}_{i}^{\xi}(t)=0$ with the time constant ${\tau}_{\mathrm{\text{a}}}^{\xi}$. The ratio ${x}_{i}^{\xi}(t)$ of releasable neurotransmitters decreases with the presynaptic activation; this decrease is proportional to ${u}_{i}^{\xi}(t)$. Then, ${x}_{i}^{\xi}(t)$ returns to its steady state ${x}_{i}^{\xi}(t)=1$ with the time constant ${\tau}_{\mathrm{\text{R}}}^{\xi}$. The variable ${u}_{i}^{\xi}(t)$ increases when a presynaptic neuron is activated, and returns to its steady state ${u}_{i}^{\xi}(t)={U}_{\mathrm{\text{se}}}^{\xi}$ with the time constant ${\tau}_{\mathrm{\text{F}}}^{\xi}$; this increase is proportional to ${U}_{\mathrm{\text{se}}}^{\xi}$. The dynamics described above is summarized as follows:

$$\begin{array}{l}{a}_{i}^{\xi}(t+1)={a}_{i}^{\xi}(t)-\frac{{a}_{i}^{\xi}(t)}{{\tau}_{\mathrm{\text{a}}}^{\xi}}+{s}_{i}^{\xi}(t){x}_{i}^{\xi}(t){u}_{i}^{\xi}(t)/{U}_{\mathrm{\text{se}}}^{\xi},\end{array}$$

(3)

$$\begin{array}{l}{x}_{i}^{\xi}(t+1)={x}_{i}^{\xi}(t)+\frac{1-{x}_{i}^{\xi}(t)}{{\tau}_{\mathrm{\text{R}}}^{\xi}}-{s}_{i}^{\xi}(t){x}_{i}^{\xi}(t){u}_{i}^{\xi}(t),\end{array}$$

(4)

$$\begin{array}{l}{u}_{i}^{\xi}(t+1)={u}_{i}^{\xi}(t)+\frac{{U}_{\mathrm{\text{se}}}^{\xi}-{u}_{i}^{\xi}(t)}{{\tau}_{\mathrm{\text{F}}}^{\xi}}+{U}_{\mathrm{\text{se}}}^{\xi}(1-{u}_{i}^{\xi}(t)){s}_{i}^{\xi}(t),\end{array}$$

(5)

where the ratio ${\tau}_{\mathrm{\text{R}}}^{\xi}/{\tau}_{\mathrm{\text{F}}}^{\xi}$ determines whether the plasticity is short-term depression or facilitation (Wang et al., 2006).

The dynamic mean field model was developed according to the following two steps: In the first step, a microscopic mean field model is derived from the stochastic model by taking the expectation (noise average), with respect to each variable. This expectation cannot be replaced with the population or time average. In the second step, a macroscopic mean field model is derived from the microscopic model by taking the average over the population. By assuming that the recurrent network is fully connected via synapses with the coupling strength of the order of 1/*N*_{ξ}, we obtain an eight-dimensional discrete-time dynamical system, not dependent on the number of neurons, *N*_{ξ}.

First, Equations (1) and (2) are converted into the following forms:

$$\begin{array}{l}\langle {s}_{i}^{\xi}(t+1)\rangle ={g}_{{\beta}^{\xi}}\left[\langle {h}_{i}^{\xi}(t)\rangle \right],\end{array}$$

(6)

$$\begin{array}{l}\langle {h}_{i}^{\xi}(t)\rangle ={\displaystyle \sum _{j\ne i}^{{N}_{\xi}}}\left[{J}_{ij}^{\xi \xi}\langle {a}_{j}^{\xi}(t)\rangle \right]+{\displaystyle \sum _{j=1}^{{N}_{\eta}}}\left[{J}_{ij}^{\xi \eta}\langle {a}_{j}^{\eta}(t)\rangle \right]+{I}^{\xi},\end{array}$$

(7)

where the notation, ·, indicates the noise average. Similarly, the following equations, corresponding to Equations (3) to (5), are obtained:

$$\begin{array}{l}\langle {a}_{i}^{\xi}(t+1)\rangle =\langle {a}_{i}^{\xi}(t)\rangle -\frac{\langle {a}_{i}^{\xi}(t)\rangle}{{\tau}_{\mathrm{\text{a}}}^{\xi}}+\langle {s}_{i}^{\xi}(t){x}_{i}^{\xi}(t){u}_{i}^{\xi}(t)\rangle /{U}_{\mathrm{\text{se}}}^{\xi},\end{array}$$

(8)

$$\begin{array}{l}\langle {x}_{i}^{\xi}(t+1)\rangle =\langle {x}_{i}^{\xi}(t)\rangle +\frac{1-\langle {x}_{i}^{\xi}(t)\rangle}{{\tau}_{\mathrm{\text{R}}}^{\xi}}-\langle {s}_{i}^{\xi}(t){x}_{i}^{\xi}(t){u}_{i}^{\xi}(t)\rangle ,\end{array}$$

(9)

$$\begin{array}{l}\langle {u}_{i}^{\xi}(t+1)\rangle =\langle {u}_{i}^{\xi}(t)\rangle +\frac{{U}_{\mathrm{\text{se}}}^{\xi}-\langle {u}_{i}^{\xi}(t)\rangle}{{\tau}_{\mathrm{\text{F}}}^{\xi}}+{U}_{\mathrm{\text{se}}}^{\xi}\langle (1-{u}_{i}^{\xi}(t)){s}_{i}^{\xi}(t)\rangle .\end{array}$$

(10)

Here, we assume that ${J}_{ij}^{\xi \xi}$ and ${J}_{ij}^{\xi \eta}$ are of the order of 1/*N*_{ξ} and 1/*N*_{η}, respectively; therefore, the correlation between ${s}_{i}^{\xi}(t)$ and ${x}_{i}^{\xi}(t)$ is negligible when *N*_{ξ} → ∞ (Igarashi et al., 2010). Likewise, the correlation between ${s}_{i}^{\xi}(t)$ and ${u}_{i}^{\xi}(t)$ approaches zero as *N*_{ξ} → ∞. Furthermore, in a previous study, it has been found that the correlation between ${x}_{i}^{\xi}(t)$ and ${u}_{i}^{\xi}(t)$ is negligible when the number of neurons is sufficiently large (Katori et al., 2012), while it has been reported that the independence of ${x}_{i}^{\xi}(t)$ and ${u}_{i}^{\xi}(t)$ is maintained, if there is no facilitation (Tsodyks et al., 1998). Consequently, we assume the following for simplicity:

$$\begin{array}{l}\langle {s}_{i}^{\xi}(t){x}_{i}^{\xi}(t){u}_{i}^{\xi}(t)\rangle =\langle {s}_{i}^{\xi}(t)\rangle \langle {x}_{i}^{\xi}(t)\rangle \langle {u}_{i}^{\xi}(t)\rangle ,\end{array}$$

(11)

$$\begin{array}{l}\langle {s}_{i}^{\xi}(t){u}_{i}^{\xi}(t)\rangle =\langle {s}_{i}^{\xi}(t)\rangle \langle {u}_{i}^{\xi}(t)\rangle .\end{array}$$

(12)

By using these relations, we obtain the following microscopic mean field model:

$$\begin{array}{l}{m}_{i}^{\xi}(t+1)={g}_{{\beta}^{\xi}}\left[{\displaystyle \sum _{j\ne i}^{{N}_{\xi}}}({J}_{ij}^{\xi \xi}{A}_{j}^{\xi}(t))+{\displaystyle \sum _{j=1}^{{N}_{\eta}}}({J}_{ij}^{\xi \eta}{A}_{j}^{\eta}(t))+{I}^{\xi}\right],\end{array}$$

(13)

$$\begin{array}{l}{A}_{i}^{\xi}(t+1)={A}_{i}^{\xi}(t)-\frac{{A}_{i}^{\xi}(t)}{{\tau}_{\mathrm{\text{a}}}^{\xi}}+{m}_{i}^{\xi}(t){X}_{i}^{\xi}(t){U}_{i}^{\xi}(t)/{U}_{\mathrm{\text{se}}}^{\xi},\end{array}$$

(14)

$$\begin{array}{l}{X}_{i}^{\xi}(t+1)={X}_{i}^{\xi}(t)+\frac{1-{X}_{i}^{\xi}(t)}{{\tau}_{\mathrm{\text{R}}}^{\xi}}-{m}_{i}^{\xi}(t){X}_{i}^{\xi}(t){U}_{i}^{\xi}(t),\end{array}$$

(15)

$$\begin{array}{l}{U}_{i}^{\xi}(t+1)={U}_{i}^{\xi}(t)+\frac{{U}_{\mathrm{\text{se}}}^{\xi}-{U}_{i}^{\xi}(t)}{{\tau}_{\mathrm{\text{F}}}^{\xi}}+{U}_{\mathrm{\text{se}}}^{\xi}(1-{U}_{i}^{\xi}(t)){m}_{i}^{\xi}(t),\end{array}$$

(16)

where we have set ${m}_{i}^{\xi}(t)\equiv \langle {s}_{i}^{\xi}(t)\rangle $, ${A}_{i}^{\xi}(t)\equiv \langle {a}_{i}^{\xi}(t)\rangle $, ${X}_{i}^{\xi}(t)\equiv \langle {x}_{i}^{\xi}(t)\rangle $, and ${U}_{i}^{\xi}(t)\equiv \langle {u}_{i}^{\xi}(t)\rangle $, respectively.

We represent the fixed point of the microscopic mean field model as ${\overline{m}}_{i}^{\xi}$, ${\u0100}_{i}^{\xi}$, ${\overline{X}}_{i}^{\xi}$, and ${\u016a}_{i}^{\xi}$. The fixed point for Equations (14) to (16) is obtained as follows:

$$\begin{array}{l}{\u0100}_{i}^{\xi}=\frac{{\tau}_{\mathrm{\text{a}}}^{\xi}{\u016a}_{i}^{\xi}{\overline{m}}_{i}^{\xi}{\overline{X}}_{i}^{\xi}}{{U}_{\mathrm{\text{se}}}^{\xi}},\end{array}$$

(17)

$$\begin{array}{l}{\overline{X}}_{i}^{\xi}=\frac{1}{1+{\tau}_{\mathrm{\text{R}}}^{\xi}{\u016a}_{i}^{\xi}{\overline{m}}_{i}^{\xi}},\end{array}$$

(18)

$$\begin{array}{l}{\u016a}_{i}^{\xi}=\frac{{U}_{\mathrm{\text{se}}}^{\xi}(1+{\tau}_{\mathrm{\text{F}}}^{\xi}{\overline{m}}_{i}^{\xi})}{1+{\tau}_{\mathrm{\text{F}}}^{\xi}{U}_{\mathrm{\text{se}}}^{\xi}{\overline{m}}_{i}^{\xi}}.\end{array}$$

(19)

Using these equations, we obtain the value ${\overline{m}}_{i}^{\xi}$ at the fixed point as follows:

$$\begin{array}{l}{\overline{m}}_{i}^{\xi}={g}_{{\beta}^{\xi}}[{\displaystyle \sum _{j\text{}\ne \text{}i}^{{N}_{\xi}}{J}_{ij}^{\xi \xi}}\left(\frac{{\tau}_{\text{a}}^{\xi}{\overline{m}}_{j}^{\xi}(1+{\tau}_{\text{F}}^{\xi}{\overline{m}}_{j}^{\xi})}{1+({\tau}_{\text{F}}^{\xi}+{\tau}_{\text{R}}^{\xi}){U}_{\text{se}}^{\xi}{\overline{m}}_{j}^{\xi}+{U}_{\text{se}}^{\xi}{\tau}_{\text{F}}^{\xi}{\tau}_{\text{R}}^{\xi}{\overline{m}}_{j}^{\xi}{\overline{m}}_{j}^{\xi}}\right)\\ \text{}+\text{}{\displaystyle \sum _{j\text{}=\text{}1}^{{N}_{\eta}}{J}_{ij}^{\xi \eta}}\left(\frac{{\tau}_{\text{a}}^{\eta}{\overline{m}}_{j}^{\eta}(1+{\tau}_{\text{F}}^{\eta}{\overline{m}}_{j}^{\eta})}{1+({\tau}_{\text{F}}^{\eta}+{\tau}_{\text{R}}^{\eta}){U}_{\text{se}}^{\eta}{\overline{m}}_{j}^{\eta}+{U}_{\text{se}}^{\eta}{\tau}_{\text{F}}^{\eta}{\tau}_{\text{R}}^{\eta}{\overline{m}}_{j}^{\eta}{\overline{m}}_{j}^{\eta}}\right)+{I}^{\xi}].\end{array}$$

(20)

We derive a macroscopic mean field model by considering a network with all-to-all connections, where the weights ${J}_{ij}^{\xi \xi}$ and ${J}_{ij}^{\xi \eta}$ are given as follows:

$$\begin{array}{l}{J}_{ij}^{\xi \xi}=\frac{{J}_{0}^{\xi \xi}}{{N}_{\xi}},\end{array}$$

(21)

$$\begin{array}{l}{J}_{ij}^{\xi \eta}=\frac{{J}_{0}^{\xi \eta}}{{N}_{\eta}}.\end{array}$$

(22)

Here, ${J}_{0}^{\xi \xi}$ and ${J}_{0}^{\xi \eta}$ are the parameters specifying the strength of the uniform connections. Because of this synaptic connection uniformity, the variables ${m}_{i}^{\xi}$, ${A}_{i}^{\xi}$, ${X}_{i}^{\xi}$, and ${U}_{i}^{\xi}$ can be substituted with their respective population averages ${m}_{0}^{\xi}=(1/{N}_{\xi})\sum _{i=1}^{{N}_{\xi}}{m}_{i}^{\xi}$, ${A}_{0}^{\xi}=(1/{N}_{\xi})\sum _{i=1}^{{N}_{\xi}}{A}_{i}^{\xi}$, ${X}_{0}^{\xi}=(1/{N}_{\xi})\sum _{i=1}^{{N}_{\xi}}{X}_{i}^{\xi}$, and ${U}_{0}^{\xi}=(1/{N}_{\xi})\sum _{i=1}^{{N}_{\xi}}{U}_{i}^{\xi}$. The macroscopic mean field model for a network with uniform connections is given as follows:

$$\begin{array}{l}{m}_{0}^{\xi}(t+1)={F}_{m}^{\xi}(\Omega (t)),\end{array}$$

(23)

$$\begin{array}{l}{A}_{0}^{\xi}(t+1)={F}_{A}^{\xi}(\Omega (t)),\end{array}$$

(24)

$$\begin{array}{l}{X}_{0}^{\xi}(t+1)={F}_{X}^{\xi}(\Omega (t)),\end{array}$$

(25)

$$\begin{array}{l}{U}_{0}^{\xi}(t+1)={F}_{U}^{\xi}(\Omega (t)),\end{array}$$

(26)

where

$$\begin{array}{l}{F}_{m}^{\xi}(\Omega (t))={g}_{{\beta}^{\xi}}\left[{J}_{0}^{\xi \xi}{A}_{0}^{\xi}+{J}_{0}^{\xi \eta}{A}_{0}^{\eta}+{I}^{\xi}\right],\end{array}$$

(27)

$$\begin{array}{l}{F}_{A}^{\xi}(\Omega (t))={A}_{0}^{\xi}-\frac{{A}_{0}^{\xi}}{{\tau}_{\mathrm{\text{a}}}^{\xi}}+{m}_{0}^{\xi}{X}_{0}^{\xi}{U}_{0}^{\xi}/{U}_{\mathrm{\text{se}}}^{\xi},\end{array}$$

(28)

$$\begin{array}{l}{F}_{X}^{\xi}(\Omega (t))={X}_{0}^{\xi}+\frac{1-{X}_{0}^{\xi}}{{\tau}_{\mathrm{\text{R}}}^{\xi}}-{m}_{0}^{\xi}{X}_{0}^{\xi}{U}_{0}^{\xi},\end{array}$$

(29)

$$\begin{array}{l}{F}_{U}^{\xi}(\Omega (t))={U}_{0}^{\xi}+\frac{{U}_{\mathrm{\text{se}}}^{\xi}-{U}_{0}^{\xi}}{{\tau}_{\mathrm{\text{F}}}^{\xi}}+{U}_{\mathrm{\text{se}}}^{\xi}(1-{U}_{0}^{\xi}){m}_{0}^{\xi},\end{array}$$

(30)

with the state vector Ω(*t*) defined as follows:

$$\begin{array}{l}\Omega (t)={\left[{m}_{0}^{\mathrm{\text{E}}}(t),{m}_{0}^{\mathrm{\text{I}}}(t),{A}_{0}^{\mathrm{\text{E}}}(t),{A}_{0}^{\mathrm{\text{I}}}(t),{X}_{0}^{\mathrm{\text{E}}}(t),{X}_{0}^{\mathrm{\text{I}}}(t),{U}_{0}^{\mathrm{\text{E}}}(t),{U}_{0}^{\mathrm{\text{I}}}(t)\right]}^{\mathrm{\text{T}}}.\end{array}$$

(31)

As shown in the Results section, the dynamic mean field model is consistent with the simulation on the stochastic model.

By modifying Equation (20), the fixed point for the macroscopic mean field model can be calculated as follows:

$$\begin{array}{l}{\overline{m}}^{\mathrm{\text{E}}}={f}_{\mathrm{\text{E}}}({\overline{m}}^{\mathrm{\text{E}}},{\overline{m}}^{\mathrm{\text{I}}}),\end{array}$$

(32)

$$\begin{array}{l}{\overline{m}}^{\mathrm{\text{I}}}={f}_{\mathrm{\text{I}}}({\overline{m}}^{\mathrm{\text{I}}},{\overline{m}}^{\mathrm{\text{E}}}),\end{array}$$

(33)

where

$$\begin{array}{l}{f}_{\xi}({\overline{m}}^{\xi},{\overline{m}}^{\eta})\\ =\text{}{g}_{{\beta}^{\xi}}[{J}_{0}^{\xi \xi}\left(\frac{{\tau}_{\text{a}}^{\xi}{\overline{m}}^{\xi}\left(1+{\tau}_{\text{F}}^{\xi}{\overline{m}}^{\xi}\right)}{1+\left({\tau}_{\text{F}}^{\xi}+{\tau}_{\text{R}}^{\xi}\right){U}_{\text{se}}^{\xi}{\overline{m}}^{\xi}+{U}_{\text{se}}^{\xi}{\tau}_{\text{F}}^{\xi}{\tau}_{\text{R}}^{\xi}{\overline{m}}^{\xi}{\overline{m}}^{\xi}}\right)\\ +{J}_{0}^{\xi \eta}\left(\frac{{\tau}_{\text{a}}^{\eta}{\overline{m}}^{\eta}\left(1+{\tau}_{\text{F}}^{\eta}{\overline{m}}^{\eta}\right)}{1+\left({\tau}_{\text{F}}^{\eta}+{\tau}_{\text{R}}^{\eta}\right){U}_{\text{se}}^{\eta}{\overline{m}}^{\eta}+{U}_{\text{se}}^{\eta}{\tau}_{\text{F}}^{\eta}{\tau}_{\text{R}}^{\eta}{\overline{m}}^{\eta}{\overline{m}}^{\eta}}\right)+{I}^{\xi}].\end{array}$$

(34)

After solving the above equations simultaneously, we obtain the values ${\overline{m}}_{0}^{\mathrm{\text{E}}}$ and ${\overline{m}}_{0}^{\mathrm{\text{I}}}$ at the fixed point. By substituting ${\overline{m}}_{0}^{\mathrm{\text{E}}}$ and ${\overline{m}}_{0}^{\mathrm{\text{I}}}$ into the following fixed point equations,

$$\begin{array}{l}{\u0100}_{0}^{\xi}=\frac{{\tau}_{\mathrm{\text{a}}}^{\xi}{\u016a}_{0}^{\xi}{\overline{m}}_{0}^{\xi}{\overline{X}}_{0}^{\xi}}{{U}_{\mathrm{\text{se}}}^{\xi}},\end{array}$$

(35)

$$\begin{array}{l}{\overline{X}}_{0}^{\xi}=\frac{1}{1+{\tau}_{\mathrm{\text{R}}}^{\xi}{\u016a}_{0}^{\xi}{\overline{m}}_{0}^{\xi}},\end{array}$$

(36)

$$\begin{array}{l}{\u016a}_{0}^{\xi}=\frac{{U}_{\mathrm{\text{se}}}^{\xi}(1+{\tau}_{\mathrm{\text{F}}}^{\xi}{\overline{m}}_{0}^{\xi})}{1+{\tau}_{\mathrm{\text{F}}}^{\xi}{U}_{\mathrm{\text{se}}}^{\xi}{\overline{m}}_{0}^{\xi}},\end{array}$$

(37)

we obtain the values ${\u0100}_{0}^{\xi}$, ${\overline{X}}_{0}^{\xi}$, and ${\u016a}_{0}^{\xi}$ at the fixed point. Because nonlinear simultaneous equations cannot generally be solved analytically due to the function ${g}_{{\beta}^{\xi}}\left[\xb7\right]$, we use Newton's method to numerically obtain the fixed point in the Results section.

We analyze the stability of the fixed point with small deviations, $\delta {m}_{0}^{\xi}(t)$, $\delta {A}_{0}^{\xi}(t)$, $\delta {X}_{0}^{\xi}(t)$, and $\delta {U}_{0}^{\xi}(t)$, around the fixed point as follows:

$$\begin{array}{l}{m}_{0}^{\xi}(t)={\overline{m}}_{0}^{\xi}+\delta {m}_{0}^{\xi}(t),\end{array}$$

(38)

$$\begin{array}{l}{A}_{0}^{\xi}(t)={\u0100}_{0}^{\xi}+\delta {A}_{0}^{\xi}(t),\end{array}$$

(39)

$$\begin{array}{l}{X}_{0}^{\xi}(t)={\overline{X}}_{0}^{\xi}+\delta {X}_{0}^{\xi}(t),\end{array}$$

(40)

$$\begin{array}{l}{U}_{0}^{\xi}(t)={\u016a}_{0}^{\xi}+\delta {U}_{0}^{\xi}(t).\end{array}$$

(41)

By neglecting the higher order components, we obtain the following locally linearized equation:

$$\begin{array}{l}\left(\begin{array}{c}\delta {m}_{0}^{\mathrm{\text{E}}}(t+1)\\ \delta {A}_{0}^{\mathrm{\text{E}}}(t+1)\\ \delta {X}_{0}^{\mathrm{\text{E}}}(t+1)\\ \delta {U}_{0}^{\mathrm{\text{E}}}(t+1)\\ \delta {m}_{0}^{\mathrm{\text{I}}}(t+1)\\ \delta {A}_{0}^{\mathrm{\text{I}}}(t+1)\\ \delta {X}_{0}^{\mathrm{\text{I}}}(t+1)\\ \delta {U}_{0}^{\mathrm{\text{I}}}(t+1)\end{array}\right)=K\left(\begin{array}{c}\delta {m}_{0}^{\mathrm{\text{E}}}(t)\\ \delta {A}_{0}^{\mathrm{\text{E}}}(t)\\ \delta {X}_{0}^{\mathrm{\text{E}}}(t)\\ \delta {U}_{0}^{\mathrm{\text{E}}}(t)\\ \delta {m}_{0}^{\mathrm{\text{I}}}(t)\\ \delta {A}_{0}^{\mathrm{\text{I}}}(t)\\ \delta {X}_{0}^{\mathrm{\text{I}}}(t)\\ \delta {U}_{0}^{\mathrm{\text{I}}}(t)\end{array}\right),\end{array}$$

(42)

where *K* denotes the following Jacobian matrix:

$$\begin{array}{l}K=\left(\begin{array}{cccccccc}{K}_{mm}^{\mathrm{\text{EE}}}& {K}_{mA}^{\mathrm{\text{EE}}}& {K}_{mX}^{\mathrm{\text{EE}}}& {K}_{mU}^{\mathrm{\text{EE}}}& {K}_{mm}^{\mathrm{\text{EI}}}& {K}_{mA}^{\mathrm{\text{EI}}}& {K}_{mX}^{\mathrm{\text{EI}}}& {K}_{mU}^{\mathrm{\text{EI}}}\\ {K}_{Am}^{\mathrm{\text{EE}}}& {K}_{AA}^{\mathrm{\text{EE}}}& {K}_{AX}^{\mathrm{\text{EE}}}& {K}_{AU}^{\mathrm{\text{EE}}}& {K}_{Am}^{\mathrm{\text{EI}}}& {K}_{AA}^{\mathrm{\text{EI}}}& {K}_{AX}^{\mathrm{\text{EI}}}& {K}_{AU}^{\mathrm{\text{EI}}}\\ {K}_{Xm}^{\mathrm{\text{EE}}}& {K}_{XA}^{\mathrm{\text{EE}}}& {K}_{XX}^{\mathrm{\text{EE}}}& {K}_{XU}^{\mathrm{\text{EE}}}& {K}_{Xm}^{\mathrm{\text{EI}}}& {K}_{XA}^{\mathrm{\text{EI}}}& {K}_{XX}^{\mathrm{\text{EI}}}& {K}_{XU}^{\mathrm{\text{EI}}}\\ {K}_{Um}^{\mathrm{\text{EE}}}& {K}_{UA}^{\mathrm{\text{EE}}}& {K}_{UX}^{\mathrm{\text{EE}}}& {K}_{UU}^{\mathrm{\text{EE}}}& {K}_{Um}^{\mathrm{\text{EI}}}& {K}_{UA}^{\mathrm{\text{EI}}}& {K}_{UX}^{\mathrm{\text{EI}}}& {K}_{UU}^{\mathrm{\text{EI}}}\\ {K}_{mm}^{\mathrm{\text{IE}}}& {K}_{mA}^{\mathrm{\text{IE}}}& {K}_{mX}^{\mathrm{\text{IE}}}& {K}_{mU}^{\mathrm{\text{IE}}}& {K}_{mm}^{\mathrm{\text{II}}}& {K}_{mA}^{\mathrm{\text{II}}}& {K}_{mX}^{\mathrm{\text{II}}}& {K}_{mU}^{\mathrm{\text{II}}}\\ {K}_{Am}^{\mathrm{\text{IE}}}& {K}_{AA}^{\mathrm{\text{IE}}}& {K}_{AX}^{\mathrm{\text{IE}}}& {K}_{AU}^{\mathrm{\text{IE}}}& {K}_{Am}^{\mathrm{\text{II}}}& {K}_{AA}^{\mathrm{\text{II}}}& {K}_{AX}^{\mathrm{\text{II}}}& {K}_{AU}^{\mathrm{\text{II}}}\\ {K}_{Xm}^{\mathrm{\text{IE}}}& {K}_{XA}^{\mathrm{\text{IE}}}& {K}_{XX}^{\mathrm{\text{IE}}}& {K}_{XU}^{\mathrm{\text{IE}}}& {K}_{Xm}^{\mathrm{\text{II}}}& {K}_{XA}^{\mathrm{\text{II}}}& {K}_{XX}^{\mathrm{\text{II}}}& {K}_{XU}^{\mathrm{\text{II}}}\\ {K}_{Um}^{\mathrm{\text{IE}}}& {K}_{UA}^{\mathrm{\text{IE}}}& {K}_{UX}^{\mathrm{\text{IE}}}& {K}_{UU}^{\mathrm{\text{IE}}}& {K}_{Um}^{\mathrm{\text{II}}}& {K}_{UA}^{\mathrm{\text{II}}}& {K}_{UX}^{\mathrm{\text{II}}}& {K}_{UU}^{\mathrm{\text{II}}}\end{array}\right),\end{array}$$

(43)

and each element of this matrix is given as follows:

$$\begin{array}{l}{K}_{mA}^{\xi \xi}=\frac{\partial {F}_{m}^{\xi}}{\partial {A}^{\xi}}={g}_{{\beta}^{\xi}}^{\prime}\left[{h}^{\xi}\right]{J}_{0}^{\xi \xi},\end{array}$$

(44)

$$\begin{array}{l}{K}_{mA}^{\xi \eta}=\frac{\partial {F}_{m}^{\xi}}{\partial {A}^{\eta}}={g}_{{\beta}^{\xi}}^{\prime}\left[{h}^{\xi}\right]{J}_{0}^{\xi \eta},\end{array}$$

(45)

where

$$\begin{array}{l}{g}_{{\beta}^{\xi}}^{\prime}\left[{h}^{\xi}\right]=\frac{{\beta}^{\xi}}{2}\left[1-{tanh}^{2}({\beta}^{\xi}{h}^{\xi})\right],\end{array}$$

(46)

$$\begin{array}{l}{h}^{\xi}={J}_{0}^{\xi \xi}{A}^{\xi}+{J}_{0}^{\xi \eta}{A}^{\eta}+{I}^{\xi}.\end{array}$$

(47)

Furthermore, the remaining elements of matrix *K* are given as follows:

$$\begin{array}{l}{K}_{Am}^{\xi \xi}=\frac{{U}^{\xi}{X}^{\xi}}{{U}_{\mathrm{\text{se}}}^{\xi}},\end{array}$$

(48)

$$\begin{array}{l}{K}_{AX}^{\xi \xi}=\frac{{m}^{\xi}{U}^{\xi}}{{U}_{\mathrm{\text{se}}}^{\xi}},\end{array}$$

(49)

$$\begin{array}{l}{K}_{AU}^{\xi \xi}=\frac{{m}^{\xi}{X}^{\xi}}{{U}_{\mathrm{\text{se}}}^{\xi}},\end{array}$$

(50)

$$\begin{array}{l}{K}_{AA}^{\xi \xi}=1-\frac{1}{{\tau}_{\mathrm{\text{a}}}^{\xi}},\end{array}$$

(51)

$$\begin{array}{l}{K}_{Xm}^{\xi \xi}=-{U}^{\xi}{X}^{\xi},\end{array}$$

(52)

$$\begin{array}{l}{K}_{XX}^{\xi \xi}=(1-\frac{1}{{\tau}_{\mathrm{\text{R}}}^{\xi}})-{m}^{\xi}{U}^{\xi},\end{array}$$

(53)

$$\begin{array}{l}{K}_{XU}^{\xi \xi}=-{m}^{\xi}{X}^{\xi},\end{array}$$

(54)

$$\begin{array}{l}{K}_{Um}^{\xi \xi}={U}_{\mathrm{\text{se}}}^{\xi}(1-{U}^{\xi}),\end{array}$$

(55)

$$\begin{array}{l}{K}_{UU}^{\xi \xi}=(1-\frac{1}{{\tau}_{\mathrm{\text{F}}}^{\xi}})-{U}_{\mathrm{\text{se}}}^{\xi}{m}^{\xi},\end{array}$$

(56)

and other elements are zeroes. We numerically analyze the stability of the fixed point by eigenvalue analysis with this Jacobian matrix.

We use the slice analysis developed by Komuro et al. (2016), to elucidate the bifurcations of quasi-periodic oscillations arising from the proposed model. Mathematically, a slice with a width of ϵ is defined here as

$$\begin{array}{l}{\Sigma}^{\u03f5}=\left\{\Omega \in {\mathbb{R}}^{8}|\mathrm{\text{dist}}(\Omega ,\Sigma )<\u03f5\right\},\end{array}$$

(57)

where dist(·, ·) denotes the Euclidean distance between a point of the state vector Ω(*t*) (Equation 31) and a codimension-one plane Σ, called “section” (Komuro et al., 2016). Here, to differentiate a trajectory in the state space from one in the section qualitatively, we introduce the following two useful terms: a *d*-dimensional torus in map (MT*d*) and a *d*-dimensional torus in section (ST*d*) (Kamiyama et al., 2014; Komuro et al., 2016). Accordingly, for example, an MT1 (a closed curve) and an MT2 (a two-dimensional torus) are converted into an ST0 (an isolated point) and an ST1 (a closed curve), respectively, via the slice (Figure (Figure2).2). Because the bifurcation of MT*d* can be interpreted as that of ST(*d* − 1) almost equivalently, we apply the conventional bifurcation theory to ST*d* in order to consider its bifurcations as well as those of MT*d* (Komuro et al., 2016).

In this section, we analyze a variety of bifurcation structures arising from the proposed model by considering ${\tau}_{\mathrm{\text{a}}}^{\xi}$, ${J}_{0}^{\xi \xi}$, ${J}_{0}^{\xi \eta}$, *I*^{ξ}, and ${\tau}_{\mathrm{\text{R}}}^{\xi}/{\tau}_{\mathrm{\text{F}}}^{\xi}$ to be bifurcation parameters and *T*^{ξ} and ${U}_{\mathrm{\text{se}}}^{\xi}$ to be constants; we set *T*^{ξ} = 0.8 and ${U}_{\mathrm{\text{se}}}^{\xi}=0.1$, respectively, because the activation function ${g}_{{\beta}^{\xi}}\left[\xb7\right]$ used here is an idealized sigmoid form different from the experimentally known frequency-current relationship (Tateno et al., 2004) and ${U}_{\mathrm{\text{se}}}^{\xi}$ merely determines the steady value of the calcium concentration. The bifurcation analysis focuses on two types of dynamic synapses, namely short-term depression (${\tau}_{\mathrm{\text{R}}}^{\xi}/{\tau}_{\mathrm{\text{F}}}^{\xi}\gg 1$) and facilitation (${\tau}_{\mathrm{\text{R}}}^{\xi}/{\tau}_{\mathrm{\text{F}}}^{\xi}\ll 1$) synapses, where ${\tau}_{\mathrm{\text{R}}}^{\xi}$ has been fixed at 70 and ${\tau}_{\mathrm{\text{F}}}^{\xi}$ 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 ${J}_{0}^{\mathrm{\text{EI}}}=0$ and ${J}_{0}^{\mathrm{\text{EE}}}\ge 0$ for the excitatory network and ${J}_{0}^{\mathrm{\text{IE}}}=0$ and ${J}_{0}^{\mathrm{\text{II}}}<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}_{\xi}=1{0}^{4}$. 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 Figures33–**12**. First, we analyze the excitatory and inhibitory networks independently, so that a parameter space generating the oscillatory state becomes clear (Figures (Figures3,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 (Figures3,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 ${J}_{0}^{\mathrm{\text{EI}}}$ and ${J}_{0}^{\mathrm{\text{IE}}}$ (**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 (Figure44).

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*_{0}, *I*) phase diagram (Figure (Figure3A).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*_{0} decreases/increases or *I* increases from the region of the OSE state. On the other hand, when *J*_{0} 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*_{0} on (*I*, τ_{a}) = (−1, 2.5) (Figure (Figure3B),3B), shows the OSE state emergent from the steady state. The mean neural activity, *m*_{0}, in the steady state, increases with an increase in *J*_{0}, whereas the steady state destabilizes via the first NS bifurcation at *J*_{0} = 1.63 and then, the OSE state appears. Because of the second NS bifurcation at *J*_{0} = 3.48, the OSE state disappears, and accordingly, the steady state reappears. In contrast, the bifurcation diagram, with respect to *J*_{0} on (*I*, τ_{a}) = (1, 2.5) (Figure (Figure3C),3C), shows the OSI state emergent from the steady state. The mean neural activity, *m*_{0}, in the steady state decreases with a decrease in *J*_{0}, whereas the steady state destabilizes via the NS bifurcation at *J*_{0} = −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*_{0}, 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 (Figure3A).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 Figure3D3D 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 Figure3E3E 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 (Figure5).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 (Figure5A).5A). In contrast, τ_{R}/τ_{F} does not have much effect on the period for the OSI state (Figure (Figure5B5B).

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 (Figure6).6). The population averages of stochastic variables, ${s}_{0}(t)\left[=(1/N)\sum _{i}^{N}{s}_{i}(t)\right]$, ${a}_{0}(t)\left[=(1/N)\sum _{i}^{N}{a}_{i}(t)\right]$, ${x}_{0}(t)\left[=(1/N)\sum _{i}^{N}{x}_{i}(t)\right]$, and ${u}_{0}(t)\left[=(1/N)\sum _{i}^{N}{u}_{i}(t)\right]$, 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 Figure77 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: $({J}_{0}^{\mathrm{\text{EE}}},{I}^{\mathrm{\text{E}}},{\tau}_{\mathrm{\text{a}}}^{\mathrm{\text{E}}})=(2,-1,2.5)$ for the excitatory subnetwork, $({J}_{0}^{\mathrm{\text{II}}},{I}^{\mathrm{\text{I}}},{\tau}_{\mathrm{\text{a}}}^{\mathrm{\text{I}}})=(-10,1,12.5)$ for the inhibitory subnetwork, and ${\tau}_{\mathrm{\text{R}}}^{\xi}/{\tau}_{\mathrm{\text{F}}}^{\xi}=11.7$ for the depression synapses. This network with $({J}_{0}^{\mathrm{\text{EI}}},{J}_{0}^{\mathrm{\text{IE}}})=(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 $({J}_{0}^{\mathrm{\text{EI}}},{J}_{0}^{\mathrm{\text{IE}}})\ne (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 $({J}_{0}^{\mathrm{\text{EI}}},{J}_{0}^{\mathrm{\text{IE}}})$ phase diagram in Figure Figure7A).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 ${\tau}_{\mathrm{\text{R}}}^{\xi}/{\tau}_{\mathrm{\text{F}}}^{\xi}$, which merely changes the period of oscillations on subnetwork ξ.

Here we use the term MT*d* 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 ${J}_{0}^{\mathrm{\text{EI}}}$ on ${J}_{0}^{\mathrm{\text{IE}}}=2$, shows a transition between the OS1C and the OS2T states, where section ${A}_{0}^{\mathrm{\text{I}}}(t)=1.7$ and width ϵ = 0.001 have been used for the slice (Figure (Figure7B).7B). As ${J}_{0}^{\mathrm{\text{EI}}}$ increases from ${J}_{0}^{\mathrm{\text{EI}}}=-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 ${J}_{0}^{\mathrm{\text{EI}}}=-0.142$ for the OS1C state, and ${J}_{0}^{\mathrm{\text{EI}}}=-0.139$ for the OS2T state (Figure (Figure7F).7F). In contrast, as ${J}_{0}^{\mathrm{\text{EI}}}$ decreases from ${J}_{0}^{\mathrm{\text{EI}}}=-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 (Figure7D).7D). The negative exponents with multiplicity of two approach zero simultaneously as ${J}_{0}^{\mathrm{\text{EI}}}$ increases from ${J}_{0}^{\mathrm{\text{EI}}}=-0.145$; this is a property of the NS bifurcation. However, only one exponent starts to decrease as ${J}_{0}^{\mathrm{\text{EI}}}$ decreases from ${J}_{0}^{\mathrm{\text{EI}}}=-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 ${J}_{0}^{\mathrm{\text{IE}}}$ on ${J}_{0}^{\mathrm{\text{EI}}}=-0.05$, shows a transition between the OS2C and the OS2T states, where section ${A}_{0}^{\mathrm{\text{E}}}(t)=0.8$ and width ϵ = 0.001 have been used for the slice (Figure (Figure7C).7C). As ${J}_{0}^{\mathrm{\text{IE}}}$ decreases from ${J}_{0}^{\mathrm{\text{IE}}}=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 ${J}_{0}^{\mathrm{\text{IE}}}$ increases from ${J}_{0}^{\mathrm{\text{IE}}}=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 ${J}_{0}^{\mathrm{\text{IE}}}=3.12$ for the OS2C state and ${J}_{0}^{\mathrm{\text{IE}}}=3$ for the OS2T state (Figure (Figure7G).7G). The bifurcation involving this feature has been recently identified as the saddle-node cycle bifurcation of ST0 (ST0SNC), within the slice (Figure (Figure8A)8A) (Komuro et al., 2016), or as the MT1SNC bifurcation, outside the slice (Figure (Figure8B)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 (Figure7E);7E); only one negative exponent approaches zero as ${J}_{0}^{\mathrm{\text{IE}}}$ decreases from ${J}_{0}^{\mathrm{\text{IE}}}=9$.

Figures Figures99–**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, ${s}_{0}^{\xi}(t)\left[=(1/{N}_{\xi})\sum _{i}^{{N}_{\xi}}{s}_{i}^{\xi}(t)\right]$, ${a}_{0}^{\xi}(t)\left[=(1/{N}_{\xi})\sum _{i}^{{N}_{\xi}}{a}_{i}^{\xi}(t)\right]$, ${x}_{0}^{\xi}(t)\left[=(1/{N}_{\xi})\sum _{i}^{{N}_{\xi}}{x}_{i}^{\xi}(t)\right]$, and ${u}_{0}^{\xi}(t)\left[=(1/{N}_{\xi})\sum _{i}^{{N}_{\xi}}{u}_{i}^{\xi}(t)\right]$, 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 (Figures9C,9C, 10C, 11C); the trajectory in the state space (Figures (Figures9D,9D, 10D, 11D); and its distribution in the slice (Figures (Figures9E,9E, 10E, 11E).

Figure Figure99 shows the appearance of the OS1C state in the network dynamics. Excitatory neurons in this state fire less coherently, such that the excitatory subnetwork does not show slow oscillations, but exhibits fast fluctuations (Figure (Figure9A).9A). In contrast, inhibitory neurons fire coherently at a high frequency, such that the inhibitory subnetwork shows the fast oscillation (Figure (Figure9B).9B). The fast fluctuation/oscillation can be generated with external inhibitory input and with a connection from the inhibitory to excitatory subnetworks. These two types of inhibition are needed for the generation of the OS1C state. The fast oscillation in this state exhibits a single representative peak (Figure (Figure9C),9C), indicating that the amplitude of the oscillation is not modulated by the phase of the excitatory oscillation (Figures (Figures9B,9B, 12A). The dynamics emergent from both the excitatory and the inhibitory neurons exhibits a closed curve (Figure (Figure9D).9D). Accordingly, its slice includes less points, forming an ST0 (Figure (Figure9E9E).

Figure Figure1010 shows the appearance of the OS2T state in the network dynamics. Both excitatory and inhibitory neurons in this state fire coherently, such that all variables can show oscillations, whereas, the frequency differs between the excitatory and the inhibitory subnetworks. Excitatory oscillations are relatively slower than inhibitory ones (Figures 10A,B). The inhibitory fast oscillation exhibits two representative peaks (Figure 10C), indicating that the oscillation amplitude is slightly modulated by the phase of the excitatory slow oscillation (Figures 10B, 12B). This modulated oscillation is a phenomenon of phase-amplitude CFC and appears continuously (Figure 10B). Therefore, we call this as continuous CFC following Hyafil et al. (2015). The dynamics composed of both the slow and the fast oscillations exhibits a two-dimensional torus (Figure 10D). Accordingly, its slice includes a closed curve, forming an ST1 (Figure 10E).

Figure Figure1111 shows the appearance of the OS2C state in the network dynamics. The firing properties of excitatory and inhibitory neurons and their macroscopic oscillations in this state are similar to those in the OS2T state (Figures 11A–C), where the slow and the fast oscillations appear in the excitatory and the inhibitory subnetworks, respectively. However, the amplitude of the fast oscillation is more evidently modulated by the phase of the slow oscillation, compared to the OS2T state (Figures 11B, 12C). This modulated oscillation is likewise a phenomenon of phase-amplitude CFC but appears intermittently (Figure 11B). Therefore, we call this as intermittent CFC following Hyafil et al. (2015). The dynamics composed of both the slow and the fast oscillations exhibits a complicated closed curve (Figure 11D), whereas, its slice includes a few points, forming an ST0 (Figure 11E).

In this study, we analyzed a stochastic network model composed of excitatory and inhibitory neurons with dynamic synapses, and converted the model into the corresponding macroscopic mean field model. The bifurcation analysis of the mean field model revealed the overall dynamical properties of the network. The excitatory and the inhibitory subnetworks represent slow and fast oscillations, respectively. The interaction between these two subnetworks generates diverse oscillatory states with two major frequency components. This oscillatory phenomenon corresponds to phase-amplitude cross-frequency coupling (CFC). The dependence of the oscillatory states on coupling strengths, mediating between the subnetworks, has been clarified by the bifurcation analysis. Furthermore, it has been found that the oscillatory states of the CFC can be classified into two subtypes, namely, an oscillatory state with two frequency components on a two-dimensional torus (OS2T), which can generate the continuous CFC, and an oscillatory state with two frequency components on a closed curve (OS2C), which can generate the intermittent CFC.

The present model is an extension of an excitatory neural network model with dynamic synapses (Katori et al., 2012). The previous model, which corresponds to the excitatory subnetwork in this study, was modified in terms of the following three aspects: First, we analyzed the dependence of the network dynamics on the coupling strength, ${J}_{0}^{\xi \xi}$, and on the external input *I*^{ξ}; these parameters were fixed in the previous study (Katori et al., 2012). The analysis revealed that these two parameters can be crucial for the generation of a variety of oscillatory states. The second point is the introduction of an additional variable, ${a}_{i}^{\xi}(t)$, corresponding to the synaptic activity, and a parameter ${\tau}_{\mathrm{\text{a}}}^{\xi}$, corresponding to the time constant of the decay process for ${a}_{i}^{\xi}(t)$. The third point is a combination of the excitatory network with the inhibitory one, where parameters ${J}_{0}^{\mathrm{\text{EI}}}$ and ${J}_{0}^{\mathrm{\text{IE}}}$ were introduced.

Depending on the synaptic properties, the network dynamics changes. Decay time constants, ${\tau}_{\mathrm{\text{a}}}^{\mathrm{\text{E}}}$ and ${\tau}_{\mathrm{\text{a}}}^{\mathrm{\text{I}}}$ of the synaptic activity, may determine the frequency of slow oscillations in the excitatory network and that of fast oscillations in the inhibitory network, respectively (Figures 3D,E). The frequency can be likewise changed by ${\tau}_{\mathrm{\text{R}}}^{\xi}/{\tau}_{\mathrm{\text{F}}}^{\xi}$, the property of short-term plasticity (Figures 5A,B). A variation of frequency bands in neural activity, such as the delta, theta, alpha, beta, and gamma bands, is often observed in the brain, and it has been suggested that this variation correlates with the brain functions (Buzsáki and Draguhn, 2004). The brain functions may be attributed to the above synaptic parameters. Indeed, aminomethylphosphonic acid (AMPA) synapses have a relatively short time constant, whereas, N-methyl-D-aspartate (NMDA) synapses have a longer time constant. This synaptic difference should affect the generation of neural oscillations and brain functions.

We have found that the generation mechanism of the OSI state on the inhibitory subnetwork is qualitatively consistent with the physiological experiments (Fisahn et al., 1998; Mann et al., 2005). Inhibitory interneurons in the rat hippocampal CA3 region show a fast oscillation, which is referred to as the gamma oscillation. This oscillation is blocked by the AMPA or the gamma-aminobutyric acid (GABA) type-A receptor antagonist. AMPA-type synapses send excitatory input to the interneurons, while GABA type-A synapses send recurrent inhibitory connections. These antagonists can be considered as the realization of input *I*^{I} and the absolute value of coupling strength ${J}_{0}^{\mathrm{\text{II}}}$, respectively (Figure (Figure3A).3A). Taken together, the OSI state generated from the present model can be related to the gamma oscillation in inhibitory interneuron networks.

The network composed of both excitatory and inhibitory neurons shows phase-amplitude CFC, in which the amplitude of the fast oscillation is modulated by the phase of the slow oscillation (Figures 10B, 11B, 12B,C). The property of this oscillatory phenomenon is similar to that of the experimentally known CFC between the theta and the gamma oscillations observed in the entorhinal cortex of the hippocampus (Chrobak and Buzsáki, 1998).

Various oscillatory phenomena, including CFC arising from our model, may contribute to information coding in the brain. The presence of distinctive oscillatory states in the model implies that a variety of information coding schemes can exist in brain networks. It has been shown that the oscillatory states with two major frequency components can be classified into two subtypes, OS2T and OS2C states (Figures (Figures10,10, ,11).11). In the OS2T state, peaks of the fast oscillation are broadly distributed in the phase of the slow oscillation (Figure 10B). In contrast, in the OS2C state, the phase of the fast oscillation is partially locked by the slow oscillation; that is, peaks of the fast oscillation appear in specific phases of the slow oscillation (Figure 11B). The neural activity phase can be utilized to encode certain information. Indeed, it has been suggested that the physiologically observed CFC provides a basis for effective communications among neurons (Chrobak and Buzsáki, 1998).

The OS2C state may contribute to multi-item representation (Hyafil et al., 2015), because this state can generate the intermittent CFC in the inhibitory subnetwork. One cycle of the fast oscillation of the CFC would correspond to one item, associated with the working memory, the spatial memory, or the visual attention. Owing to intermittency of this oscillation, external multi-items could be effectively encoded in the brain; i.e., the slow oscillation, generated from the excitatory subnetwork, would play a role in optimizing storage capacity for the items. The typical fast oscillation, ${m}_{0}^{\mathrm{\text{I}}}(t)$, generating the CFC in the OS2C state depicted in Figure 11B, shows that approximately 11 items are possibly encoded. The number of items to be stored may increase or decrease depending on parameters, *I*^{E} and *I*^{I}, related to, e.g., visual input, because these parameters affect the frequency of the slow and the fast oscillations, respectively (Figures 3D,E). In contrast, the OS2T state would not be suitable for multi-item representation because this state generates continuous CFC. Moreover, the encoding scheme for the multi-items can be likewise observed on the closed curve underlying the fast oscillation ${m}_{0}^{\mathrm{\text{I}}}(t)$ in the OS2C state (Figure 11D); that is, the number of items to be stored would be limited in order for the brain to avoid producing wasted storage capacity. In contrast, the attractor underlying the OS2T state is a two-dimensional torus (Figure 10D); that is, more items would be encoded in the torus than the closed curve. However, the torus may not be efficient where only a few items are stored, because the torus consists of a dense orbit and may produce wasted storage capacity.

Our main finding is that the MT1SNC bifurcation may underlie a switching phenomenon between the continuous and the intermittent CFCs; this result supports the study of Fontolan et al. (2013). Fontolan et al. have reported that these two CFCs are switched via the saddle-node on invariant circle (SNIC) bifurcation, on a simplified Pyramidal Interneuron Network Gamma (PING) model (Fontolan et al., 2013). The SNIC bifurcation mediates between a limit cycle and a two-dimensional torus in flow, whereas the MT1SNC bifurcation mediates between a closed curve and a two-dimensional torus in map; here both bifurcations occur via a saddle-node cycle. If we assume that the MT1 arising from the proposed discrete-time model is the limit cycle in the corresponding continuous-time model, these two bifurcations will become consistent. Thus, the present study implies that phase-amplitude CFC in the brain can be interpreted in a discrete-time model.

Overall, the bifurcation analysis revealed that oscillatory dynamics, arising from the proposed model, qualitatively changes depending on parameters, which would be one origin of characterizing cell assemblies. Although the present study focused only on one cell assembly receiving inhibitory and external input, in fact, there exist many assemblies (Yoshimura et al., 2005), which could differ in their roles for neural information processing. Each of assemblies, in layer 2/3 of the cortex, selectively receives its recurrent connections and excitatory input from layer 4, possibly based on environmental change, while input from layer 5 might modulate activity between assemblies (Yoshimura et al., 2005). Such selectively interconnected neurons would play a crucial role for utilizing phase-amplitude CFC in the brain.

The mechanisms and functions of oscillatory phenomena must be further explored in the future. The oscillatory phenomena observed in the proposed model, a binary-state discrete-time neuron model, should be evaluated with a more realistic network model that reflects the detailed properties of the cerebral cortex.

YK and TS modeled neural networks with dynamic synapses. MK checked the details of the bifurcation analysis. KA supervised the overall research. TS analyzed the details of the model and did all other things.

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.

This research was supported by JSPS KAKENHI Grant Numbers 16K00246, 26280093, 15H05707, and CREST, JST.

- Aru J., Aru J., Priesemann V., Wibral M., Lana L., Pina G., et al. . (2015). Untangling cross-frequency coupling in neuroscience. Curr. Opin. Neurobiol. 31, 51–61. 10.1016/j.conb.2014.08.002 [PubMed] [Cross Ref]
- Benita J. M., Guillamon A., Deco G., Sanchez-Vives M. V. (2012). Synaptic depression and slow oscillatory activity in a biophysical network model of the cerebral cortex. Front. Comput. Neurosci. 6:64. 10.3389/fncom.2012.00064 [PMC free article] [PubMed] [Cross Ref]
- Bibitchkov D., Herrmann J. M., Geisel T. (2002). Pattern storage and processing in attractor networks with short-time synaptic dynamics. Netw. Comput. Neural Syst. 13, 115–129. 10.1080/net.13.1.115.129 [PubMed] [Cross Ref]
- Bonachela J. A., Franciscis S. D., Torres J. J., Munoz M. A. (2010). Self-organization without conservation: are neuronal avalanches generically critical? J. Stat. Mech. Theor. Exp. 6:P02015 10.1088/1742-5468/2010/02/p02015 [Cross Ref]
- Brunel N., Hakim V. (1999). Fast global oscillations in networks of integrate-and-fire neurons with low firing rates. Neural Comput. 11, 1621–1671. 10.1162/089976699300016179 [PubMed] [Cross Ref]
- Brunel N., Wang X. (2003). What determines the frequency of fast network oscillations with irregular neural discharges? I. Synaptic dynamics and excitation-inhibition balance. J. Neurophysiol. 90, 415–430. 10.1152/jn.01095.2002 [PubMed] [Cross Ref]
- Buzsáki G., Draguhn A. (2004). Neuronal oscillations in cortical networks. Science 304, 1926–1929. 10.1126/science.1099745 [PubMed] [Cross Ref]
- Canolty R. T., Edwards E., Dalal S. S., Soltani M., Nagarajan S. S., Kirsch H. E., et al. . (2006). High gamma power is phase-locked to theta oscillations in human neocortex. Science 313, 1626–1628. 10.1126/science.1128115 [PMC free article] [PubMed] [Cross Ref]
- Canolty R. T., Knight R. T. (2010). The functional role of cross-frequency coupling. Trends Cogn. Sci. 14, 506–515. 10.1016/j.tics.2010.09.001 [PMC free article] [PubMed] [Cross Ref]
- Chrobak J. J., Buzsáki G. (1998). Gamma oscillations in the entorhinal cortex of the freely behaving rat. J. Neurosci. 18, 388–398. [PubMed]
- Connors B. W., Gutnick M. J. (1990). Intrinsic firing patterns of diverse neocortical neurons. Trends Neurosci. 13, 99–104. 10.1016/0166-2236(90)90185-D [PubMed] [Cross Ref]
- Connors B. W., Gutnick M. J., Prince D. A. (1982). Electrophysiological properties of neocortical neurons
*in vitro*. J. Neurophysiol. 48, 1302–1320. [PubMed] - Cortes J. M., Torres J. J., Marro J., Garrido P. L., Kappen H. J. (2006). Effects of fast presynaptic noise in attractor neural networks. Neural Comput. 18, 614–633. 10.1162/neco.2006.18.3.614 [PubMed] [Cross Ref]
- Csicsvari J., Jamieson B., Wise K. D., Buzsáki G. (2003). Mechanisms of gamma oscillations in the hippocampus of the behaving rat. Neuron 37, 311–322. 10.1016/S0896-6273(02)01169-8 [PubMed] [Cross Ref]
- Engel A. K., Fries P., Singer W. (2001). Dynamic predictions: oscillations and synchrony in top-down processing. Nat. Rev. Neurosci. 2, 704–716. 10.1038/35094565 [PubMed] [Cross Ref]
- Fisahn A., Pike F. G., Buhl E. H., Paulsen O. (1998). Cholinergic induction of network oscillations at 40 Hz in the hippocampus
*in vitro*. Nature 394, 186–189. 10.1038/28179 [PubMed] [Cross Ref] - Fontolan L., Krupa M., Hyafil A., Gutkin B. (2013). Analytical insights on theta-gamma coupled neural oscillators. J. Math. Neurosci. 3:16. 10.1186/2190-8567-3-16 [PMC free article] [PubMed] [Cross Ref]
- Freeman W. J., Rogers L. J., Holmes M. D., Silbergeld D. L. (2000). Spatial spectral analysis of human electrocorticograms including the alpha and gamma bands. J. Neurosci. Methods 95, 111–121. 10.1016/S0165-0270(99)00160-0 [PubMed] [Cross Ref]
- Hyafil A., Giraud A. L., Fontolan L., Gutkin B. (2015). Neural cross-frequency coupling: connecting architectures, mechanisms, and functions. Trends Neurosci. 38, 725–740. 10.1016/j.tins.2015.09.001 [PubMed] [Cross Ref]
- Ibarz B., Casado J. M., Sanjuan M. A. F. (2011). Map-based models in neuronal dynamics. Phys. Rep. 501, 1–74. 10.1016/j.physrep.2010.12.003 [Cross Ref]
- Igarashi Y., Oizumi M., Okada M. (2010). Mean field analysis of stochastic neural network models with synaptic depression. J. Phys. Soc. Jpn. 79:84001 10.1143/JPSJ.79.084001 [Cross Ref]
- Jansen O., Colgin L. L. (2007). Cross-frequency coupling between neuronal oscillations. Trends Cogn. Sci. 11, 267–269. 10.1016/j.tics.2007.05.003 [PubMed] [Cross Ref]
- Kamiyama K., Komuro M., Endo T., Aihara K. (2014). Classification of bifurcations of quasi-periodic solutions using lyapunov bundles. Int. J. Bifurcation Chaos 24:1430034 10.1142/S0218127414300341 [Cross Ref]
- Katori Y., Igarashi Y., Okada M., Aihara K. (2012). Statility analysis of stochastic neural network with depression and facilitation synapses. J. Phys. Soc. Jpn. 81:114007 10.1143/JPSJ.81.114007 [Cross Ref]
- Katori Y., Otsubo Y., Okada M., Aihara K. (2013). Stability analysis of associative memory network composed of stochastic neurons and dynamic synapses. Front. Comput. Neurosci. 7:6. 10.3389/fncom.2013.00006 [PMC free article] [PubMed] [Cross Ref]
- Katori Y., Sakamoto K., Saito N., Tanji J., Mushiake H., Aihara K. (2011). Representational switching by dynamical reorganization of attractor structure in a network model of the prefrontal cortex. PLoS Comput. Biol. 7:e1002266. 10.1371/journal.pcbi.1002266 [PMC free article] [PubMed] [Cross Ref]
- Kawaguchi Y., Kubota Y. (1997). Gabaergic cell subtypes and their synaptic connections in rat frontal cortex. Cereb. Cortex 7, 476–486. 10.1093/cercor/7.6.476 [PubMed] [Cross Ref]
- Klimesch W. (1999). EEG alpha and theta oscillations reflect cognitive and memory performance: a review and analysis. Brain Res. Rev. 29, 169–195. 10.1016/S0165-0173(98)00056-3 [PubMed] [Cross Ref]
- Komuro M., Kamiyama K., Endo T., Aihara K. (2016). Quasi-periodic bifurcations of higher-dimensional tori. Int. J. Bifurcation Chaos 26:163001 10.1142/S0218127416300160 [Cross Ref]
- Kopell N., Ermentrout G. B., Whittington M., Traub R. D. (2000). Gamma rhythms and beta rhythms have different synchronization properties. Proc. Natl. Acad. Sci. U.S.A. 97, 1867–1872. 10.1073/pnas.97.4.1867 [PubMed] [Cross Ref]
- Kramer M. A., Roopun A. K., Carracedo L. M., Traub R. D., Whittington M. A., Kopell N. J. (2008). Rhythm generation through period concatenation in rat somatosensory cortex. PLoS Comput. Biol. 4:e1000169. 10.1371/journal.pcbi.1000169 [PMC free article] [PubMed] [Cross Ref]
- Kuznetsov Y. A. (2004). Elements of Applied Bifurcation Theory (Applied Mathematical Sciences Vol. 112), 3rd Edn. New York, NY: Springer-Verlag.
- Lux H. D., Pollen D. A. (1966). Electrical constants of neurons in the motor cortex of the cat. J. Neurophysiol 29, 207–220. [PubMed]
- Malerba P., Kopell N. J. (2012). Phase resetting reduces theta-gamma rhythmic interaction to a one-dimensional map. J. Math. Biol. 66, 1361–1366. 10.1007/s00285-012-0534-9 [PubMed] [Cross Ref]
- Mann E. O., Suckling J. M., Hajos N., Greenfield S. A., Paulsen O. (2005). Perisomatic feedback inhibition underlies cholinergically induced fast network oscillations in the rat hippocampus
*in vitro*. Neuron 45, 105–117. 10.1016/j.neuron.2004.12.016 [PubMed] [Cross Ref] - Markram H., Tsodyks M. (1996). Redistribution of synaptic efficacy between neocortical pyramidal neurons. Nature 382, 807–810. 10.1038/382807a0 [PubMed] [Cross Ref]
- Markram H., Wang Y., Tsodyks M. (1998). Differential signaling via the same axon of neocortical pyramidal neurons. Proc. Natl. Acad. Sci. U.S.A. 95, 5323–5328. 10.1073/pnas.95.9.5323 [PubMed] [Cross Ref]
- Marro J., Torres J. J., Cortes J. M. (2007). Chaotic hopping between attractors in neural networks. Neural Netw. 20, 230–235. 10.1016/j.neunet.2006.11.005 [PubMed] [Cross Ref]
- Matsumoto N., Ide D., Watanabe M., Okada M. (2007). Retrieval property of attractor network with synaptic depression. J. Phys. Soc. Jpn. 76:084005 10.1143/JPSJ.76.084005 [Cross Ref]
- Mejias J. F., Hernandez-Gomez B., Torres J. J. (2012). Short-term synaptic facilitation improves information retrieval in noisy neural networks. EPL 97:48008 10.1209/0295-5075/97/48008 [Cross Ref]
- Mejias J. F., Kappen H. J., Torres J. J. (2010). Irregular dynamics in up and down cortical states. PLoS ONE 5:e13651. 10.1371/journal.pone.0013651 [PMC free article] [PubMed] [Cross Ref]
- Mejias J. F., Torres J. J. (2008). The role of synaptic facilitation in spike coincidence detection. J. Comput. Neurosci. 24, 222–234. 10.1007/s10827-007-0052-8 [PubMed] [Cross Ref]
- Mejias J. F., Torres J. J. (2009). Maximum memory capacity on neural networks with short-term synaptic depression and facilitation. Neural Comput. 21, 851–871. 10.1162/neco.2008.02-08-719 [PubMed] [Cross Ref]
- Mejias J. F., Torres J. J. (2011). Emergence of resonances in neural systems: the interplay between adaptive threshold and short-term synaptic plasticity. PLoS ONE 6:e17255. 10.1371/journal.pone.0017255 [PMC free article] [PubMed] [Cross Ref]
- Mongillo G., Barak O., Tsodyks M. (2008). Synaptic theory of working memory. Science 319, 1543–1546. 10.1126/science.1150769 [PubMed] [Cross Ref]
- Nowak L. G., Azouz R., Sanchez-Vives M. V., Gray C. M., McCormick D. A. (2003). Electrophysiological classes of cat primary visual cortical neurons
*in vivo*as revealed by quantitative analyses. J. Neurophysiol. 89, 1541–1566. 10.1152/jn.00580.2002 [PubMed] [Cross Ref] - Otsubo Y., Nagata K., Oizumi M., Okada M. (2011). Influence of synaptic depression on memory storage capacity. J. Phys. Soc. Jpn. 80:084004 10.1143/JPSJ.80.084004 [Cross Ref]
- Pantic L., Torres J. J., Kappen H. J. (2003). Coincidence detection with dynamic synapses. Network 14, 17–33. 10.1088/0954-898X/14/1/302 [PubMed] [Cross Ref]
- Pantic L., Torres J. J., Kappen H. J., Gielen S. C. (2002). Associative memory with dynamic synapses. Neural Comput. 14, 2903–2923. 10.1162/089976602760805331 [PubMed] [Cross Ref]
- Pinamonti G., Joaquin J. M., Torres J. (2012). Stochastic resonance crossovers in complex networks. PLoS ONE 7:e51170. 10.1371/journal.pone.0051170 [PMC free article] [PubMed] [Cross Ref]
- Roopun A. K., Kramer M. A., Carracedo L. M., Kaiser M., Davies C. H., Traub R. D., et al. . (2008). Period concatenation underlies interactions between gamma and beta rhythms in neocortex. Front. Cell. Neurosci. 2:1. 10.3389/neuro.03.001.2008 [PMC free article] [PubMed] [Cross Ref]
- Rulkov N. F. (2002). Modeling of spiking-bursting neural behavior using two-dimensional map. PRE 65:041922. 10.1103/PhysRevE.65.041922 [PubMed] [Cross Ref]
- Sirota A., Csicsvari J., Buhl D., Buzsáki G. (2003). Communication between neocortex and hippocampus during sleep in rodents. Proc. Natl. Acad. Sci. U.S.A. 100, 2065–2069. 10.1073/pnas.0437938100 [PubMed] [Cross Ref]
- Steriade M. (2001). Impact of network activities on neuronal properties in corticothalamic systems. J. Neurophysiol. 86, 1–39. [PubMed]
- Tateno T., Harsch A., Robinson H. P. C. (2004). Threshold firing frequency-current relationship of neurons in rat somatosensory cortex: type 1 and type 2 dynamics. J. Neurophysiol. 92, 2283–2294. 10.1152/jn.00109.2004 [PubMed] [Cross Ref]
- Thomson A. M. (2000). Facilitation, augmentation and potentiation at central synapses. Trends Neurosci. 23, 305–312. 10.1016/S0166-2236(00)01580-0 [PubMed] [Cross Ref]
- Torres J. J., Cortes J. M., Marro J., Kappen H. J. (2007). Competition between synaptic depression and facilitation in attractor neural networks. Neural Comput. 19, 2739–2755. 10.1162/neco.2007.19.10.2739 [PubMed] [Cross Ref]
- Torres J. J., Marro J. (2015). Brain performance versus phase transitions. Sci. Rep. 5:12216. 10.1038/srep12216 [PMC free article] [PubMed] [Cross Ref]
- Torres J. J., Marro J., Mejias J. F. (2011). Can intrinsic noise induce various resonant peaks? New J. Phys. 13:053014 10.1088/1367-2630/13/5/053014 [Cross Ref]
- Torres J. J., Pantic L., Kappen H. J. (2002). Storage capacity of attractor neural networks with depressing synapses. Phys. Rev. E 66:06190. 10.1103/PhysRevE.66.061910 [PubMed] [Cross Ref]
- Tsodyks M., Pawelzik K., Markram H. (1998). Neural networks with dynamic synapses. Neural Comput. 10, 821–835. 10.1162/089976698300017502 [PubMed] [Cross Ref]
- Wang Y., Markram H., Goodman P. H., Berger T. K., Ma J., Goldman-Rakic P. S. (2006). Heterogeneity in the pyramidal network of the medial prefrontal cortex. Nat. Neurosci. 9, 534–542. 10.1038/nn1670 [PubMed] [Cross Ref]
- White J. A., Banks M. I., Pearce R. A., Kopell N. J. (2000). Networks of interneurons with fast and slow-aminobutyric acid type A (GABA
_{A}) kinetics provide substrate for mixed gamma-theta rhythm. Proc. Natl. Acad. Sci. U.S.A. 97, 8128–8133. 10.1073/pnas.100124097 [PubMed] [Cross Ref] - Yoshimura Y., Dantzker J. L. M., Callaway E. M. (2005). Excitatory cortical neurons form fine-scale functional networks. Nature 433, 868–873. 10.1038/nature03252 [PubMed] [Cross Ref]

Articles from Frontiers in Computational Neuroscience are provided here courtesy of **Frontiers Media SA**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |