Inferring network properties of cortical neurons with synaptic coupling and parameter dispersion

Computational models at different space-time scales allow us to understand the fundamental mechanisms that govern neural processes and relate uniquely these processes to neuroscience data. In this work, we propose a novel neurocomputational unit (a mesoscopic model which tell us about the interaction between local cortical nodes in a large scale neural mass model) of bursters that qualitatively captures the complex dynamics exhibited by a full network of parabolic bursting neurons. We observe that the temporal dynamics and fluctuation of mean synaptic action term exhibits a high degree of correlation with the spike/burst activity of our population. With heterogeneity in the applied drive and mean synaptic coupling derived from fast excitatory synapse approximations we observe long term behavior in our population dynamics such as partial oscillations, incoherence, and synchrony. In order to understand the origin of multistability at the population level as a function of mean synaptic coupling and heterogeneity in the firing rate threshold we employ a simple generative model for parabolic bursting recently proposed by Ghosh et al. (2009). Further, we use here a mean coupling formulated for fast spiking neurons for our analysis of generic model. Stability analysis of this mean field network allow us to identify all the relevant network states found in the detailed biophysical model. We derive here analytically several boundary solutions, a result which holds for any number of spikes per burst. These findings illustrate the role of oscillations occurring at slow time scales (bursts) on the global behavior of the network.


INTRODUCTION
The neuronal spike-burst activity is characterized by recurrent transitions between rest state and firing state where bursts are temporal groupings of multiple spikes. Certain cells in the mammal brain, for example, neurons in the thalamus during periods of drowsiness, attentiveness, and sleep are known to exhibit this type of spike-burst behavior (Sherman and Koch, 1986;Steriade and Llinás, 1988;McCormick and Feeser, 1990;Steriade et al., 1993;Amzica and Steriade, 1998). Autonomously bursting neurons are found in a variety of neural systems, from the mammalian cortex (Morris and Lecar, 1981;Dhamala et al., 2004a,b) to brainstem (Hindmarsh and Rose, 1984;Wang, 1994;Izhikevich, 2007;Jirsa and McIntosh, 2007;Jirsa, 2008). When neurons are coupled with each other, they produce different modes of behavior, including synchrony and phase-locking, which have been implicated in memory, cognition, sensory processing, motor planning, and execution (McCormick and Feeser, 1990;Wang, 1994;Jirsa and McIntosh, 2007). Many neurological diseases, on the other hand, including Parkinson, schizophrenia, and epilepsy, are the result of abnormal synchronization (Uhlhaas and Singer, 2006;Jensen et al., 2007), which suggests that a better understanding of the basic mechanisms producing synchrony and phase locking will be a stepping stone toward the repair of brain function. Modeling attempts using large scale networks to understand emergence of cognitive states rely heavily on the approximation of the dynamics as a neural ensemble. The concept of a neural mass like abstraction (Hebb, 1949;Beurle, 1956) designates a group of Co-activated neurons capable of acting like a closed system when performing a certain function. A small scale network of this kind is sometimes referred to as a "neurocomputational unit." In large scale brain networks, these mesoscopic units of operation serve as the network nodes (see for instance, Deco et al., 2008Deco et al., , 2011Ghosh et al., 2008). On intermediate spatial scales of few cm, neural activations along the spatially continuous cortical sheets are described by neural fields, for which the connectivity is assumed to be translationally invariant (see, Wilson and Cowan, 1972;Nunez, 1974;Amari, 1977;Jirsa and Haken, 1997;Feng et al., 2006;Jirsa, 2009;Robinson, 2011). To define such small neurocomputational units, simplified neuron models, known as phase models, offer an attractive tool for the study of network modes, since they allow for detailed mathematical analysis of network dynamics (Breakspear et al., 2010). As an example, Carbal et al. have explored the role of local network oscillations in resting-state functional connectivity by using such phase oscillators in the respective nodes of the simulated network. They have shown when these oscillatory units are integrated in the network, they behave as weakly coupled oscillators. Moreover, for a set of network parameters they found subsets of nodes tend to synchronize although the network is not globally synchronized (Cabral et al., 2011). For the present work we use a recently proposed phenomenological model that admits parabolic bursting in one dimension, which is a type of bursting observed in the R-15 neuron in abdominal ganglion of aquatic mollusc Aplysia Californica (Ermentrout and Kopell, 1986;Izhikevich, 2000;Ghosh et al., 2009). This type of bursting can arise even without bistability in the generation of spikes. The investigation carried out in this work with a detailed neuron model capable of displaying spiking and bursting behavior and a minimal model that not only reproduces the mean field amplitude of the original networks but also capture the most important temporal features of its dynamics. The detailed model used here is extensively discussed in Rinzel and Ermentrout (1989). On the other hand, our phase model is a minimal model that captures the generality of the mechanism of bursting present in the detailed model. As we vary network parameters including mean field coupling strength and dispersion, both networks display various temporal dynamics. In order to understand these states in mathematically tractable terms we take advantage of the mean field coupled network of phase model. Our goal is to identify to what degree this mean field model serves as a reliable neurocomputational unit and captures the qualitative features of temporal dynamics of the full network as a function of the investigated network parameters. Mean field analysis for singleton burst reveals solutions such as incoherence and partial oscillation which can be completely described analytically. However, as we are interested in a multispike system where analytical calculation is rather non-trivial and therefore, we combine semi-analytical approach with numerics to derive the stability diagram. Mean field phase network allow us to identify the mechanism of transitions between various network states that appear as solutions of the full network. Stability diagram is independent of number of spikes per burst and qualitatively commensurates well with the findings in our full network. The paper is structured as follows. In the next section, we introduce the Rinzel-Ermentrout model (Rinzel and Ermentrout, 1989) for parabolic bursting and describe the model in details. In the following section, we couple individual neurons via global coupling and present our analysis of this network model. In the subsequent section, we set up a generic network of bursters coupled to their mean field and derive semianalytically all the network states and corresponding phase transition boundaries. In the next section, we derive numerically a stability diagram using global phase coherence measure. In the final section, we summarize the results obtained from mean field descriptions and link them systematically with the network states obtained from biophysical model network.

