An Efficient Population Density Method for Modeling Neural Networks with Synaptic Dynamics Manifesting Finite Relaxation Time and Short-Term Plasticity

Abstract When incorporating more realistic synaptic dynamics, the computational efficiency of population density methods (PDMs) declines sharply due to the increase in the dimension of master equations. To avoid such a decline, we develop an efficient PDM, termed colored-synapse PDM (csPDM), in which the dimension of the master equations does not depend on the number of synapse-associated state variables in the underlying network model. Our goal is to allow the PDM to incorporate realistic synaptic dynamics that possesses not only finite relaxation time but also short-term plasticity (STP). The model equations of csPDM are derived based on the diffusion approximation on synaptic dynamics and probability density function methods for Langevin equations with colored noise. Numerical examples, given by simulations of the population dynamics of uncoupled exponential integrate-and-fire (EIF) neurons, show good agreement between the results of csPDM and Monte Carlo simulations (MCSs). Compared to the original full-dimensional PDM (fdPDM), the csPDM reveals more excellent computational efficiency because of the lower dimension of the master equations. In addition, it permits network dynamics to possess the short-term plastic characteristics inherited from plastic synapses. The novel csPDM has potential applicability to any spiking neuron models because of no assumptions on neuronal dynamics, and, more importantly, this is the first report of PDM to successfully encompass short-term facilitation/depression properties.


Introduction
The population density method (PDM) appears as a time-saving alternative to direct numerical simulations [i.e., Monte Carlo simulation (MCS)] of large-scale neural networks of spiking neurons (Knight et al., 1996;Brunel and Sergi, 1998;Brunel, 2000;Nykamp and Tranchina, 2000;Omurtag et al., 2000;Casti et al., 2002;Huertas and Smith, 2006;Augustin et al., 2013). In this method, biophysically similar neurons are grouped into a population, and a probability density function for each population is defined to describe the variability of state variables of such neurons. The evolution of this density function is governed by a master equation (generally a partial differential equation such as a Fokker-Planck equation). A prominent advantage of PDM is that if the master equation can be easily solved in analytical or numerical ways, precise estimates of the quantities of interest can be obtained much faster than MCS.
Because of the significant impact of synaptic dynamics on network dynamics, PDMs must ultimately include accurate descriptions of synaptic kinetics to be widely applicable. Unfortunately, incorporating more synaptic dynamics increases the dimension of density functions, leading to a catastrophic decline in computational efficiency because of the increasing difficulty of solving highdimensional master equations. How to incorporate realistic synaptic dynamics in a computationally efficient manner therefore becomes an outstanding problem in the field of PDMs (Apfaltrer et al., 2006;Tranchina, 2009).
Although many previous studies addressed this issue (Haskell et al., 2001;Nykamp and Tranchina, 2001;Apfaltrer et al., 2006;Rangan and Cai, 2006;Ly and Tranchina, 2007;Rangan et al., 2008), they left out some unsolved problems (section 7.10 in Tranchina, 2009). Recently, Ly (2013) developed a principled dimensionreduction method, called modified mean-field method (MMFM), that could deal with those unsolved problems. Nevertheless, we found that the MMFM might not be a robust method and not suitable for nonlinear spiking neuron models. The previous studies concentrated, to the best of our knowledge, only on synaptic dynamics with finite decay and neglected another important synaptic attribute, the short-term plasticity (STP; Markram and Tsodyks, 1996;Markram et al., 1998;Zucker and Regehr, 2002) that also played an important role in certain specific network dynamics, such as reverberatory bursting activity in in vitro and in vivo neuronal networks (Volman et al., 2007;Gritsun et al., 2011), self-organized criticality (Levina et al., 2007), anticipative neural responses (Fung et al., 2012), and sustained population activity related to working memory (Barak and Tsodyks, 2007;Mongillo et al., 2008). Therefore, our goal in the present study is to develop a new solution that allows the PDM to include realistic synaptic dynamics of both finite decay and STP in a computationally efficient manner.
In this study, we present a novel and efficient PDM, called colored-synapse PDM (csPDM), to achieve our goal. The strategy behind the csPDM is also to take advantage of dimension-reduction methods, relying on the following scheme: (1) synaptic dynamics are assumed to operate in the diffusive limit so that the dynamic equations of neuronal dynamics are converted into Langevin equations with colored driving noise and all synapseassociated state variables are transformed into input terms (Moreno-Bote and Parga, 2004;Richardson and Gerstner, 2005; Destexhe and Rudolph-Lilith, 2012), and (2) the probability density function method for this dynamic system is then employed (Wang et al., 2013;Barajas-Solano and Tartakovsky, 2016). As such, the synapse-associated state variables are not viewed as the system state variables of the resulting Langevin equations, leading to the reduction of the dimension of density functions. As a consequence, incorporating realistic synaptic dynamics does not decrease the computational efficiency when adopting csPDM, because it does not increase the dimension of master equations.
We compare the simulated network activities from the csPDM with those from MCSs, original full-dimensional PDMs (fdPDMs), and MMFMs. The numerical results show that the results of the csPDM reveal good agreement with those of the MCS and fdPDM in both the steady-state and dynamic regimes. The csPDM provides more accurate simulations than the MMFM and remains more computationally efficient than the fdPDM. We hope that the csPDM will be valuable to computational neuroscientists who seek accurate and efficient network modeling tools to explore effects of synaptic dynamics with/ without STP properties on macroscopic behaviors of neural networks.

Materials and Methods
Here, we describe the equations and parameters used for neuronal, synaptic and network modeling. We present the model (master) equations of the fdPDM, csPDM, and MMFM. The corresponding closed-form solutions for evaluating firing rate responses (network activities) are also described. At the end of this section, we present the numerical method used to solve the model equations of the csPDM as well as the methods used to quantify simulation errors. A summary of the notation for all the main dynamical variables and physiological parameters is given in Table 1.

