Bifurcations of Emergent Bursting in a Neuronal Network

Complex neuronal networks are an important tool to help explain paradoxical phenomena observed in biological recordings. Here we present a general approach to mathematically tackle a complex neuronal network so that we can fully understand the underlying mechanisms. Using a previously developed network model of the milk-ejection reflex in oxytocin cells, we show how we can reduce a complex model with many variables and complex network topologies to a tractable model with two variables, while retaining all key qualitative features of the original model. The approach enables us to uncover how emergent synchronous bursting can arise from a neuronal network which embodies known biological features. Surprisingly, the bursting mechanisms are similar to those found in other systems reported in the literature, and illustrate a generic way to exhibit emergent and multiple time scale oscillations at the membrane potential level and the firing rate level.


Introduction
In neural systems, oscillatory rhythms have essential roles in sensory, cognitive, and motor functioning; in many experimental conditions [1][2][3], diverse physiological information can be encoded by the oscillatory activity of neuronal ensembles. However, the mechanisms by which rhythmic dynamics are produced vary considerably, from single pacemaker neurons, which can be mathematically described by voltage threshold models such as the integrate-and-fire model [4,5], or the more biophysical Hodgkin-Huxley type model [6], to large cortical networks, where interactions between neurons are responsible for the rhythmic behaviors (see [7][8][9] and the references therein).
Thus rhythmic oscillations can be observed and studied at different levels in neural systems, from the single neuron level, to the neuronal population level. Synchronous spikes in a neuronal population, which is a special case of population oscillating dynamics, may play an essential role in neuronal computation in cognition [14], and attention selection [15][16][17][18]. Synchronization is a population behavior, and accordingly has to be studied at the network level, and as shown in [19,20], synaptic interactions can be one cause of synchronous dynamics. Synchronous bursting emerges periodically in neuronal networks at a time scale of minutes, much longer than the millisecond time scale of individual neuronal spikes. Synchronous behavior can also be characterized as metastability, i.e. a transmission between different patterns [21,22], rather than attractors.
Some neuronal networks can exhibit rhythmic oscillations at multiple time scales. An interesting example is reported in a recent paper [23], in which a neuronal network model was developed to reproduce paradoxical phenomena observed from recordings of oxytocin-secreting neurons. Oxytocin is a hormone that is released by neuroendocrine neurons into the blood where it can trigger milk let-down in lactation, and it is also released within the brain, where it has powerful behavioral effects. Notably, in humans it is reported that oxytocin can increase the bonding and trust between individuals. These effects have made oxytocin a key drug target for new therapies aimed at mental disorders of social behavior such as autism.
The oxytocin network model in [23] was developed to explain the observed activity of oxytocin neurons in response to suckling. When young suckle, they are rewarded intermittently with a letdown of milk that results from reflex secretion of oxytocin; without oxytocin, newly born young will die unless they are fostered [24].
Oxytocin is made by magnocellular hypothalamic neurons, and is secreted from their nerve endings in the pituitary in response to action potentials (spikes) that are generated in the cell bodies and which are propagated down their axons to the nerve endings. Normally, oxytocin cells discharge asynchronously at 1-3 spikes/s, but during suckling, every 5 min or so, each discharges a brief, intense burst of spikes that release a pulse of oxytocin into the circulation [23]. The near-synchronous bursting is the consequence of vesicles of oxytocin released from the dendrites of oxytocin neurons as a result of spike activity, and this release of oxytocin can activate other oxytocin neurons via its effects on neighboring dendrites. The model revealed how emergent synchronous bursting at a very low frequency could arise from a neuronal network which implements all known features of the physiology of oxytocin cells. In that model, bursting is an emergent behavior of a complex system, involving both positive and negative feedbacks, between many sparsely connected cells. The oxytocin cells are regulated by independent random afferent inputs, but they are also excited by the dendritic release of oxytocin and inhibited by endocannabinoids, which are also produced by oxytocin neurons as a result of spike activity. The oxytocin that is released from the dendrites does not only have a local role; so much is released that it can act at distant sites, where it is believed to mediate one of the benefits of breast feeding: increasing the bonding between mother and baby.
A simple version of the network model is illustrated in Fig. 1. This model network has 48 cells and 12 bundles, and each cell has two dendrites ended up in different bundles, and two cells can interact if they share a common bundle. Each bundle contains the same number of dendrites, which we refer to as a 'homogeneous arrangement of the connections' (Fig. 1A, B). In the model, the dendritic stores of readily-releasable vesicles are continuously incremented by the suckling-related 'priming' input. Their level increases relatively steadily between bursts despite activitydependent depletion, and synchronous bursts tend to occur when the oxytocin level at the store is relatively high (Fig. 1C, D). In addition to the synchronicity, the bursts possess the characteristic that the inter-burst intervals are almost constant. More interestingly, we observed a number of paradoxical behaviors also observed in experimental studies. For example, increased spike activity between bursts enhances depletion of the stores and so can delay or even suppress bursting (Fig. 1E). Conversely, an increase in inhibitory inputs can promote the reflex in a system which fails to express bursting because of insufficient priming. For example, injections of the inhibitory neurotransmitter GABA into the supraoptic nucleus of a suckled, lactating rat can trigger milkejection bursts (Fig. 1F).
This neuronal network illustrates a hierarchical rhythmic oscillation dynamics. Each neuron discharges spikes periodically in a way that can be regarded as oscillating dynamics at the neuron level (the msec time scale for the inter-spike-interval); the network population synchronizes and exhibits bursting dynamics periodically in a way that can be regarded as oscillating dynamics at the network level (the minute time scale of the inter-burst intervals), comparing Fig. 1C with Fig. 1D. In general, a network system can have diverse oscillation dynamics at different levels, because of the interactions between individual units. Each node oscillates and exhibits a faster rhythmic dynamics, but the network synchronization also oscillates, with slower, rhythmic dynamics.
Different approaches have been proposed to deal with hierarchical rhythmic dynamics, as exemplified by neuronal bursting. The theory of slow-fast dynamical systems was introduced to explain how a neuron model can demonstrate coexistence of tonic spiking and bursting [25,26]. Abundant bifurcation behaviors in oscillations including spiking and bursting were detected in various neuron models [27][28][29][30][31][32], and are thought to be biophysically plausible. Moreover, reduction of complex neuronal networks to models with a few variables was performed, and mean field models were constructed to describe the average activity of the neuron systems [26,[33][34][35].
In this study, we aim to explain why and how emergent bursting occurs in the oxytocin network, and to reveal the underlying mechanisms of a particular puzzle: how increasing excitatory inputs can sometime stop the burst and increasing inhibitory inputs can promote the burst. Despite the many published papers in this area, we find that a novel approach is required. Most theories only deal with deterministic dynamics, but in the more biologically realistic oxytocin model, each neuron receives stochastic (Poisson) inputs, so an approximation to simplify each single neuron model is needed [36]. We approximate the system by a two-dimensional slow-fast dynamical system, where the variables used are threshold and oxytocin store level. This simplification is achieved after elaborately testing different model variables. The original model included many variables that were needed to match the physiological data quantitatively, including a hyperpolarizing after potential (HAP) and a slow afterhyperpolarising potential (AHP), different delays in the systems, and variables to model endocannabinoids actions. This complexity makes the original model hard to deal with mathematically. After eliminating non-essential variables, we conclude that a model incorporating just the dynamics of the threshold and oxytocin store level can be used to mimic the original model. Using this two-dimensional model, we can then apply bifurcation theory to explore the hierarchical rhythmic dynamics. We find there exists a critical value of the input rate beyond which bursting can emerge. This phenomenon can be described by a saddle-node bifurcation of limit cycles. As excitatory inputs increase in frequency, synchronized bursts arise in such a manner that the intervals between bursts are constant. More interestingly, and counter intuitively, the bursts disappear when the excitatory input frequency passes a larger critical value corresponding to another saddle-node bifurcation of limit cycles. We also detect occurrences of the subcritical Hopf bifurcation as the input frequency varies between the above two critical values. The saddle-node bifurcation plays a more significant role corresponding to the generation and ending of the bursting activity in the network.
Neuronal oscillation is a concept with many facets. The oxytocin model is a typical example among them and has its own specificities in comparison with other oscillating systems such as cortical oscillations relevant for cognition, as reviewed above. The response of the oxytocin network to an external stimulus is much slower, in comparison with cortical oscillations and therefore it involves very different intracellular and extracellular mechanisms. We emphasize here that the underlying chemical and physical mechanisms leading the oxytocin network to oscillate have very little to do with mechanisms for the development of cortical oscillations relevant for cognitions. For example, slow oscillations such as the theta-rhythm could arise from the GABAslow current, as recently modeled in the decision making [5] and theta-nested gamma oscillations [37]. Nevertheless, it is interesting to note that the mathematical approach developed here could be useful and serve as a general purpose tool to tackle a model, no matter it is a simplified integrate-and-fire model network or a biophysically realistic Hodgkin-Huxley type model network. Finally, synchronization in the oxytocin network is bursting-tobusting and such a mechanism of population spikes in neocortical networks for cognition has long been postulated in the literature [38].

Oxytocin Neuronal Network
In [23], a neuronal network model, based on leaky integrate-fire neurons with adaptive thresholds, dependent on the store level of oxytocin, was proposed and shown to exhibit emergent bursting dynamics. In this model, neurons receive random synaptic inputs from other, external neurons (excitatory and inhibitory post synaptic potentials, EPSPs and IPSPs). Changes in excitability of the neurons were modeled as changes in the membrane potential threshold for triggering spikes, and depend on previous spike activity and on dendritic oxytocin release, which is non-linearly related to spike activity and proportional to the size of the readilyreleasable store of oxytocin in the dendrites. As a closed loop, the stores of oxytocin that are available for release decrease when oxytocin is released from the dendrites but are increased as a result of the suckling stimulus. In the current paper we consider a simplified version of this model, which still preserves the synchronous bursting behavior. We refer the model in [23] as the original neuronal network (ONN) and our simplified model as the simplified neuronal network (SNN).
The core step of our simplification is the topology of the network. We consider a neuronal network with N neurons and n b bundles, where each neuron has two dendrites in different bundles. We assume that the network is homogeneously arranged, i.e. each of the n b bundles contains the same number of dendrites. In [23], we modeled the individual oxytocin neurons using the leaky integrate-and-fire model, modified to incorporate activity-dependent changes in excitability. The membrane potential v i of cell i obeys where t is the membrane time constant, v rest is the resting potential, N j E,i ,N j I,i are independent Poisson processes with the varied excitatory input ratel j E,i and the fixed inhibitory input rate are the magnitude of single EPSPs and IPSPs at v rest , and v E , v I are the excitatory and inhibitory reversal potentials. A spike is produced in cell i at time t~t s i , s~1,2,:: is the spike threshold at time t. After a spike, v i is reset to v rest . Activity-dependent changes in excitability and the effects of oxytocin are modeled by effects on spike threshold. Different from the model for the dynamical threshold in [23], we eliminate the effects of HAP and AHP in the spike threshold, so, Where T 0 is a constant. The increase in excitability due to oxytocin is modeled by T OT , Where t OT , k OT are constants, r m j (t) is the instantaneous release rate from dendrite m of cell j, and the sums pick up all the cells whose dendrites share the same bundle as cell i. The network topology is represented by matrices C k~f c k ij g, k~1, :::, n b ; c k ij~1 if dendrite j of cell i is in bundle k, and zero otherwise.
The readily-releasable store of oxytocin in dendrite j of cell i is represented by r j i , where Where t r is a time constant, k p is the rate of priming due to the suckling input, and r j i is the instantaneous release rate from dendrites j. In [23], the release of oxytocin is proportional to the readily-releasable stores: Where k r is the maximum fraction of the stores that can be released by a spike, D is a fixed delay before release, and the summation extends over the set ft s i vt,t s i {t s{1 i vt rel g, with t rel a constant. This ensures that only spikes occurring at intervals of less than t rel ('doublet' spikes) induce any release from dendrites. Here, we neglect the delay term D in (4) and the doublet effects by letting t rel~z ?, which means that spikes occurring at intervals of any length can induce release.
The model in [23] also took the inhibitory effects of endocannabinoids into consideration, but here we neglect it for simplicity.
The parameter values for simulations are as in Table 1. The ONN in [23] displays the transition between spiking and bursting (Fig. 2). The spiking rate is recorded on a network of 48 neurons and 12 bundles in Fig. 2A, and the voltage trace and store level of oxytocin are shown in Fig. 2C and E. The bursting events are essentially attributed to the drop of the spike threshold (red line) and store level. Our simplification of the ONN does not destroy such basic behaviors of the network in the sense that the SNN displays similar network activity in Fig. 2B, 2D and 2F as the ONN in Fig. 2A, 2C and 2E. As expected, the SNN fires faster than the ONN even though the input rate l j E,i in the SNN (50 Hz) is smaller than in the ONN (80 Hz), because we have discarded all bursting terminating mechanisms related to the negative feedback effects of the HAP and AHP on the spike threshold, the doublet effects in the impulsive release of oxytocin and the feedback inhibition by endocannabinoids.
Next we regard t s i as a series of random variables, and use Brownian motion to approximate the discrete spiking series, resulting in the following approximated release rate: activity of all 48 cells in the network through the first simulated milk-ejection burst. Where m output,i is the spiking rate and s output,i is the variance of the correlated Brownian motions B i (t),i~1,:::,N.
Because of the assumption that the network is homogeneously arranged and the observation that the neuronal population is activated synchronously, we can make a useful approximation by employing the mean field method. Explicitly, let T OT ,r, and m output denote the corresponding dynamical variables averaged over the entire population, and suppose that the number of entities in the summation in (2) is n(n~4N=n b ). As a first approximation, we can ignore the random effect and then omit all the subscripts in T OT, i , r j i , m output, i in (2), (3) and (5). A two-dimensional determinant dynamical system that describes the behavior of the averaged neuronal activity is as follows: We make a further simplification by removing the limit on the maximal value of the reduction of the spike threshold, which is set to be 25 mV in [23].

Firing Rate Map Approximation
In the system (6), m output (t) is an unknown term varying with time, which makes (6) a non-autonomous system. To overcome this difficulty, we present a method to evaluate the mean firing rate m output (t) of the network activity so that the system (6) becomes a mathematically tractable autonomous system. Intuitively, the firing rate m output (t) varies in response to the fluctuation of the spike threshold T(t) and the frequency of the afferent input l E . If we write m output (t) as a function of T and l E : a firing rate map, and substitute (7) in (6), we obtain the following two-dimensional system with two parameters l E ,n: To find the analytical expression of the firing rate map, we adopt a numerical approach by simulating the leaky integrate-fire model. Simulations of equation (1) for a single cell are conducted by fixing T on each trial. Fig. 3A shows the relationship between m output and T corresponding to varied excitatory inputs. Note that there are three discontinuities in the output rate for varied threshold shown in Fig. 3A, resulting from the altered mechanism of triggering a spike. The reason can be found in the parameter setting of the integrate-and-fire model: (1) and Table 1), which coincides with the gap between two consecutive discontinuities. To illustrate this, consider the discontinuity at T~{58mV. Suppose the neuron is initially in a resting state (v~v rest~{ 62mV) and the threshold T is fixed below {58mV. A single excitatory input to a neuron will result in a membrane potential increment of 4 mV, triggering a spike. In such a case, the neuron partially loses the dynamical behavior modeled by (1), which indicates distinct mechanism from the case that the threshold is fixed above {58mV. Similar explanations can be given for the other two discontinuities.
It seems straightforward to use the discontinuous firing rate map in the mean field model (8), but when we simulate the trajectories of (8) (see the results section), we face the challenge that the simulation is either computationally expensive if the output rate is derived from each desired T and l E in Figure 3A, or far from accurate, especially when the parameter is near the bifurcation point if the T-axis is partitioned and the output rate is derived from a neighboring point of the desired threshold. Furthermore, a discontinuous vector field is intractable in the bifurcation analysis (see the results section). Therefore, we need a continuous surrogate for the discontinuous version of the firing rate map.
Given the shape of the firing rate map in Fig. 3A, we use a sigmoid-like function to fit it: Here a(l E ) is the center of the curve, and b(l E ) is a tunable factor that controls the sharpness, c(l E ) is the term to describe the spike activity when the spike threshold is at the initial level T 0~{ 50Hz. By numerical experiments, we find a(l E ){ 66z0:02l E , b(l E )~ffi ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 0:02(l E z20) p and c(l E )~35 l E 200 To summarize all procedures above, we include a flow chart (Fig. 4). In the first step, we simplify a network model with a single neuron of 10 variables by discarding the negative feedbacks in the spike threshold and the doublet effects on the impulsive release of oxytocin, and obtain a simplified network model with four variables for each neuron. After evaluating the firing rate map, we derive the reduced deterministic autonomous system (8) (the mean field model) in the second step, which enables us to perform the bifurcation analysis. A similar approach could be employed generally to deal with other complex and stochastic neuronal networks.

Bifurcation Analysis
The value of the spike threshold T is closely related to the appearance of bursting behavior. In particular, a lower value of Tcan trigger a burst. Therefore, a systematic exploration of the dynamical properties of the system (8) enables us to understand the mechanism of the entire network activities.
Allowing l E (resp. n) to vary while keeping n (resp. l E ) fixed, the system (8) displays two types of bifurcations: the saddle-node bifurcation of limit cycles and the subcritical Hopf bifurcation. To exemplify this conclusion of the bifurcations in Eq. (8), we first investigate the system's dynamical behavior by fixing n~n 0~2 2 and varying l E .
When l E is small, the unique fixed point equilibrium in the r-T OT plane is asymptotically stable. Thus, from the asymptotical convergence of the trajectory if ½r(t),T OT (t) (Fig. 5A) we conclude that there is no bursting activity.
When l E increasingly exceeds a critical value l c E &60:1386343160437030, the saddle-node bifurcation of limit cycles occurs. To demonstrate the existence of this bifurcation and verify the stability of the bifurcated limit cycle, we construct the Poincaré map of (8)| lE . Denote the equilibrium point by x 0 (l E )~½r 0 (l E ),T OT,0 (l E ) and set L l E~f ½r 0 (l E ),T OT D T OT ƒT OT,0 (l E ), T OT [Rg, so that L l E is a half line transversal to the vector field in the neighborhood of the equilibrium x 0 (l E ). Here we introduce a new coordinate system along L l E , where x 0 (l E ) is an origin and n l E is a unit vector parallel to L l E . Hence, a becomes the coordinate of a point x on L l E if x~an l E zx 0 (l E ) for some a §0. Now, suppose that w l E (t,x) represents the solution of (8)| l E with the initial point x. Mathematically, it can be validated that there exists a number t t(x)). In other words, w l E ( t t(x),x) is the point at which the trajectory intersects with L l E for the first time after it departures from the initial point x. Thus, the coordinate a a of w l E ( t t(x),x) can be uniquely determined through w l E ( t t(x),x)~ a an l E zx 0 (l E ), and consequently the Poincaré map, denoted by P : R z ?R z , is established by P(a)~ a a for aw0. Fig. 6B shows the curves of the constructed Poincaré map Pfor different values of l E , where, clearly, each intersection between the curves and the black line P 1 (a)~a is a fixed point of P. When l E is smaller, P has no fixed point for a[R z . When l E~l c E , it has a unique fixed point. Since the quantity 1{DP'(a)D at the two sides of the fixed point has different signs, this fixed point is attracting on the right side and repelling on the left. Whenl E becomes slightly larger than l c E , two fixed points branch off: one is stable and the other is unstable. These stabilities can be derived from the sign of the above quantity at different fixed points. For example, when l E~6 1Hz, the quantities at the two fixed points are 0:94 and {4:06 respectively. Because the fixed points of P correspond to limit cycles, the system (8)D l E~l c E has a semi-stable limit cycle and the system (8)D l E §l c E has two bifurcated limit cycles: the one with a larger amplitude is stable and the other in the interior is unstable. In the simulation, the two bifurcated limit cycles can be numerically observed (Fig. 5B).
As shown in Fig. 5C-D, the interior limit cycle gradually shrinks to the equilibrium as l E increasingly departs from l c E to l h1 E &64:9Hz. When l E passes through l h1 E , a subcritical Hopf bifurcation occurs. The stable limit cycle is preserved, but the   shrinking interior limit cycle coincides with the equilibrium, and this makes the equilibrium unstable (Fig. 5E). The stabilities of the equilibrium and the limit cycle attributed to the Hopf bifurcation can be validated by calculating the first Lyapunov coefficient (FLC). The FLC for the bifurcation point l h1 E is 0:4721w0, which validates the existence of the subcritical Hopf bifurcation.
Interestingly, aside from the above two bifurcations, the other two bifurcations appear almost symmetrically and consecutively. When l E passes through l h2 E &90:9Hz, the other subcritical Hopf bifurcation of the system (8) emerges with a positive FLC, (0.6262), which changes the stability of the originally-unstable equilibrium and brings an unstable limit cycle (Fig. 5F-G). Moreover, the amplitude of the bifurcated unstable limit cycle grows until l E increasingly approaches l n E &99:6Hz, where the other saddlenode bifurcation occurs. This bifurcation leads to the coalescence and annihilation of the two limit cycles (Fig. 5H-I). The aboveexpatiated bifurcation procedure of the system (8) is illustrated in Fig. 7A.
As mentioned above, a burst is triggered if the spike threshold T is sufficiently low, and, because of the bifurcated stable limit cycle, there exists a stable periodic obit fluctuating between the two critical excitation levels. This indicates that bursts can occur continuously with an inter-burst interval that is equal to the period of the stable periodic orbit. Fig. 7B dynamically shows the spiking threshold and store level of the system (8) as the parameters are taken as l~90:9Hz and n~n 0 . Since t r is always set much larger than t OT , the sharp peaks in Fig. 7B (left) and the sharp valleys in Fig. 7B (right) reflect the characteristics of the slow-fast dynamical system.
For l E w99:6, the system (8) has no limit cycle but only one stable fixed point. In such a case, the dynamical behavior is analogous to that of the system with a small l E . Therefore, bursts disappear as excitation is beyond the critical level. From the perspective of the ONN, oxytocin is released so frequently that the stores are not replenished fast enough to reach the critical level required to trigger a burst. Under such conditions, bursts are rarer and less predictable, until eventually over-excitation disrupts the reflex secretion of oxytocin [23].
In Fig. 7C, the phase trajectories of the store level and T OT are plotted to show the bifurcation transition regulated by the input rate l E . Here, we fix nat 22. As shown in the inner plot of Fig. 7C, we start with stable attractor (the green star) with l E~5 7Hz. By increasing l E to 62 Hz, a value located in the bifurcation region as shown in Fig. 7A, the phase trajectory goes to a limit cycle (the purple curve). The bifurcation of bursting is generated. Further increasing l E so that the rate enters the high-rate stable attractor Hz; (E) l E~6 4:9 Hz; (F) l E~9 0:9 Hz; (G) l E~9 2:5 Hz; (H) l E~9 9:6 Hz; (I) l E~9 9:7 Hz. The r-nullcline and T OT -nullcline are colored in blue and green respectively. The red circles represent the unstable limit cycles, and the black curves stand for the orbits with the initial point (0,0). doi:10.1371/journal.pone.0038402.g005 region as shown in Fig. 7A, the system becomes stable again (the red curve and star), as the bursting activity is destroyed by overwhelming excitatory inputs. Decreasing l E to the bifurcation region (or, equivalently, increasing the inhibitory input rate) the system goes to a limit cycle (the blue curve) so that the bursting is induced. This transition coincides with the phenomena shown in Fig. 1F, that injections of inhibitory substances can paradoxically trigger bursts.

Comparing the Mean Field Model with the SNN
Based on the bifurcation analysis, we return to the network bursting dynamics and compare the SNN and the mean field model (Fig. 6B). In the SNN, a burst is recorded if the firing rate exceeds 30 Hz. For a network of 48 neurons and 12 dendrites, bursting emerges when the excitatory input frequency is between 15 Hz and 70 Hz. In the mean field model, we can replicate the network bursting dynamics from the traces of voltage and store level. For a givenn and l E , we say that there is network bursting if a stable limit cycle exists in the reduced system, and the inter-burst interval is the period of the limit cycle. Therefore, the burst frequency in the mean field model (Fig. 7D) is the reciprocal of the period of the limit cycle. For the SNN with 48 neurons and 12 dendrites, the value of n in the corresponding mean field model should be n~4N=n b~1 6. To compare with the SNN, we pick n~22 in the mean field model for reasons stated below. Fig. 7D shows that the replicated 'network' possesses similar bursting dynamics to the SNN.

Discussion
In the current paper, we present a general approach to tackle a complex neuronal network dynamics which exhibits oscillations at multiple time scales. Under the homogenous topological assumption of the network, the neurons display spiking activities induced by afferent inputs at the neuronal level, while the global network demonstrates synchronous oscillation at the network level. The ONN showed paradoxical network behaviors that the bursting events occur continually when the excitatory input rate is in a certain range, but disappear when the excitatory input rate is sufficiently large.
We developed a simplified version (SNN) of this model which preserves these basic behaviors. Then, we used the mean field approach and reduced the SNN to a mean field model, in which the bursting activity corresponds to a limit cycle. The critical step is the firing rate map approximation. We obtained the map via numerically simulating the leaky integrate-fire model with fixed threshold in each trial. A sigmoid-like function is then constructed to approximate the firing rate map.

Generality of the Approach
The main purpose of a mathematical model is to reveal the mechanism of a complex biological system, while retaining its main features. It is certainly unsatisfactory if we can only replace one complex (biological) system with another equally complex (mathematical) system. The ultimate aim of a mathematical model is to capture the essence of the system so that we can understand, interfere and control the system.
Our approach allows us to simplify a network model with a single neuron of 10 variables to a simple two-dimensional model: the mean field model. The generality of the approach is based on the facts that: 1. The original network has a routine configuration incorporating leaky integrate-and-fire model cells, spiking series represented by random processes, and a topological structure composed of coupled neural units, which is intensively used in many other biological modeling. 2. Using the oxytocin network as a vivid example, the framework of the procedure is common in the sense that the techniques, such as cutting off minor variables, approximating discrete spiking series, and employing bifurcation theory, are ubiquitous and often inevitably used in other network analyses. 3. It can be easily applied to other similar neuronal networks. The oxytocin model has attracted considerable interest, and other groups have tried to investigate it analytically as well [26]. Their approach is interesting, but does not address the actual mechanisms of the model. Another closely-related model is presented in [38] and its dynamical behavior should be very Here, x 0 (l E ) is the equilibrium point and L lE is transversal to the vector field in the neighborhood of x 0 (l E ). Note that both x 0 (l E ) and L lE depend continuously on l E ; (B) Curves of the Poincaré map P. Each intersection between the curves and the black line P 1 (a)~a corresponds to a fixed point of P as well as to a limit cycle of the system (8). For l E~5 9Hz, the curve has no intersection with the black line, so that there is no limit cycle. At higher values of l E , the curve moves upward; it first intersects with the black line at l E~l c E , where a single semi-stable limit cycle emerges. As l E increases to 61Hz, two bifurcated limit cycles appears. Here, one cycle is stable characterized by the quantity 1{DP'(a)D~0:94w0 at one fixed point, and the other cycle is unstable with the quantity {4:06v0 at the other fixed point. doi:10.1371/journal.pone.0038402.g006 similar to ours, as pointed out in our early paper [23], although we have not seen published work on it [39].

Parameter Choices
The mean field model is a two-dimensional dynamical system with two parameters, l E and n, where l E is the excitatory input frequency and n denotes the connection strength. For a fixed n §22, the dynamical system (8) displays two types of bifurcations as l E varies: a saddle-node bifurcation of limit cycle and a subcritical Hopf bifurcation. The former bifurcation accounts for the generation and ending of the bursting events and the identical inter-burst intervals, which is more significant in the network behaviors.
In the preceding investigations, bifurcations are studied with n fixed at 22. Actually, n is determined by the scale and the connection of the network. Numerical investigations of the system (8) show that, if n is small, there is no bifurcation for any possible values of l E . Indeed, for n,n o = 22, the system (8) has no limit cycles, but a unique stable equilibrium, i.e. no bursting activity appears in such a case. For a given n$n 0 and with the variation of l E , a pair of conjugate eigenvalues of the linearized system transversally cross the imaginary axis twice, so that the limit cycle generated by the Hopf bifurcation emerges. This makes it possible to generate stable limit cycles coexistent with unstable limit cycles, and explains why we picked n~22 for the mean field model in the comparison with the SNN with n~16 (Fig. 7D).
The bifurcations with respect to the input rate and network size can be summarized as the emergence of a codimension two bifurcation, namely the Bautin bifurcation [40], by regarding the mean field model as a member in the two-parameter family of autonomous ordinary differential equations. This bifurcation is beyond the scope of this paper.
We should point out that the phenomena based on the bifurcations of network size described above for the mean field model are not consistent with the ONN or SNN. As for the SNN as well as the ONN, even with a small network population or weak connection (i.e., each neuron is connected with few other neurons), bursting events still exist. Actually the mean field model might be more reasonable and closer to the underlying mechanisms of the real neuronal system in the sense that bursts could hardly be triggered for a single or few neurons. The discrepancy between the mean field model and the ONN tells us the shortcomings of the ONN model, despite successfully fitting of the model with experimental data.  Fig. 5F). The sharp peaks in the left plot and the sharp valleys in the right plot reflect the characteristics of the slow-fast dynamical system. (C) The bifurcation transition regulated by the input rate l E with n~n 0 . The inner plot indicates the dynamics of the input rate with respect to time. We set l E (t)~57Hz for t[½0,500(in seconds), l E (t)~62 Hz for t[½500,1100, l E (t)~200 Hz for t[½1100,1600 and l E (t)~90 Hz for t[½1600,2100. (D) Network bursting dynamics in: (blue line) the SNN composed of 48 neurons and 12 dendrites. (red line) the 'network' replicated from the traces of voltage and store level in the mean field model with n~22. Note that bursting events are recorded if the firing rate is .30 Hz in the SNN, while the burst frequency in the mean field model is the reciprocal of the period of the stable limit cycle. doi:10.1371/journal.pone.0038402.g007