SINGLE NEURON MODEL
A dynamical system with multiple time scales (for example, a neuron with spiking-bursting behaviors) can be written in a singularly perturbed form:ẋ = f (x, y),ẏ = rg (x, y), where x is the vector of fast variables, y the vector of slow variables that modulate the fast activity, and r 1 is a ratio of fast/slow time scales. A system which has been proposed to describe parabolic bursting behavior is known as Rinzel model (1989). Single neuron model parameters used here are exactly as described in Rinzel and Ermentrout (1989).
where V is the membrane potential, w is associated with the fast current, Na + or K + , Ca and n are the two slow currents, Model parameters which are held fixed throughout our simulations are, .0005, and g kCa = 1, g sCa = 1. I is the applied input current. The ionic currents are given by an ohmic leak current, determined by the leak conductance g l and leak reversal potential V l , and a Na + current which is responsible for the generation of spikes. The dynamics of this model which is relevant to our study is outlined as follows. When the input current I exceeds a critical value I c a single neuron described by Equation (1) undergoes a Saddle-node bifurcation on an invariant circle (SNIC). This same system for two different parameterization of I and in the presence of the slow currents can exhibit both spiking as well as parabolic bursting behavior. Spiking behaviors are elicited for a slightly higher value of the external drive. For example, to observe a typical burst-like pattern in this system we held the input current to the values I = 68 and for spikes I ≥ 70. Figure 1 displays the relationship between the applied input current and a parabolic bursting pattern that is observed in the single neuron dynamics.

PHASE MODEL
The generality of the underlying mechanism for parabolic bursting is investigated in details by numerous authors (Ermentrout and Kopell, 1986;Baer et al., 1995;Izhikevich, 2000). In many such formulations, parabolic bursting neurons are typically in their canonical form described as: where function f (x, y) in the above equation couples to spike generative mechanism depending on the slow variables x, y dynamics, respectively. The function f (x, y) is a smoothly varying periodic function alternating signs such that the system undergoes a SNIC to generate parabolic burst at the single neuron level.
Recently Ghosh et al. (2009) has also proposed a simpler model that in principle captures the underlying mechanism of parabolic bursting involving only a circular phase variable θ and moreover, involve only one slow term to allow the fast dynamics to enter or get out of repetitive firing. Motivation for using such a model is primarily mathematical tractability. Parameter space of this model cannot be directly linked to the biophysical parameters, however, qualitatively it may account for the transient and longterm behavior of more detailed biophysical models. In this model a single neuron is described by the following equation, In Equation (4) a slow variable activation term is represented by a modulation term cos( θ n ) which mimics the entire slow subsystem instead of describing it as a separate dynamical system, I is the applied input current and n is an integer, which determines the number of spikes per burst. In our simulation with this model all the results are for n = 5 spikes per burst unless otherwise specified. Golomb and Rinzel (1993) considered a heterogeneous network of all-to-all coupled inhibitory bursting neurons and found regimes of synchronous, anti-synchronous and asynchronous behavior when the width of the heterogeneity was changed (Golomb and Rinzel, 1993;Jirsa, 2008, 2011;Smeal et al., 2010;Jirsa and Stefanescu, 2011). We describe our network equations via a fast instantaneous coupling. N synaptically coupled (all-to-all) parabolic bursting neurons are described by a similar set of non-linear differential equations with fast chemical synapse. To this end we formally describe:

FULL NETWORK MODEL
where all the parameters and the gating variables inherit from the single neuron model Equation (1, 2) and b is a rescaling factor to applied drive to cross the threshold and elicit spike/burst in the uncoupled system. Stimulus that all the neurons see I i > 0 are drawn from a uniform distribution assumed to be symmetrically distributed over the interval I i ∈ [2.1 − I, 2.1 + I].
Where I is the spread of the applied stimulus parameter. I introduces a heterogeneity in the spike threshold. The synaptic coupling appears as an ensemble average given by S = 1 (1 + exp(−V i /2)) is a sigmoidal activation function.
The synaptic strength K is the same for all the neurons. For the entire simulation, we fixed the reversal potential of potassium ions to v th ≈ 0.0 (for purely excitatory connectivity).
Analysis is carried out for a fast synapse (AMPA-type glutamate receptors), such as those found in the auditory system, the rise time is instantaneous, and post-synaptic responses commence almost instantaneously after the start of presynaptic action potential (Nunez, 1974;Morris and Lecar, 1981). This brisk communication is a consequence of rapid calcium-channel kinetics, which allows significant calcium entry during the upstroke of the presynaptic action potential (Sabatini and Regehr, 1996). Under the fast synapse approximation the variable s i relaxes much more rapidly than V i , in which case we may apply a quasi-static approximation to (Equation 5) (e),ṡ i ≈ 0, allowing us to adiabatically eliminate the synaptic variable via s i = β (1 + β + exp(−V i /2)) . The time course of the postsynaptic conductivity caused by an activation of AMPA receptors can be captured by a rise time β rise = 0.09 ms and decay time β decay = 1.5 ms (Gabbiani et al., 1994;Parnas and Parnas, 1994). Numerical results in Figure 3 provides a good approximation for β in the range between [0.01 ms, 0.5 ms]. Although, we have provided here the details about the fast excitatory synaptic connectivity, our approach can be readily extended to inhibitory connectivity as well. In the continuum limit, a mean field formulation with inhibitory synaptic coupling is provided in details in Appendix.

MEAN FIELD COUPLED PHASE MODEL
Each generic neuron is coupled to this mean field and typically their response to the mean field expressed as R(θ) explicitly dependent on θ, and implicitly on time. In absence of any coupling, their vector field flow on a real line is governed by F(θ) = ω − cos(θ) − cos(θ)/n. In the absence of the term cos(θ)/n expression reduces to a mathematical description used in Roy et al. (2011). Together, we can write for N (still finite) such neurons:θ Recently, we have proposed a formulation for mean synaptic activation term under fairly general setting and taking advantage of instantenous activation, deactivation between pre and postsynaptic events. It allows one to describe synaptic activation variable as a non-linear transfer function of membrane voltage (Roy et al., 2011). Moreover, we have described how the mean field coupled spiking neurons can be described mathematically with this synaptic coupling. Details of this formulation is described elsewhere, (Roy et al., 2011). Collective activity of synapses is described by a mean field. For a given population of neurons is expressed more formally as, where is the mean field influence function. Coupling K is the same for all the neurons. In our previous work, response to such mean field coupling explicitly described as R(θ i ) = sin θ i (cos θ i − v th ), where O( ) contains non-circular deviations of the order that results due to perturbations. v th ≈ 0.0 for all simulations and analytical calculations unless mentioned otherwise. It is important to note that the couplings in the phase descriptions retain their mathematical expression in the full model plus some linearly added correction terms, which scale with the degree of order of deviation from the circle (Roy et al., 2011). Hence, in application it is rather suitable when phase perturbations are close to the circular orbit. The above equation further can be written combining the terms containing a single Fourier harmonic in the coupling plus the higher order Fourier terms.
See for details (Roy et al., 2011). Where, in Equation (6) the frequencies I i ≥ 0 are assumed to be symmetrically distributed over the interval I i ∈ [I − I, I + I] according to a uniform probability distributions.

CHARACTERIZATION OF SPIKE/BURST COHERENCE IN BIOPHYSICAL NETWORK MODEL WITH MEAN FIELD COUPLING
The bursting coherence and incoherence is quantitatively characterized in terms of a statistical-mechanical spike-based measure. We consider an excitatory population of neurons coupled to a common mean field drive and heterogeneity in their threshold for spikes/bursts. By varying the strength of the coupling K and the stimulus spread I we investigate the emergence of spike/burst coherence. Emergence of collective spiking/bursting coherence may be well described by the (population-averaged) global potential, In the thermodynamic limit (N → ∞), a collective state becomes is non-stationary (i.e., an oscillating global potential V mean appears for a coherent case), where the overbar represents the time average, and also, the correlated mean field (t) activity appears oscillatory. Otherwise (i.e., when V mean is time independent or stationary), it either becomes incoherent (IN) or partial oscillatory (PO). In N → ∞ limit both these states converges to a stationary solution. Thus, the mean square deviation of the global potential is a global marker for mean burst coherence for the entire population described here. More formally one can write it as (i.e., time-averaged fluctuations of V mean ), plays the role of an order parameter used for describing the coherence-incoherence transition (Manrubia et al., 2004). For the coherent (IN) state, the order parameter R(t) approaches a nonzero (zero) limit value as N goes to the infinity. We compute R(t) in Equation (12) as a function of mean field coupling strength K and dispersion parameter I for the full system. We vary both K, I from 0 to 1 in a step size of 0.01. Subsequently, computed values of R(t) is plotted in grid size of 100 × 100. Contour plot is colorcoded from low values at zero (blue) to high values at 1 (red). Nearly (in phase or anti phase) synchronized population spike/burst activity is lumped into a regime with labeled as SR and IN population spike/burst activity is lumped into a regime called IN activity. In the IN regime as described above R(t) values stays close to zero with substantial subthreshold fluctuations. Partial bursty regime is labeled as PO observed for R(t) values stationary and close to values other than zero. This regime displays dynamical behaviors far from synchrony, such as multi-clustering (some of the neurons are firing incoherently while others are not firing at all) in the phase for instance. Depending on the heterogeneity in stimulus spread we get random distribution of phases such that individual members can exhibit cluster hopping. Multiclustering in our model can reliably be captured using an ensemble average quantity rotation number ρ i given by Equation (14).