Network models
The large-scale network considered here is a neuronal population of N uncoupled exponential integrate-and-fire (EIF) neurons (Fourcaud-Trocmé et al., 2003) receiving only synaptic inputs from the outside. We choose the EIF model mainly because it can offer better prediction of spiking times for a given current input than the leaky integrate-and-fire models and quadratic ones (Fourcaud-Trocmé et al., 2003). In addition, its adaption version, i.e., the adaptive EIF model (Brette and Gerstner, 2005), is capable of reproducing various firing patterns observed in in vivo recordings of cortical cells Touboul and Brette, 2008). In fact, the csPDM derived in the present study can be directly applied to simulate neural networks of adaptive EIF neurons.
In this network, the external synaptic inputs to each neuron are assumed to be mediated by m different types of neurotransmitter receptors. For the sth type of receptors, there are presumably c s synaptic connections. In other words, each neuron has totally ͚ sϭ1 sϭm c s external synaptic connections. The synaptic input on each connection is a set of unitary spiking events whose random arrival times are assumed to be governed by an inhomogeneous Poisson process. s (t) denotes the mean rate of the Poisson process for the input through the sth-type receptor. to satisfy basic premises of the PDM, two arbitrary synaptic inputs are assumed to be statistically independent if they are mediated by different types of receptors and to be statistically identical if mediated by the same type.

Single neuron models
For each neuron, the dynamics of its membrane voltage, V (t), is described by where I syn ͑V, t͒ is the total synaptic current generated by the external synaptic inputs, C is the membrane capacitance, g l is the leak conductance, E l is the leak reversal potential, is the sharpness factor, and V T is the threshold voltage. The spiking mechanism is the following: a spike is triggered at t s when V (t s ) reaches a cutting voltage, V c (in general, V c Ն V T ). Afterward, V is immediately reset to a resetting voltage, V r , and, then, clamped for a refractory period, ref .

Synapse models
Specifically, I syn ͑V, t͒ is the sum of all postsynaptic currents mediated by various types of neurotransmitter receptors: where g s (t) is the (collective) synaptic conductance for the input mediated by the sth type of receptors, and E s is the corresponding synaptic reversal potential. Equation 2 refers to conductance-based synaptic models.
The synaptic dynamics means the dynamics of the g s (t). In our treatment, the synaptic conductance follows first order kinetics, where the rise in conductance at the arrival of a unitary synaptic event is instantaneous and the decline is exponential. Therefore, the dynamic equation for the synaptic conductance has the following form: where s is the decay or synaptic time constant and ⌫ s is the total variation of the synaptic conductance induced by a single event. The size of the instantaneous jump in conductance equals to ⌫ s / s . t j s is the arrival time of the jth synaptic event, and ͚ j ␦͑t Ϫ t j s ͒ refers to the spike train consisting of synaptic events summed from all synaptic connections through the sth-type receptor.
In our numerical examples, the value of ⌫ s is set according to the following equation: where ⌬V s is a free parameter defined as the instant voltage modulation provoked by a single synaptic event when the membrane voltage is E l just before the event arrives in cases where (in Eq. 1) is in the limit of ¡ 0 and the synaptic dynamics has a zero-order kinetics, i.e., s ¡ 0.

Short-term plastic synapses
What stated above does not consider STP property on synaptic connections. In the following, we concentrate on how to add it into the synaptic dynamics. A typical dynamic model of STP was the one proposed by Barak and Tsodyks (2007). They introduced two running parameters to describe the change of the connection weight of a synaptic coupling in a phenomenological way. Following the same idea, the ⌫ s in Equation 3 is replaced by where the running parameters, u s and x s , are governed by at the level of a single synapse. The jump size of the synaptic conductance caused by a spike now depends on u s and x s . u s is the running utilization parameter and x s is the running fraction of available neurotransmitters. The dynamic value ranges from U s to 1 for u s and from 0 to 1 for x s . U s refers to the base level of u s . u s ϭ 1 means a presynaptic spike is allowed to use all available neurotransmitters, and x s ϭ 1 means all neurotransmitters are available. Briefly, the Equations 6, 7 depict a phenomenon, where when a synaptic spike arrives, u s increases by an amount of U s ·͑1 Ϫ u s ͒, and, afterward, decays to its baseline level, U s , with the facilitation time constant, s,f , and, meanwhile, x s decreases by an amount of x s u s and recovers to its baseline value of 1 with the recovery time constant, s,r , afterward. In general, s,f and s,r are short, ranging from hundreds of milliseconds to seconds (Markram and Tsodyks, 1996;Markram et al., 1998). Besides s , E s , and ⌫ s , each synaptic connection with the STP property is also characterized by three additional parameters: U s (initial utilization parameter) and s,f as well as s,r (facilitation and recovery time constants, respectively), which control the type of the STP property from strong depression ( r Ͼ f and relatively large values of U s ) to strong facilitation ( f Ͼ r and small values of U s ; Barak and Tsodyks, 2007).
Note that we keep the subscript s in the Equations 5-7 above, because it is likely that the synaptic connections through different types of transmitter receptors have different STP types.

MCS
We are interested in the network activity, i.e., the output population firing rate, of such a neuronal population in response to a given set of s ͑t͒ (sʦ͕1, 2, ѧ, m͖) as a whole.
This quantity can be computed by the direct use of the MCS, where states of all neurons and synapses are explicitly traced according to Equations 1-3 or Equations 1-7 if STP properties are involved. For all our simulation examples, MSC is performed with the brian2 simulator (Stimberg et al., 2014). In all simulations, the network composes of N ϭ 10,000 EIF neurons, and c s ϭ 200 for all receptor types.
In the MCS, the output population firing rate is calculated by r͑t͒ ϭ n s ͑t, t ϩ ⌬t͒ N⌬t , where n s ͑t, t ϩ ⌬t͒ is the total number of spikes produced by all neurons in the population within the time period of ͑t, t ϩ ⌬t͒, and ⌬t is a small time window. ⌬ t is always set to 1 ms in our simulation examples.