CHARACTERIZATION OF SPIKE/BURST COHERENCE IN PHASE NETWORK MODEL WITH MEAN FIELD COUPLING
The bursting coherence and incoherence is quantitatively characterized in terms of statistical mechanical order parameter coherence measure. As an alternative to storing and plotting many time series data θ i (t), i = 1, . . . , N for all N = 1000 variables, we define an order parameter Equation (13) measures the population dynamics. The advantage of using such a formulation becomes apparent immediately. Let's say our model system has periodic orbit then θ i (t) θ i (t + T), where T periodic pacing spikes or bursts (latency). Then in order parameter space one can can detect this state in a straight forward manner as a solution R θ (t)R θ (t + T). This result holds for all i, t. In this case, R θ dynamics is dominated mostly by the x co-ordinate dynamics. Absolute values of mean order parameter mod R θ ≤ 1. There is a mathematical relationship of macroscopic global phase measure with macroscopic V mean (t) in Equation (11). The interval between each microscopic spike/burst in an arbitrary ith stripe of spike/burst can be determined in a statistical-mechanical way by taking into consideration its contribution to the macroscopic global membrane potential V mean (t). In this interpretation, the time series of the global potential V mean (t) has a local maxima and minima, respectively and strictly bounded between [0,1]. The global cycle in the suprathreshold regime starting from the minimum of V mean (t) which appears first after the transient time is regarded as the first global cycle, which is denoted by G 1 . The 2nd global cycle G 2 begins from the next following right minimum of G1, and so on. Then, we can introduce an instantaneous global phase measure θ(t) of V mean (t) via a linear interpolation in the two successive subregions forming a complete global cycle (Lim and Kim, 2011). A microscopic spike makes the most constructive (inphase) contribution to V mean when the corresponding global phase θ k for kth cycle of spikes/burst is 2nπ (n = 0, 1, 2, . . .), while it makes the most destructive (anti-phase) contribution to V mean (t) when θ i for an arbitrary ith cycle of burst is 2(n1/2)π. By averaging the contributions of all microscopic spikes within a burst in the ith burst stripe to V mean , we can obtain the following degree of ordering of spikes/bursts. Hence, the contribution of kth microscopic burst occurring at the time t k is ordered by R θ (t k ). If the degree of synchrony is high between the bursts/spikes then R θ (t k ) → 1. We quantify the average firing frequency to compare the long-term behavior of individual neurons in the population model. We compute the average frequency (also known as the rotation number) of population of neurons using Averaging is carried out over about 1000 neurons starting from random initial conditions after the transient have died out.
Collective states of ensemble of N = 1000 neurons with spikes per burst n = 5 as indicated by their rotation numbers with uniform distribution of frequency I in the interval [2.1 − I, 2.1 + I].
Different branch of rotation index indicate different dynamical states of the network as a function of mean field coupling strength K, I. We carry out a grid search in the 2D parameter space K, I. Our goal is to obtain a phase transition diagram to understand long-term collective behavior of Equation (8) for large N, as a function of the coupling strength K ≥ 0 and the stimulus spread I ∈ [0, 1). Global order parameter R θ (t) is computed for different parameterization of K, I and embedded on a contour plot. Color spectrum is the same as the one used for displaying phase diagram in the full network. The values which are high and close to 1 are indicated by red and the values which are close to zero are indicated by blue.

CLUSTERING ANALYSIS IN N COUPLED FULL AND PHASE NETWORK MODEL
We describe firing patterns in large networks (finite N) with excitatory mean field coupling in terms of array diagrams. Array diagrams are obtained by simulating a coupled system consisting of mean field coupled biophysical neurons (N = 100) governed by the Equation (5). All the coupling coefficients are the same K where i = 1, . . . , N. In the arrays the intensity of the voltage variables V 1 , . . . , V i have been encoded in color spectrum. Two different color spectrums are used for the biophysical network (see Figure 4). In the first color spectrum blue part of the array values implies the quiescent activity of the spikes where the voltage variables have relatively lower values. All the other colors in the spectrum indicates the higher values for the voltage variables, consequently these pixels in the array imply the spike activity.
The horizontal line of the array shows the time with increasing epoches of activity. The second color spectrum used here shows burst depiction in the nearly coherent parameter regime. Green colors in the array indicate completely silent neurons. Purple pixels on the green background shows burst activity. On the vertical axis neuron index are aligned and again, on the horizontal axis gives the direction of time. These diagrams were obtained from a phase network by monitoring phases of individual neurons i = 1, . . . , N and aligning them on the vertical axis. The choice of the color spectrum used for phases is given by a colorbar with uniformly distributed phase values. In Figure 9 red color index in the spectrum corresponds to higher phase values of θ (close to π) and orange color index are for lower phase values (close to −π). First initial conditions θ i (0) is generated randomly and then they are sorted according to their neuron index and subsequently distributed uniformly about [−π, π]. The parameters K, I, for both realizations are chosen from SR, IN regime of the respective phase diagrams.

NUMERICAL PROCEDURES AND VISUALIZATION OF THE SYSTEM DYNAMICS
Two network models were implemented in Matlab, numerically integrated using second order Runge Kutta routine and Euler-Maruyama (EM) method (Higham, 2001). The simulations were performed with a fixed time step of dt = 0.05. The first 200 time points of the simulation are disregarded to set the network to a steady state. Thus, the results within this time were ignored. The membrane potential V(t), standard deviation of membrane potential std V, mean field (t), order parameters R(t), R θ (t) are captured for the entire population. For full network, simulation is carried out for N = 100 neurons and for the phase network for N = 1000 neurons. Numerical Phase diagrams are obtained using parallel for loops implemented in Matlab. Coupled mean field Phase model represented in Equation (8) can be visualized as a collection of N points rotating around the unit circle, where the estimated phase for each neuron θ i (t) denotes their position on a ring or a circle at time t. This alternate representation of the dynamical system (as N points moving along a circular reference frame, instead of a single point tracing out a trajectory in an N-dimensional phase space) is possible because the system's state space, the N-torus, is equivalent to N copies of the unit circle. It is worth noting that for most other N-dimensional state spaces such a reduced representation is not feasible. In order to distinguish between oscillators with different natural frequencies, we color the dots according to the standard color spectrum: the neurons correspond to the low end of the spectrum (close to -π)(red), neurons at the high end (close to π) (blue), while those in between occupy the middle part of the spectrum (orange/yellow/green).
To show how the system evolves from one instant to the next, we plot a series of snapshots of the system at different times (see Figures 11B-D, for example). This allows us to observe the behavior of individual neurons at the same time as we witness the collective evolution of the system toward an attractive state.

SINGLE NEURON BURST DYNAMICS
We first examine the behavior of single neuron model Equation (2.1) as the applied input current I is brought close to the threshold for generating spikes or bursts. For the given parameters In Equation (2.1) a neuron is excitable. Figure 1 depicts the relationship between applied input and parabolic bursting pattern. We are only interested in the behavior of this system for low current values where the resting state of membrane voltage is sufficiently depolarized below −40 mV. When the applied input current I is below a critical value membrane potential V(t) maintains their steady state value and for values greater than the threshold exhibits bursting behavior. For the parameterization used here we find that at I ≥ 60 steady state destabilizes exhibiting multispikes. When the applied input current is further increased a neuron make transition from bursting to spiking behavior. In order to observe a typical spiking behavior we set I ≥ 70. To get an intuitive understanding about the relationship between slow and fast subsystems, Rinzel and Lee analyzed this model by varying ca (a variable in the slow subsystem) as a bifurcation parameter to report that parabolic bursting is obtained from an oscillation in the slow subsystem that periodically moves the ca variable back and forth across the SNIC bifurcation, to link the steady state solution of this system to (quiescence state in Figure 1A) the branch of periodic solutions ( Figure 1B) and vice versa. Time series of fast variable shows that the interspike interval is relatively longer at the beginning and end of each burst. As has been shown by numerous authors oscillation for the fast dynamics is obtained when the slow variables are held fixed; it is where the saddle-node-loop bifurcation occurs. There is a clear threshold below which there is a unique stable fixed point. Parabolic bursting can occur without having any bistability in the spike generating process. One way to achieve parabolic bursting behavior without requiring any bistability in the generating process and moreover, mathematically tractable would require a generic description like the one shown in Equation (2.2) (see section 2). From numerical results we find that as the applied input current I → 2, time period T → ∞. Applied input current can be tuned such that it is possible to obtain parabolic bursts of desired interburst gap. The time evolution of a single neuron activity is shown in Figure 2, where a membrane potential like variable V(t) = − cos[θ(t)] is plotted by numerically integrating Equation (2.2). Temporal dynamics shows regular parabolic bursting behavior. For I < I c a neuron fires few spikes before it settles into a steady state.
For I > I c (I = 2.01, n = 5) neurons exhibits parabolic bursting behavior. Based on the qualitative similarity in the burst pattern with parabolic bursting neurons (At the start and the end of the active phase the spike frequency is smaller compared to the middle of the active phase as can be seen in Figure 2 detailed model is substituted to investigate the network effects.

NETWORK DYNAMICS OF PARABOLIC BURSTING NEURONS WITH HETEROGENEITY
To understand the influence of heterogeneity in the applied input current and the coupling strength in a network of single neurons exhibiting parabolic bursting we use Equation (5) and parameters as described in section 2.4. We use fast excitatory synapses to couple these units. When the synaptic coupling is sufficiently fast, the coupling tends to push the neurons toward anti-synchrony (Wang and Rinzel, 1992;Friesen, 1994;Van Vreeswijk et al., 1994). Moreover, several studies have observed emergence of multistable solutions in their mean field network with parameter heterogeneity (Assisi et al., 2005;Jirsa and Stefanescu, 2011). Our motivation is to go toward this particular direction to capture the relevant network dynamics at the population level. In particular to understand the combined effect of heterogeneity in the firing rate threshold (biophysical model) with the fast time scale of activation-deactivation of synapses in the coupling; the interplay between these two critical factors in spike/burst timing at the population level is largely unknown. In our formalism their individual and combined influence on the network dynamics become clearly visible. Typical time course of such responses of synaptic variable in our model simulation are shown in Figure 2). Fast synapse approximation holds as long as the variable s i relaxes much more rapidly than V i , in which case we may apply a quasistatic approximation to reduce s i further in Equation (5)  states which fires in a wave-like pattern. Corresponding time series of mean quantities such as the membrane potential V mean (t) in Equation (11), mean field (t) shows periodic activity in Figure 4C. Membrane potential spiking activity is nearly synchronized across population of neurons in Figure 4. On the other hand, for the IN state mean membrane potential fast decays to zero and shows subthreshold fluctuations about mean zero. Response of mean membrane potential is more suppressed compared to their mean field oscillations between [0,1]. Amplitude of mean field (t) changes in time systematically but fluctuates about the mean value of 0.5 instead of approaching zero values as can be seen in Figure 4D. Population burst synchrony is observed for many different parameterization, for one such choice of parameter I = 0.002, K = 0.8, an array diagram is computed and plotted in Figure 5A. As can be seen in the figure a wavelike spread of activity. In Figure 5B various time series plots of population burst synchrony is shown across 10 neurons. In order to identify different network states for all possible combination In the incoherent state mean field amplitude systematically decrease in time.
of two parameters K, I we carry out a grid search and compute the values of R in Equation (12). Global order parameter measure identifies three distinct network states in the parameter space as shown in Figure 6A.  Figure 6C all neurons are showing nearly synchronized parabolic bursting behavior.