fdPDM
The assumptions and conditions, under which the fd-PDM provides an exact description of network activity, have been previously discussed (Nykamp and Tranchina, 2000;Omurtag et al., 2000).
Specifically, the fdPDM is a density-based modeling approach to directly estimate r(t) with a known distribution of states of neurons, which is described by a joint probability density function, ͑V, g , t͒, defined by ͐ ⍀ ͑V, g , t͒d⍀ ϭ Pr͕͑V͑t͒, g ͑t͒͒ ʦ ⍀͖ , (9) where g ϭ ͓g 1 , g 2 , ѧ, g m ͔ is the conductance state vector composed of m synaptic conductance variables. ⍀ is a subdomain within the state space bounded by the definition domains of all state variables. In general, the definition domain of ͕g sԽ sʦ͕1, 2, ѧ, m͖͖ ranges from 0 to ϱ. About the membrane voltage V, its upper bound is V c , and its lower bound is given by min͑E l , V r , E s ͒ Խ sʦ͕1, 2, ѧ, m͖.
Equation 9 states that the integration of the density function over ⍀ is the probability of finding neurons whose states are within that subdomain in a large neuronal population. Note that, here, we do not consider STP properties on synaptic connections to illustrate our main ideas in the present study in a simple way. We will consider these properties only for the csPDM later.
Unlike in the MCS, where one has to track every neuron individually, in the fdPDM, one just track one density function for each population. As explained in details in previous papers (Nykamp and Tranchina, 2000;Haskell et al., 2001;Apfaltrer et al., 2006;Ly, 2013), the master equation governing the evolution of ͑V, g , t͒ forms a (1ϩ m)-dimensional (in space) partial-integral differential equation: where the voltage probability flux, i.e., the probability flux along the voltage-direction, J V, is given by and J g ϭ ͓J g 1 , J g 2 , ѧ, J g m ͔ is the conductance probability flux vector composed of m conductance probability fluxes. The conductance probability flux, J g s , is given by The population firing rate, r(t), is obtained by integrating the voltage probability flux across V c over conductance, yielding ⍀ gЈ denotes the subdomain where J V ͑V c , g , t͒ Ͼ 0, meaning that neurons at V c intend to cross the boundary from below to generate spikes in this domain. As seen, for a given set of inputs, s (t) (sʦ͕1, 2, ѧ, m͖), we need to solve the Equation 10 to obtain the output response of the neuronal population, i.e., the output population firing rate, r(t). For solving that equation, the boundary conditions are: where V lb denotes the lower bound of the membrane voltage. The former boundary condition results from the resetting mechanism of spiking, meaning that the voltage probability flux across V c , accounting for the generation of action potentials, re-enters the state space on V r after the refractory period. The latter ones mean no probability fluxes flow outward the domain through other boundaries. As stated, the boundary conditions result in the conservation of the number of neurons. That is to say, the integration of the density function over the domain always equals to one in the absence of refractory period.
The computational time required for solving the master equation depends on the number of its dimensions. Generally, solving a high-dimensional master equation is highly time consuming. Previous studies (Apfaltrer et al., 2006;Tranchina, 2009) have suggested that the fdPDM cannot be considered as a time-saving alternative to the MCS if the dimension of the master equation exceeds three. In the following, two other methods used to tackle this issue are described. MMFM is a documented method, and csPDM is a new method proposed in this paper.

MMFM
MMFM is proposed by Ly (2013) for speeding up the estimation of the population firing rate. In this method, g in the Equation 10 is viewed as a deterministic parameter so that the master equation becomes a one-dimensional partial differential equation: In this case, a steady-state firing rate for a given g is obtained by using the mean-field model, given by Recall that the synaptic inputs for different types of neurotransmitter receptors are assumed to be independent mutually; therefore, synaptic conductances for different receptor types are also independent with each other. Furthermore, they only depend on their individual inputs. So, it is realized that the value of g s evolves isolatedly. For a large neuronal population, one can define a marginal conductance density function, s ͑g s , t͒, for g s by analogy with ͑V, g , t͒. Also, the integration of s ͑g s , t͒ over a certain space on g s is the probability of finding neurons whose values of g s are located in that space. According to what derived by Ly (2013), the master equation governing the time evolution of s ͑g s , t͒ is yielded as where The boundary conditions for solving Equation 16 are J s ͑0, t͒ ϭ J s ͑ϱ, t͒ ϭ 0. Assuming that the values of all elements of g are drawn from the marginal conductance density functions individually, MMFM estimates the timevarying firing rate, r͑t͒ , with expected value of r ͑g ͒ conditioned on g , that is, As such, instead of solving a (1 ϩ m)-dimensional partial-integral differential equation in the fdPDM, one only needs to solve m one-dimensional partial differential equations simultaneously for estimating r(t) in the MMFM. Ly suggested that MMFM had higher computational efficiency than fdPDM, especially when the values of r at all possible points of g could be computed before simulations.