NETWORK DYNAMICS OF PARABOLIC BURSTING PHASE MODEL WITH HETEROGENEITY
In this section, we use a phase network with mean field coupling to get some insights about the novel network states observed in (K, I) the parameter space of the full network model. Coupling between each phase neuron via a mean field is formulated in section 2.4. Numerically we integrate Equation (6) to compute time averaged membrane potential, mean field as in Equation (10) (see section 2), global measure of coherence R θ as a function of K, I a parameter combination which is used in the detailed network model. Time evolution of the above quantiles are shown in Figure 7 for a parameterization K = 0.8, I = 0.001. The parameter choice is the same as the full network investigations. With this combination of parameters all the neurons synchronously spikes. Mean membrane potential-like quantity V mean (t) oscillates in phase with synchronized spike activity as plotted in Figure 7B. Here, n a quantity which determines the number of spikes per burst is kept at n = 1. Mean field also shows up and down states (Locked in time) and act as an oscillating drive to each individual neurons. The time series of the global order parameter R θ (t) for synchronized spiking is periodic in Figure 7E. Next, we show in Figure 8 temporal evolution of the mean quantities for the choice of K = 0.8, I = 0.5. For medium values of mean field coupling strength and stimulus spread network shows PO behavior, where some of the neurons are firing incoherently and others are completely silent. This means for some parameterization network has two stable branches of solutions. It is important to note PO state of the network was observed in the full network for a comparable parameterization (see Figure 6). Time series for 10 neurons and their order parameter evolution in time is plotted in Figure 8. Three neurons are completely silent while other seven neurons are bursting with variable inter-burst intervals. As there is no noise in this system and coupling magnitude is set at high values as in the case of sync, this variability must be introduced by the heterogeneity in their individual firing rate threshold via stimulus spread. Figure 6 shows the parameter space diagram for the full and phase models presented in Equations (5) and (6-9). Phase boundaries are calculated by computing the mean field for both full and the phase model for different combination of (K, I) values on a two-dimensional grid. In the following subsection we would lay out the details for obtaining the phase transition boundaries semi analytically. Interestingly, Over a wide range of (K − I) values the collective dynamics of the two networks primarily show three distinct regions of interest which are close to each other in the parameter space. For sufficiently large K values holding I fixed to a narrow range of values near zero, the system converges to a state of partial oscillations in which the some of the neurons are not firing altogether, while the others display IN oscillations. Very large K values result in damping of oscillation activity and all the neurons stops firing altogether. The stability state of locking is much more difficult to achieve and infact we found distinct branches in their rotation number, these states should all be regarded as variants of 1:1 locking, and therefore we lump them together in the locked region of the stability diagram. With further increase in I, parameter heterogeneity, successively more neurons peel away until eventually the entire population is IN.