Colored-synapse population density method (csPDM)
This section presents the derivation of the csPDM. csPDM is inspired by the probability density method for Langevin equations with colored noise (Wang et al., 2013;Barajas-Solano and Tartakovsky, 2016). So, we start by presenting how to transform the Equation 1 into a stochastic Langevin equation through the diffusion approximation of g s .
In the diffusion limit, under the condition that the mean rate of synaptic input received by each neuron, c s s , is sufficiently high and the jump size ⌫ s / s for each synaptic spike is small enough, the random synaptic conductance g s ͑t͒ can be treated as the Ornstein-Uhlenbeck process (Uhlenbeck and Ornstein, 1930;Risken, 1996), whose dynamics is given by (Richardson, 2004) where g s is the mean value of g s ͑t͒ and g s is the standard deviation of g s ͑t͒. ͑t͒ is a ␦-correlated white-noise process of unit variance. Introducing a new variable g s ϭ 1 / g s ͑g s Ϫ g s ͒ and substituting for g s in the equation above, we obtain a new Ornstein-Uhlenbeck process for g s : g s has zero mean and unit variance (Risken, 1996); and its autocorrelation function is As a result, g s is a colored noise due to the exponential form of C s ͑͒ and finite synaptic time constant, i.e., s Ͼ 0. Replacing g s ͑t͒ in the Equation 1 with g s ϩ g s g s , we obtain as the new dynamic equation for the membrane voltage, where g ϭ ͓ g 1 , g 2 , ѧ, g m ͔, g ϭ ͓ g 1 , g 2 , ѧ, g m ͔, and g ϭ ͓g 1 , g 2 , ѧ, g m ͔. So, as shown in Equation 22, the fluctuation of V͑t͒ now is characterized by a Langevin equation with m independent colored noises ͕g sԽ sʦ ͕1, 2, ѧ, m͖͖ resulting from synaptic dynamics with finite synaptic time constants (Destexhe and Rudolph-Lilith, 2012). That is why we use "colored-synapse" as the prefix of this method. Barajas-Solano and Tartakovsky (2016) proposed a closed-form quasi-Fokker-Planck equation as the master equation for Langevin equations driven by colored noise. They considered a dynamic system characterized by the stochastic Langevin equation in n dimensions Each v i ͑x, t͒ is decomposed into a deterministic function or "mean-field velocity" ͗v i ͑x, t͒͘ 0 and a stochastic fluctuation term v i Ј͑x, t͒. The authors proposed a quasi-Fokker-Planck equation: as the master equation governing the temporal evolution of the joint probability density function of the system states, p͑x, t͒. The stochastic diffusion tensor, D, is calculated by if fluctuation velocity components are exponentially autocorrelated and mutually uncorrelated, i.e., J T ͑x͒ is the Jacobian of the mean-field velocity with components J ij ͑x, t͒ ϭ Ѩ͗v i ͑x, t͒͘ 0 /Ѩx j . Considering V, Ᏼ 0 ͑V, g ͒, and Ᏼ 1 ͑V, g , g ͒ in Equation 22 as x, ͗v i ͑x, t͒͘ 0 , and v i Ј͑x, t͒ in Equation 23, respectively, we can yield the master equation for the marginal voltage density function, V ͑V, t͒, by using Equation 24, as the following: where (27) By using Equation 25, we obtain and ͗V͘ denotes the average membrane voltage across the population. What remains in Equation 26 are equations for evaluating g s and g s . To derive the equation for evaluating g s ϭ ͐ 0 ϱ g s s ͑g s , t͒dg s , we multiply Equation 16 by g s and integrate g s from 0 to ϱ: Similarly, one can derive the equation for evaluating For g s , we use the identity, g s ϭ ͙ g s 2 Ϫ g s 2 . In the csPDM, the population firing rate is directly estimated by the probability flux across V c : The average membrane voltage, ͗V͑t͒͘, is computed by using its definition: They are similar to what used in the fdPDM except for the last boundary condition, which means that no neuron locates at V c because the neurons whose voltages reach V c are reset to V r immediately.
When the STP property is included, we replace ⌫ s and ⌫ s 2 with in Equation 31, Equation 32 for evaluating g s and g s , where u s and x s are mean values of u s and x s across the neuronal population, respectively. Following mean-field equations proposed by Barak and Tsodyks (2007), we use the same equations for tracking u s and x s (i.e., Eq. 1 in Barak and Tsodyks, 2007): To sum up, unlike in the fdPDM, the network dynamics in the csPDM is described by a system consisting of a one-dimensional quasi-Fokker-Planck equation for tracking the marginal density, V ͑V, t͒, and 2m ordinary differential equations for tracking the first two statistical moments of all synapse-associated variables, ͕g s ͖, or 4m equations for statistical moments of ͕g s , u s , x s ͖ if STP is included. csPDM is expected to be more computationally efficient than fdPDM because solving one-dimensional quasi-Fokker-Planck equation undoubtedly takes less time than high-dimensional partial differential equations. Basically, they are solved numerically. In the next section, we present the numerical method used for solving the quasi-Fokker-Planck equation.

Numerical method for solving quasi-Fokker-Planck equation
To solve Equation 26, we use local Galerkin method (LGM; Cockburn and Shu, 1998;Xu and Shu, 2010). This method belongs to the discontinuous Galerkin methods and focuses on the solutions of partial differential equations with high order derivatives. We select it as the proposed numerical method because it has superior ability to handle discontinuous solutions (Huang et al., 2015) and parallelizability in computations (Biswas et al., 1994). First, in the LGM, an auxiliary variable, q͑V, t͒, is introduced to rewrite the Equation 26 as followings: As a consequence, the Equation 26 is transformed to a conservative hyperbolic equation.

Discretization of space domain and basic notations
For a given bounded V-domain I ϭ ͓V lb , V c ͔, we divide it into M meshes with an identical length as follows: so that V r and E s are exactly certain grid points. We denote the subspace I k ϭ ͑v kϪ 1 / 2 , v kϩ 1 / 2 ͒ (k ϭ 1, 2, ѧ, M) and its length L ϭ ͑V c Ϫ V lb ͒/M. Using first-order polynomials as shape functions, the approximated values of V ͑V, t͒ and q͑V, t͒ within I k are defined by: in which 1 and 2 are shape functions. kϩ 1 / 2 Ϫ ͑t͒ and kϩ 1 / 2 ϩ ͑t͒ refers to the value of V at v kϩ 1 / 2 from the left mesh I k and from right mesh I kϩ1 , respectively. Due to the discontinuity at the interface of adjacent meshes, kϩ 1 / 2 Ϫ kϩ 1 / 2 ϩ is thus possible. q kϩ 1 / 2 Ϫ ͑t͒ and q kϩ 1 / 2 ϩ ͑t͒ are defined in the same way.

Slope limiters
A slope limiter is employed to guarantee the positivity of the density function (Huang et al., 2015). After progressing one time step with the backward Euler method, the across all meshes go through a slope limiter ⌳⌸, which is defined by (47) where k is set as 1 / 2 ͑ kϪ 1 / 2 ϩ ϩ kϩ 1 / 2 Ϫ ͒. The is the "minmod function" defined as: Generally, Equations 10, 16 are also solved numerically. We adopt the so-called discontinuous Galerkin method to solve them. The details of this method are not stated here. Please refer to Huang's paper (Huang et al., 2015) for details about the numerical method.

Quantification of simulation errors
MCSs are considered as the ground-truth of all simulation examples in this study; thereby, the differences of the fdPDM, MMFM, and csPDM in simulation results with the MSC are used to explore their performances. The following two quantities are used to quantify the simulation errors of the csPDM on the estimations of marginal conductance and voltage density functions, given by to check the validity of diffusion approximation on synaptic dynamics. g csPDM describes the average error ratio on the marginal conductance density function, which is estimated only over the interval ͑ g s Ϫ 6 g s , g s ϩ 6 g s ͒. It is noted that the marginal conductance density function in the csPDM is characterized by a Gaussian distribution, here denoted by N͑g sԽ g s , g s ͒, whose mean and standard deviation are g s and g s . V csPDM is the average error ratio on the marginal voltage density function.
The third quantity is

Results
We will present four simulation examples: (1) steadystate analyses, (2) dynamic population responses to excitatory input only, (3) dynamic population responses to multiple inputs, and (4) STP population responses. Using the former two examples, we aim to compare the performances of fdPDM, MMFM, and csPDM and focus on the accuracy and efficiency of csPDM. We show the applicability of csPDM to more complicated neuronal activities by the latter two examples. Thus, in both the examples, only the simulation results of MCS and cdPDM are presented. The parameters required for simulations are listed inTable 2. The values of neuronal properties are assigned according to previous studies (Fourcaud-Trocmé et al., 2003;Destexhe, 2009;Hertäg et al., 2014), and all of them are within in the physiologic ranges. The parameters relevant to synaptic dynamics and STP are regarded as free parameters which are freely tuned to let s correspond to the time scales of common transmitter receptors or to make synaptic interactions exhibit short-term facilitation/depression if necessary. The input rate, s , shown in this section is normalized by the minimum constant input rate that drives EIF neurons to produce spikes under fluctuation-free conditions (i.e., mean-field models), except for the third and fourth simulation examples. Such a minimum input rate is g l ͑V T Ϫ E l Ϫ ͒/͓ ͑E s Ϫ V T ͒⌫ s c s ͔. Next, we start to show simulation results by first checking the validity condition of diffusion approximation of synaptic dynamics.

Validation of diffusion approximation of synaptic dynamics
For the correct use of the csPDM, to explore the validity condition of diffusion approximation of synaptic dynamics is necessary. To do so, we check steady-state values of g csPDM and V csPDM of a single neuronal population receiving only constant excitatory inputs with respect to different s and ⌬V s . The synaptic time constant s is set as 5 or 100 ms to match the time scale of AMPA-receptors or NMDAreceptors (Dayan and Abbott, 2001), and E s is set as 0 mV for both cases. Basically, small g csPDM reflects the validity of diffusion approximation. As shown in Figure 1A, larger s or smaller ⌬V s results in smaller g csPDM for both s ϭ 5 and s ϭ 100 ms, in agreement with the fact that the diffusion approximation is valid when the input rate is large and the jump size is small. As seen, increasing s produces small g csPDM because of the decreased variation of g s ( g s ϰ 1/ ͙ s at steady states). Three distributions of g s under different sets of ( s , ⌬V s ) are shown in Figure  2A-C. As shown, the distribution of g s , s ͑g s ͒, approaches a Gaussian distribution when g csPDM is Ͻ0.2 (Fig. 2C). That is, the diffusion approximation of synaptic dynamics is valid if g csPDM Ͻ 0.2. Figure 1B shows V csPDM . The amplitude of V csPDM is proportional to that of g csPDM , meaning that the error of csPDM in V ͑V͒ in large part comes from the failure of the diffusion approximation. But, surprisingly, the value of V csPDM is much lower than that of g csPDM under the same set of ( s , ⌬V s ). As seen in Figure 2A,D, under the case where s ϭ 0.1, ⌬V s ϭ 1 mV, and s ϭ 5 ms, V csPDM is only 0.174 whereas g csPDM is as high as 0.653. However, the reason why V csPDM is smaller than g csPDM is unclear. Through the observation of V ͑V͒ (Fig. 2D-F), it is found that csPDM gives an adequately accurate estimation of V ͑V͒ comparable to MCS when V csPDM is Ͻ0.2. To offer a quantitative and universal validity condition for diffusion approximation of synaptic dynamics, we also check the value of g s / g s (i.e., the coefficient of variation of g s ; Fig. 1C). In fact, it has been argued that the diffusion approximation is valid only when g s / g s Ͻ Ͻ1, meaning that the conductance mean g s should be much larger than the standard deviation g s (Richardson and Gerstner, 2005;Cai et al., 2012). It is however an impractical condition. Based on our observations of g csPDM and V csPDM , it can be conclusively said that the validity condition of diffusion approximation for the application of csPDM is g s / g s Յ 0.6, under which g csPDM and V csPDM are approximately Ͻ0.2. As a result, ⌬V s should be Ͻ1.5 mV for the case of s ϭ 5 ms such that this criterion can be satisfied in the range of s Ն 0.6. In the following simulations, we choose ⌬V s as 1 mV for the excitatory inputs.