PHASE DIAGRAM USING SEMI ANALYTICAL METHODS FOR MEAN FIELD PHASE MODEL
Mean field coupled neurons in phase model is described in Equation (6). Let's rewrite the mean field equation explicitly.
In a semianalytical approach we would like to understand the phase transitions between three distinct network states discovered in two networks. For the IN states where the average firing frequency increases monotonically plotted in Figure 10, the θ i are all distributed across the closed orbit in a unit circle. This leads to the following phase evolution equation SR state may undergo instability either through parameter changes of K or I and make phase transition to either IN or PO state. Mean field approaches a stationary density as the number of neurons are increased in both PO and IN state (see Figure 4D). Hence, approaches some positive real number for these two states. When varying K, we consider small perturbations μ to the SR solution θ = θ i = 0. With θ = θ i = 0 + μ Equation (15) becomesθ i =μ = F(μ) + 2 sin (2μ) and linearization yieldṡ θ i = (F (0) + )μ. Moreover, SR state may gets phase locked at θ = π (subpopulation clusters). Hence, θ = π may get destabilized as we changed the width of heterogeneity by changing I or the coupling strength K. Similarly, we consider small perturbations μ about solution θ = π. Hence, we can write θ i = π + μ,θ i =μ = F(π + μ) + 2 sin(2π + 2μ) and linearization yieldsθ i = (F (π) + )μ. With F(θ) ≈ I − cos θ − cos θ n for the SR state, we find that F (θ) ≈ − sin(θ) − 1 n sin θ n will be generally small for θ = 0, π. SR state solution hence becomes unstable when F (0, π) + = 0, which suggest almost a vertical critical line between SR and PO, IN state. The bifurcation route from PO (multistable state) to IN solutions as the parameter I increases is less conclusive in the framework of the circular approximation, since in the previous stability analysis the only I-dependent term is F (0, π), which is very small, hence higher orders of the approximation must be considered. We use the following ansatz: If r is the radius of a unit circle, any smooth deformation from a unit circle can be approximated as, r(θ i ) = 1 + h(θ i ). Hence we can compute the non-linear flow contribution with the above first order correction term as F(θ i ) + H(θ i ). It is possible to explicitly determine H(θ i ) for a certain choice of h(θ i ) and moreover, H(θ i ) has a periodicity of π, that is H(θ + π) = H(θ). Thus the linear stability analysis about the fixed point θ From the above equation with F (0) ≈ 0 and the π-periodicity of H(θ i ), we find that the two fixed points at 0, π lose stability at the same time for increasing I and as a result leads directly to the IN state. Since H (θ i ) ∼ I, scales linearly for fixed μ, we can also estimate the critical line of transition in the parameter space in Figure 11 which separates PO state from IN state.  I c ), which serves as a convenient guide to numerically compute the stability line separating PO region from IN. Next we try to obtain analytically the stability boundary between INC and PO oscillation states in the infinite-N limit. it turns out that the IN and partial oscillation states can be made steady in our system. The possibility of doing so was suggested by the numerical results. In numerics we observed that as the number of neurons N is increased, the order parameter R t approaches a constant for both these states Figure 8 and the oscillators tend to arrange themselves in a stationary distribution around the circle Figure 11. The way to approach these two states analytically, therefore, is to first write down the appropriate infinite-N analog of our model.
The above equation is the infinite-N analog of continuity equation introduced earlier. It is a non-linear partial integrodifferential equation for the number density f (θ, t, ω). In addition we demand f to be non-negative, 2π periodic in θ, and we impose the normalization 2π 0 f (θ, t, I)dθ = 1, For incoherence and partial oscillation the above system tends toward a stationary distribution of phases in time. The above two states are the fixed points of the stationary density in the continuum limit. To solve for the fixed points we set ∂ ∂t f = 0 in Equation (10). let's assume that f 0 (θ, ω) be the stationary phase density and v 0 = [F(θ) − (( 2π 0 I + I I − I f (θ , t, I )g(I )dI dθ )) sin(2θ)] be the velocity field. Then one can write where L(I) is a constant which is determined exactly by using normalization condition. Depending on it's applied drive I, neuron's steady state behavior falls in the following two categories: Case (ii) When I >> neuron fires incoherently and typically individual phases follows an uniform distribution about the unit circle. In this case the velocity field turns out to be, Fixed point solution demands that the density must be inversely proportional to the velocity: In the IN state, neurons driven by different external drives are firing at different phases, however, their collective state is close to being stationary. Every neuron belong to Case (ii) as described above. Further, it is possible to derive nearly an exact relationship between K, I that gives the transition from case (i) to case (ii) as described above. As shown before, in case of a finite size network such a relationship in the first order perturbation m(I − I c ) p + = 0 does exist. In this scenario those neurons with a minimum bound on their applied drive I min reach cessation of firing as we find from numerical simulations. They then fall into the Case (i) above where mean field exerts much bigger influence on the dynamics and overall effect is damping of firing activity. The first neurons to stop firing are the ones which do not cross the threshold for firing which in this case I > 2. Then the boundary that separates IN from PO in the phase diagram is almost a straight line given by, solution we can express a relationship between K and I as follows, Equation (24) gives us an empirical relationship between parameters upto fourth order perturbations for the bifurcation of a limit cycle. Optimization of the above equation gives coefficients a 0 = 8 π , a 1 = 0, a 2 = 128 π 3 , respectively. Next, we substitute the amount of dispersion I into the perturbative term to obtain the boundary between IN and SR state. Taken together we can write, Results are shown in Figure 12 in the (K, I) plane using Equations (25) and 23. Critical lines obtained semi analytically qualitatively agrees well with the numerical results that captures various network states in both these models with purely excitatory coupling. In Appendix, we show a stability calculation for an inhibitory coupled mean field network in the infinite analog limit. From numerical simulations we find that the results are independent of the number of spikes n per burst.