Steady-state analyses
The performances of fdPDM, MMFM, and csPDM are first examined via steady-state analyses of population firing rates in response to fixed excitatory inputs. Here, we also consider two cases, s ϭ 5 and s ϭ 100 ms, and set ⌬V s ϭ 1 mV and E s ϭ 0 mV. Figure 3 displays input-output curves computed by MSC, fdPDM, csPDM, and MMFM. Results show that csPDM gives an estimation of inputoutput curves that are close to those from MSC and fdPDM with errors Ͻ1 Hz for both s ϭ 5 and s ϭ 100 ms. These results indicate that csPDM can accurately estimate steady-state output population firing rates although it just tracks the marginal voltage density function. In contrast, MMFM gives accurate estimations of steadystate output firing rates only under the condition of e ϭ 100 ms. It overestimates the actual firing rate in the fluctuation-driven regime (i.e., s Ͻ 1) and underestimates the actual firing rate in the mean-driven regime (i.e., s Ͼ 1) in the case of s ϭ 5 ms (Fig. 3A, bottom panel, pink dashed line with triangle markers). The overestimation or underestimation of the MMFM in the case of e ϭ 5 ms was also observed by Ly (2013). However, the demon-

Dynamic population responses to excitatory input only
In addition to the steady-state analyses described in the above, the comparison between csPDM and MMFM is also made by investigating their ability to capture dynamic population firing rates in response to time-varying inputs. This is the most essential test to reveal whether they are adequate as dimension-reduction methods. Similar to the previous cases, excitatory inputs with s ϭ 5 and s ϭ 100 ms are considered here. The time-varying input rate is shown in Figure 4A, which varies in time and takes values that guarantee the validity of diffusion approximation, i.e., s Ն 0.6. Figure 4B shows the output population firing rates r(t) computed by MCS, fdPDM, csPDM, and MMFM and the corresponding errors in r(t) for s ϭ 5 ms. Figure  4C is similar to Figure 4B except s ϭ 100 ms. As shown in these two panels, whatever synaptic time constant is, csPDM accurately captures all the qualitative features of the output firing rates, leading to small average error ratios, which are 0.048 and 0.043 for s ϭ 5 and s ϭ 100 ms. It also gives comparable simulation results to those of fdPDM (similar average error ratios between fdPDM and csPDM). However, like in the steady-state analyses, MMFM gives accurate output firing rates only when s ϭ 100 ms. It overestimates the low output firing rates (at about t ϭ 820 ms) and underestimates the high output firing rates (at about t ϭ 730 or t ϭ 850 ms) when s ϭ 5 ms, leading to a large error ratio r MMFM of 0.179.
To show how robust the csPDM is, we explore whether the accuracy of csPDM depends on the values of model parameters. Four parameters, ⌬V s , g l , and V r , are cho- sen for this test. The same input rate as shown in the Figure 4A is considered as the input. As shown in Figure  5A, augmenting the jump size ⌬V s increases r csPDM , meaning that the accuracy of csPDM is decreased. This is due to the fact that diffusion approximation becomes invalid when ⌬V s is too large (as shown in Fig. 1). As expected, r fdPDM is less sensitive to the change of ⌬V s than r csPDM because ⌬V s is arbitrary in the fdPDM. Surprisingly, although MMFM does not have limitations on ⌬V s , in the case of s ϭ 5 ms, r MMFM severely fluctuates over a range of 0.2 when changing ⌬V s . Figure 5B-D unravel that changing g l , and V r almost does not affect the accuracy of csPDM ( r csPDM is below 0.1 in most parameter sets), implying that csPDM is a robust method. Importantly, simulation results from csPDM are comparable to those from fdPDM. MMFM can give comparable results to csPDM only when s ϭ 100 ms. However, in the case of s ϭ 5 ms, reducing g l enlarges the error ratio r MMFM with an increasing magnitude of 0.3 (from 0.07 to 0.38). So does increasing (from 0.1 to 0.4). In other words, MMFM is not a robust method for EIF models. Figure 6 shows the computational time spent by fdPDM and csPDM subject to different numbers of meshes employed in numerical methods for a simulation of 1 s. The numerical simulation used in this test is the same as shown in Figure 4. Remarkably, the computational efficiency of csPDM is ϳ1000 times better than fdPDM. Low computational efficiency of fdPDM certainly comes from the existence of the g s -dimension in the master equation because, except for the requirement of more grid meshes along this di-mension (We set 120 meshes along this dimension), its existence leads to the necessity of extremely small time steps for satisfying the Courant-Friedrichs-Lewy condition to ensure numerical stability of discontinuous Galerkin methods (Huang et al., 2015). In this case, the time step is 0.02 ms for fdPDM but 0.2 ms for csPDM.
In the next two examples, we only use csPDM to estimate the population dynamics and compare the results with those of MCS to evaluate its performance because, via the examples above, it has been demonstrated that, compared to fdPDM and MMFM, only csPDM can provide simulation results efficiently and accurately.

Dynamic population responses with excitatory and inhibitory inputs
To highlight the outstanding ability of csPDM to achieve dimension reduction, here, we consider a real case where a population of cortical pyramidal neurons receives excitatory inputs from the neighboring pyramidal neuronal population mediated by AMPA receptors and, meanwhile, inhibitory inputs from the interneuron population mediated by two types of receptors, GABA A and GABA B (Suffczynski et al., 2004;Cona et al., 2014). To mimic this situation, we consider a single uncoupled population of EIF neurons, representing pyramidal neurons, which receive external excitatory inputs mediated by AMPA-type receptors, which have s ϭ 5 ms and ⌬V s ϭ 1 mV, and inhibitory inputs mediated by GABA A -type receptors, which have s ϭ 10 ms and ⌬V s ϭ 0.25 mV, as well as GABA B -type receptors, which have s ϭ 100 ms and ⌬V s  Figure 7A, top panel, shows the time-varying input rates on excitatory and inhibitory synaptic connections. It should be noted that the inputs gated by GABA A and GABA B receptors share the common input rate in statis-tical sense. The Figure 7A, second panel from top, is the raster plot of output spikes from 200 neurons in MCS. With respect to the output population firing rate, as seen in Figure 7A, third panel, csPDM is able to capture the salient features of the population firing rate. The average error ratio r csPDM is ϳ 0.067 (below 0.1). It is ϳ 4 Hz of error in the high output firing rate (such as at t ϭ 1900 ms; Fig.  7A, forth panel) on average. csPDM also can accurately estimate the average membrane voltage across the population (Fig. 7A, fifth panel). Figure 7B shows the snapshots of the marginal density function V ͑V, t͒ at t ϭ 1700, 1800, and 1900 ms. As seen, the density functions computed by MCS and csPDM match well in all three snapshots. The corresponding values of V csPDM are 0.065, 0.07, and 0.122, respectively. These results indicate that V ͑V, t͒ can be correctly captured by csPDM.

STP population responses
In this part, we test the effects of the synaptic STP on population responses to see if csPDM can capture these effects. We consider a single population of EIF neurons receiving only excitatory inputs through plastic synapses with short-term facilitation or depression. For the facilitative synapses, we adopt s ϭ 5 ms, ⌬V s ϭ 2 mv E s ϭ 0 mV, s,f ϭ 700 ms, s,r ϭ 100, and U s ϭ 0.05. They are the same except s,f ϭ 50 ms, s,r ϭ 300, and U s ϭ 0.2 for the depressive synapses. Figure 8A shows the time-varying input rate, which consists of a fixed baseline of 500 ms and, afterward, a sinusoidal fluctuation of 1 s. The population responses and ensemble averages of u s and x s obtained from MSC and csPDM for the case of short-term facilitation are plotted in Figure 8B. As seen, the population response gradually increases in the latter cycles in the sinusoidal fluctuation (800 -1500 ms). The facilitative response actually stems from the sharp increase in u s from 0.05 to ϳ0.29. Meanwhile, the value of x s remains at a high level ( x s Ϸ 0.76, eventually). Figure 8C is similar to Figure 8B except that depressive synapses are used. As shown, the population response is depressed along the sinusoidal input, resulting from the large decrease in x s roughly by 0.4, during which u s remains almost fixed. Note that, different to the above example, where csPDM can give quantitatively accurate estimates of population responses, in this case, csPDM just captures STP property of population responses qualitatively. As seen in Figure 8C, the output firing rate obtained from csPDM drops but does not vanishes, while it vanishes in MCS. This difference is caused by the mismatch of x s between MSC and csPDM. In summary, the population response can inherit the STP property from synapses, and csPDM captures such property qualitatively.