DISCUSSION
One of the most frequent assumption employed in simulations of large neural networks is that the whole network can be lumped into small aggregates of collective unit (sometimes called a "neurocomputational unit") exhibit a sufficiently similar dynamical behavior. Consequently, the network that instantiates this ensemble, consisting of thousands of excitatory and inhibitory neurons, it is considered to display a synchronized behavior with no other significant temporal features for the dynamics of the large scale network. The main reason for this assumption, is the impractical large computational time arising from too many details considered in the large network properties. In this paper, we have analyzed the behavior of a neural network that serves as a good example of such a unit, namely a mean field coupled bursting ensemble. First, we have investigated a Hodgkin-Huxley type detailed biophysical model widely employed in theoretical and computational neuroscience with global coupling. We found that the dynamical features of the network are far more complex then the ones corresponding to synchronized or rest state behavior. The network dynamics depends critical on the balance between firing rate threshold dispersion and mean field synaptic coupling strength; in fact, the synchronized state can be found only for a specific range of parameters typically involving a large or medium values for the coupling strength and low values of dispersion. On the other hand, for large dispersion and weak coupling strength values both networks display purely IN behavior. In the IN state, individual neurons are driven by different external drives results in firing at different phases, however, their collective state is close to being stationary. This stationarity in the density distribution led us to formulate scaling relationship between coupling strength and dispersion parameter. One interesting finding is that, when mean field exerts a greater influence than parameter dispersion; it causes shutting down of the neural activity in some neurons. In this parameter range, we find interesting dynamical behavior such as partial activity. In order to address the problem of the high computational cost of such an implementation, we have further developed a selfconsistent mathematically tractable mean field coupled phase model following (Assisi et al., 2005;Ghosh et al., 2009;Jirsa and Stefanescu, 2011), but incorporating a higher degree of realism. Rather than finding the most appropriate type and number of dimensions that could minimize certain error functions or capture statistical variance in the full network, we have focused our attention on understanding a phenomenological burst generation model system which captures the most important network dynamics of bursting units at the population level. Collective activity of synapses is described by a mean field which relies on instantaneous rise and decay time (Roy et al., 2011). This mean field is then employed in the coupling to individual neurons to describe phase network. Together, we investigate this population of neurons coupled to a common mean field drive and heterogeneity in their threshold for spikes/bursts. Our detailed analysis demonstrated that the reduced representation manages to recreate correctly the topology of the mean field amplitudes of the original system for various parameter scenarios. In the full network, In the thermodynamic limit (N → ∞), a collective state becomes coherent if δV mean (t) ≡ [V mean (t) − V mean (t)] is non-stationary (i.e., an oscillating global potential V mean appears for a coherent case) and also, the correlated mean field (t) activity appears oscillatory. In the phase network, global order parameter is computed by averaging the contributions of all microscopic spikes within a burst in order to obtain a similar degree of ordering of spikes/bursts as in the full model for identical parameterization. Hence, for a dynamical behavior such as coherence-incoherence transition macroscopic order parameter gives us a crude approximation of burst timing. From a more general perspective, despite its limitations because of the consideration of purely excitatory or inhibitory network, it can be readily extended to study networks with mixed coupling. Moreover, the analytical approach to estimate the scaling relationship and transition boundaries between the IN-PO-SR states is not restricted to small scale network only. With global coupling, as the size of the network grows the boundaries may shift leading to a different parameterization than the one displayed here; however, underlying bifurcations remain the same. We have demonstrated this in our work by analytically deriving a low dimensional mean field amplitude reduction for a inhibitory coupled mean field network in the continuum limit. In this case, all the relevant dynamics of an infinite dimensional network in Equations (29) and (30) is captured by a two dimensional representation of the reduced mean field population given by Equation (40). Thus, using this approach, we derive analytically a low dimensional representation of the network dynamics and we show that the main features of the neural population's collective behavior can be captured well by the dynamics of a few cortical nodes exhibiting spiking as well as bursting behavior. While it is true that strong reductionist assumptions are common (sacrificing many of the biological realism of a network node's dynamics) in large-scale network modeling, these assumptions are usually made ad-hoc on the network node's dynamics and limit the network dynamics to a small range. We emphasize here that because of the "near to synchrony" assumption, neural mass models cannot capture complex dynamical features such as multiclustering, oscillator death or multi-time scale synchronization. Evidently a reduced small scale network model is desirable to serve as a node in a large scale network simulation whereby displaying a sufficiently rich dynamic repertoire. Here it is of less importance to find a quantitatively precise reduced description of a neural population; rather more importantly, we seek a computationally inexpensive population model (this means typically low-dimensional) which is able to display the major qualitative dynamic behaviors (synchronization, rest state, multi-clustering, etc.) for realistic parameter ranges as observed in the total population of neurons. Our approach may offer a viable alternative to the neural mass models currently used in the literature. By comparison, our model offers the possibility to account for such features (temporal details of their spiking activity considered irrelevant for the dynamics of the large network) at a very low computational cost. Therefore, the type of reduced representation discussed in this paper qualifies as a good candidate for a "neural unit" in computational simulations of large scale neural networks.
The above equation can be expressed in polar coordinates if we substitute z1 = ρ 1 exp(iφ 1 ) giving evolution equations for ρ 1 and φ 1 dρ 1 dt = 2 ρ 1 1 − ρ 2 1 − ρ 1 + F 2 1 − 2ρ 2 1 cos φ 1 + 1 2 cos (φ 1 ) + Fρ 1 2 ρ −1/n 1 − ρ 1/n 1 cos(φ 1 /n) For the Lorentzian distribution function the above equation is exact, However, we do not find any deviation of the above results for any other unimodal distributions of our firing threshold (such as uniform distribution) such as the one considered in the numerical simulations with excitatory coupling. The above two dimensional system can be solved numerically to identify the full network states and the corresponding transition boundaries.