Discussion
In this study, we present a principled and straightforward dimension-reduction method for PDM to handle realistic synaptic dynamics and compare it with another dimension-reduction method, called MMFM proposed by Ly (2013). We name the newly proposed reduction method csPDM. The csPDM does not assume specific limits on synaptic time constants so that we can consider synaptic dynamics mediated by all kinds of common receptors, including AMPA, GABA A , GABA B , and even NMDA. Through our examples, it is demonstrated that csPDM can accurately capture the firing rate responses in both the steady-state and dynamic regimes over a large range of synaptic time constants from milliseconds to hundreds of milliseconds.
As seen in Equation 24, the resulting quasi-Fokker-Planck equation leaves out all synapse-associated dimensions and just tracks the marginal density function of membrane voltages across the population. As a result, csPDM is extremely computationally efficient even when many types of receptors are included in the network model. The computational speed of csPDM is much faster than that of the original fdPDM by an order of three under the case where one type of receptors is incorporated (Fig.  6). If three types of receptors are incorporated, like in the third example in this study, fdPDM can never be considered as a computationally efficient modeling tool to simulate neural networks because of the inherent difficulty of solving a four-dimensional master equation. Actually, previous studies have argued that solving a master equation with more than three dimensions may spend much more time than MCS (Apfaltrer et al., 2006;Ly and Tranchina, 2007). However, csPDM can efficiently and correctly provides simulations of network activities in this case (Fig. 7). On the other hand, increasing the number of receptor types in the network model does not reduce the accuracy of csPDM considerably. r csPDM is 0.048 when one type of receptors is included (Fig. 4) and is 0.067 when three types are included (Fig. 7).
When csPDM is compared with MMFM, some interesting findings are observed. First, it is found that MMFM gives quantitatively accurate simulation results only for long synaptic time constants in the steady-state and dynamic regimes (Figs. 3, 4). A possible explanation is that the membrane time constant (C/g l ) is longer than the synaptic time constant ( s ) so that the use of the expected value of the steady-state firing rates for estimating the output firing rate (Eq. 18) becomes unreasonable because the neuronal dynamics cannot reach a steady state before synaptic dynamics does. This idea is also evidenced by the decline in r MMFM when g l is increased, i.e., reducing the membrane time constant (Fig. 5B).  By robustness analysis, we found that the accuracy of MMFM is sensitively influenced by the value of (Fig. 5C), which controls the nonlinearity of EIF neuronal dynamics by adjusting the strength of exponential currents (Eq. 1).
Increasing , i.e., increasing the nonlinearity, decreases the accuracy of MMFM. Actually, MMFM relies on the assumption that the membrane voltage inherits randomness entirely from the stochastic conductance variables,  The population response is facilitated during the sinusoidal input, resulting from the sharp increase in the running utilization parameter, u s As shown, the csPDM can capture the STP property of the population responses not only qualitatively but also quantitatively. C, Population responses under short-term depression. The population response is depressed, resulting from the sharp decrease in the running fraction of available neurotransmitters, x s . In this case, csPDM only captures the STP property of the population responses qualitatively. The larger error in population responses is due to the overestimation of x s and the conductance variables evolve independently of the membrane voltage (Ly, 2013). However, the author claimed that such an assumption might only be reasonable for linear models such as leaky integrate-and-fire neurons. So, MMFM may not be suitable to simulate network dynamics of nonlinear spiking neurons, like the EIF neurons used in this study. On the contrary, csPDM only assumes that the conductance variables evolve isolatedly, which is true no matter what neuron models are used. Accordingly, it is more general for the type of neuronal models. Furthermore, in the csPDM, the membrane voltage is allowed to be influenced by other random sources, for example, the random release of neurotransmitters (Faisal et al., 2008;Destexhe and Rudolph-Lilith, 2012), which may just add extra noise terms into the stochastic diffusion tensor in the quasi-Fokker-Planck equation without increasing its dimension.csPDM relies on the diffusion approximation of conductance variables. Physiologically, the assumption of diffusion approximation is quite reasonable in that the number of synaptic contacts covering a typical cortical pyramidal cell is indeed very high, in the range of 10 3 -10 4 (Braitenberg and Schüz, 1991). Each spike at a synaptic contact only alters a post-synaptic potential of Ͻ1 mV on average (Markram and Tsodyks, 1996;Markram et al., 1998). So, in our cases, each excitatory or inhibitory spike is designed to cause a maximum change of Ϸ1 mV in the membrane potential.
Our study is similar to the work of Rudolph and Destexhe (2003), in which the authors employed Itoĉ alculus to formulate the Fokker-Planck equation for Langevin equations and, then, to obtain an analytic expression for the membrane potential distribution at the steady state. They found that steady-state solutions are only available under a limited range of parameter sets. To correct this problem, they used an extended analytic expression for the membrane potential distribution (Rudolph and Destexhe, 2005), by which new solutions became available for a wider range of parameter sets, although the parameters resided in the nonphysiological extreme limits. The major difference of our approach with theirs is the method used to derive the Fokker-Planck equation. In our study, it is based on a modified largeeddy-diffusivity closure method proposed by Barajas-Solano and Tartakovsky (2016), and a quasi-Fokker-Planck equation, i.e., Equation 24, governing the temporal evolution of the membrane potential distribution, is then derived. Although we do not test whether Equation 24 is valid for any parameter sets, by means of our examples in this study, it is found that it at least provides a good solution for parameter sets in physiologic ranges. However, surprisingly, when we apply the same derivation procedures to the neuronal model adopted by them (i.e., Eq. 4.1 in Rudolph and Destexhe, 2005), we obtain the same Fokker-Planck equation as what they derived, in which the diffusion coefficient is also rescaled by 2 m e / ͑ m ϩ e ͒ ( m is the membrane time constant at rest). This fact implies that Equation 24 might be adequate for an even wider range of parameter sets.
It is theoretically possible to extend our method to any neuron models because no premises are put on the neuronal dynamics. We have tested it on adaptive EIF models (Brette and Gerstner, 2005) and leaky integrate-and-fireor-burst models (Smith et al., 2000) and find that csPDM works well (data not shown). Such an extension to adaptive EIF neurons is important for the development of PDM in computational neuroscience because the adaptive EIF model has been shown to reproduce a variety of spiking patterns observed in in vivo thalamic or cortical neurons Touboul and Brette, 2008) and well predict spike times in response to external inputs (Jolivet et al., 2008). csPDM with adaptive EIF models therefore will be a critical research subject in computational neuroscience. We suggest that the synaptic dynamics should be restricted to the first-order kinetics because the calculation of the diffusion tensor, D or D͑V, g s ͒, is based on an assumption that fluctuation velocity components are exponentially autocorrelated. The calculation becomes more difficult if higher-order synaptic dynamics is used. csPDM only captures STP property of population responses in a qualitative way due to incorrect estimation of x s . This poor estimation is caused by a wrong assumption that u s and x s are statistically independent for deriving mean-field equations for u s and x s . In fact, u s is independent of x s because x s does not appear in the Equation 6, but, by contrast, x s has a negative correlation with u s . As a consequence, the mean-field equation for x s should include higher-order statistics of both u s and x s . However, such an equation is unclosed. Before a closure method is available, the assumption of independence between u s and x s is a compromised solution . Note that it is possible to replace the phenomenological model of STP used here with another elaborating model where state variables are independent with each other. If so, the STP property of population responses can be quantitatively captured by csPDM. This is one of the directions of our future work for improving csPDM.
The derivations of csPDM and numerical examples used for testing are based on a single uncoupled population. Of course, what we derived in this study can be extended to the applications of simulating larger neural networks consisting of multiple coupled populations. However, for the correct use of csPDM, synaptic connections must be sparse, because a large number of synaptic connections can result in the violation of basic assumptions of PDM-Poisson and statistically uncorrelated inputs received by each neuron (Brunel, 2000;Nykamp and Tranchina, 2000;Huertas and Smith, 2006;Augustin et al., 2013). These basic premises can be violated when two neurons share many common pre-synaptic neurons, that is, when there are a huge number of connections between presynaptic and postsynaptic populations. We found that the number of synaptic connections depends on the size of ⌬V s . If ⌬V s is large, the number of connections has to be small because postsynaptic neurons have to receive statistically independent inputs. Postsynaptic neurons are more likely to receive exactly the same inputs if ⌬V s is large and the number of synaptic connections is huge.

Conclusions
When including many state variables, such as synaptic conductance variables for different types of receptors and synaptic facilitation and depression state variables, the original PDM becomes impractical because of the incredible increase in the computational load. In this study, this critical issue is solved by the probability density function method for Langevin equations with colored driving noise. The newly proposed method, termed csPDM, gives quantitatively accurate firing rate responses in both the steadystate and nonequilibrium regimes and highly qualified estimations of time-varying marginal density functions of membrane voltages. Importantly, csPDM is able to qualitatively manifest STP in an easy way almost without increasing computation demands. It appears that csPDM is generally applicable as a time-saving tool for modeling large-scale neural networks. We hope that our work will inspire further progress in the development of PDM and benefit computational and theoretical studies of synaptic dynamics in network dynamics.