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.
Significance Statement
Our study successfully solve an outstanding problem, how to reduce the dimension of population density equations when realistic synaptic dynamics is incorporated. With the newly proposed Fokker–Planck formalism, population density method (PDM) is conferred short-term plasticity (STP) properties and therefore becomes more widely applicable. As such, our method offers an opportunity to use the PDM to gain new insights into neural mechanisms of brain functions that are strongly dependent on STP synapses. This is the first step toward macroscopic description of large-scale neural network activities, reflected in some commonly used neurophysiological measurements, e.g., EEG, MEG, fMRI, and voltage-sensitive dye imaging data.
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 high-dimensional 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 dimension-reduction 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 synapse-associated 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 (Naud et al., 2008; 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 cs synaptic connections. In other words, each neuron has totally 
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
(1)where 
is the total synaptic current generated by the external synaptic inputs, C is the membrane capacitance, gl is the leak conductance, El is the leak reversal potential, κ is the sharpness factor, and VT is the threshold voltage. The spiking mechanism is the following: a spike is triggered at ts when V (ts) reaches a cutting voltage, Vc (in general, 
). Afterward, V is immediately reset to a resetting voltage, Vr, and, then, clamped for a refractory period, τref.
Synapse models
Specifically, 
is the sum of all postsynaptic currents mediated by various types of neurotransmitter receptors:
(2)
where gs(t) is the (collective) synaptic conductance for the input mediated by the sth type of receptors, and Es is the corresponding synaptic reversal potential. Equation 2 refers to conductance-based synaptic models.
The synaptic dynamics means the dynamics of the gs(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:
(3)
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 
. 
is the arrival time of the jth synaptic event, and 
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:
(4)
where ΔVs is a free parameter defined as the instant voltage modulation provoked by a single synaptic event when the membrane voltage is El just before the event arrives in cases where κ (in Eq. 1) is in the limit of 
and the synaptic dynamics has a zero-order kinetics, i.e., 
.
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
(5)
where the running parameters, us and xs, are governed by
(6)
(7)
at the level of a single synapse. The jump size of the synaptic conductance caused by a spike now depends on us and xs. us is the running utilization parameter and xs is the running fraction of available neurotransmitters. The dynamic value ranges from Us to 1 for us and from 0 to 1 for xs. Us refers to the base level of us. us = 1 means a presynaptic spike is allowed to use all available neurotransmitters, and xs = 1 means all neurotransmitters are available.
Briefly, the Equations 6, 7 depict a phenomenon, where when a synaptic spike arrives, us increases by an amount of 
, and, afterward, decays to its baseline level, Us, with the facilitation time constant, τs,f, and, meanwhile, xs decreases by an amount of xsus 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, Es, and Γs, each synaptic connection with the STP property is also characterized by three additional parameters: Us (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 (
and relatively large values of Us) to strong facilitation (
and small values of Us; 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 
(
) 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 cs = 200 for all receptor types.
In the MCS, the output population firing rate is calculated by
(8)
where 
is the total number of spikes produced by all neurons in the population within the time period of 
, and 
is a small time window. Δ t is always set to 1 ms in our simulation examples.
fdPDM
The assumptions and conditions, under which the fdPDM 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, 
, defined by
(9)
where 
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 
ranges from 0 to ∞. About the membrane voltage V, its upper bound is Vc, and its lower bound is given by 
. 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 
forms a (1+ m)-dimensional (in space) partial-integral differential equation:
(10)
where the voltage probability flux, i.e., the probability flux along the voltage-direction, JV, is given by
(11)
and 
is the conductance probability flux vector composed of m conductance probability fluxes. The conductance probability flux, 
, is given by
(12)
The population firing rate, r(t), is obtained by integrating the voltage probability flux across Vc over conductance, yielding
(13)
denotes the subdomain where 
, meaning that neurons at Vc intend to cross the boundary from below to generate spikes in this domain. As seen, for a given set of inputs, νs(t) (
), 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 Vlb 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 Vc, accounting for the generation of action potentials, re-enters the state space on Vr 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, 
in the Equation 10 is viewed as a deterministic parameter so that the master equation becomes a one-dimensional partial differential equation:
(14)
In this case, a steady-state firing rate for a given 
is obtained by using the mean-field model, given by
(15)
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 gs evolves isolatedly. For a large neuronal population, one can define a marginal conductance density function, 
, for gs by analogy with 
. Also, the integration of 
over a certain space on gs is the probability of finding neurons whose values of gs are located in that space. According to what derived by Ly (2013), the master equation governing the time evolution of 
is yielded as
(16)
where
(17)
The boundary conditions for solving Equation 16 are 
. Assuming that the values of all elements of 
are drawn from the marginal conductance density functions individually, MMFM estimates the time-varying firing rate, 
, with expected value of 
conditioned on 
, that is,
(18)
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 
at all possible points of 
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 gs.
In the diffusion limit, under the condition that the mean rate of synaptic input received by each neuron, 
, is sufficiently high and the jump size 
for each synaptic spike is small enough, the random synaptic conductance 
can be treated as the Ornstein–Uhlenbeck process (Uhlenbeck and Ornstein, 1930; Risken, 1996), whose dynamics is given by (Richardson, 2004)
(19)
where 
is the mean value of 
and 
is the standard deviation of 
. 
is a δ-correlated white-noise process of unit variance. Introducing a new variable 
and substituting for gs in the equation above, we obtain a new Ornstein–Uhlenbeck process for 
:
(20)
has zero mean and unit variance (Risken, 1996); and its autocorrelation function is
(21)
As a result, 
is a colored noise due to the exponential form of 
and finite synaptic time constant, i.e., 
. Replacing 
in the Equation 1 with 
, we obtain
(22)as the new dynamic equation for the membrane voltage, where 
, 
, and 
. So, as shown in Equation 22, the fluctuation of 
now is characterized by a Langevin equation with m independent colored noises 
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
(23)
for 
, where 
. Each 
is decomposed into a deterministic function or “mean-field velocity” 
and a stochastic fluctuation term 
. The authors proposed a quasi-Fokker–Planck equation:
(24)
as the master equation governing the temporal evolution of the joint probability density function of the system states, 
. The stochastic diffusion tensor, D, is calculated by
(25)
if fluctuation velocity components are exponentially autocorrelated and mutually uncorrelated, i.e.,
is the Jacobian of the mean-field velocity with components 
. Considering V, 
, and 
in Equation 22 as x, 
, and 
in Equation 23, respectively, we can yield the master equation for the marginal voltage density function, 
, by using Equation 24, as the following:
(26)
where
(27)
By using Equation 25, we obtain
(28)
(29)
where
(30)
and 
denotes the average membrane voltage across the population.
What remains in Equation 26 are equations for evaluating 
and 
. To derive the equation for evaluating 
, we multiply Equation 16 by gs and integrate gs from 0 to ∞:
(31)
Similarly, one can derive the equation for evaluating 
, yielding
(32)
For 
, we use the identity, 
.
In the csPDM, the population firing rate is directly estimated by the probability flux across Vc:
(33)
The average membrane voltage, 
, is computed by using its definition:
(34)
The boundary conditions for solving Equation 26 are assigned as:
They are similar to what used in the fdPDM except for the last boundary condition, which means that no neuron locates at Vc because the neurons whose voltages reach Vc are reset to Vr immediately.
When the STP property is included, we replace Γs and 
with
(35)
(36)
in Equation 31, Equation 32 for evaluating 
and 
, where 
and 
are mean values of us and xs across the neuronal population, respectively. Following mean-field equations proposed by Barak and Tsodyks (2007), we use the same equations for tracking 
and 
(i.e., Eq. 1 in Barak and Tsodyks, 2007):
(37)
(38)
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, 
, and 2m ordinary differential equations for tracking the first two statistical moments of all synapse-associated variables, 
, or 4m equations for statistical moments of 
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, 
, is introduced to rewrite the Equation 26 as followings:
(39)
(40)
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 
, we divide it into M meshes with an identical length as follows:
(41)
so that Vr and Es are exactly certain grid points. We denote the subspace 
(
) and its length 
. Using first-order polynomials as shape functions, the approximated values of 
and 
within Ik are defined by:
(42)
(43)
in which ϕ1 and ϕ2 are shape functions. 
and 
refers to the value of 
at 
from the left mesh Ik and from right mesh Ik+1, respectively. Due to the discontinuity at the interface of adjacent meshes, 
is thus possible. 
and 
are defined in the same way.
Element equations and numerical fluxes
Substituting 
and 
with approximated values 
and 
, respectively, multiplying shape functions and integrating with respect to v over the mesh Ik, we then obtain element equations for cell Ik, given by
(44)
(45)
in which the notation “^” means numerical fluxes at the interfaces between cells. We set 
and 
, and employ upwind fluxes for 
, yielding
(46)
According to the boundary conditions, we set 
, corresponding to 
, 
, respectively, and enforcedly assign 
as 
if 
. The firing rate r(t) is given by 
. We use the backward Euler method (Cheng and Shu, 2007) to solve element equations as well as Equations 31–38 to ensure numerical stability. Note that one needs to re-calculate element equations at each time step because 
and D are functions of time-dependent parameters, 
and 
, respectively.
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 
is set as 
. The ψ is the “minmod function” defined as:
(48)
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
(49)
(50)
to check the validity of diffusion approximation on synaptic dynamics. 
describes the average error ratio on the marginal conductance density function, which is estimated only over the interval 
. It is noted that the marginal conductance density function in the csPDM is characterized by a Gaussian distribution, here denoted by 
, whose mean and standard deviation are 
and 
. 
is the average error ratio on the marginal voltage density function.
The third quantity is
(51)
where 
. It is used for comparing performances of the fdPDM, MMFM, and csPDM. 
means the average error ratio on the population firing rate obtained from fdPDM, MMFM, or csPDM. 
is calculated from Equation 8. 
is calculated from Equation 13 for fdPDM, Equation 18 for MMFM, and Equation 33 for csPDM.
Results
We will present four simulation examples: (1) steady-state 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 
. Next, we start to show simulation results by first checking the validity condition of diffusion approximation of synaptic dynamics.
List of notation and symbols
List of values of model parameters
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 
and 
of a single neuronal population receiving only constant excitatory inputs with respect to different νs and 
. The synaptic time constant τs is set as 5 or 100 ms to match the time scale of AMPA-receptors or NMDA-receptors (Dayan and Abbott, 2001), and Es is set as 0 mV for both cases. Basically, small 
reflects the validity of diffusion approximation. As shown in Figure 1A, larger νs or smaller 
results in smaller 
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 
because of the decreased variation of gs (
at steady states). Three distributions of gs under different sets of (νs, 
) are shown in Figure 2A–C. As shown, the distribution of gs, 
, approaches a Gaussian distribution when 
is <0.2 (Fig. 2C). That is, the diffusion approximation of synaptic dynamics is valid if 
. Figure 1B shows 
. The amplitude of 
is proportional to that of 
, meaning that the error of csPDM in 
in large part comes from the failure of the diffusion approximation. But, surprisingly, the value of 
is much lower than that of 
under the same set of (νs, ΔVs). As seen in Figure 2A,D, under the case where 
mV, and τs = 5 ms, 
is only 0.174 whereas 
is as high as 0.653. However, the reason why 
is smaller than 
is unclear. Through the observation of 
(Fig. 2D–F), it is found that csPDM gives an adequately accurate estimation of 
comparable to MCS when 
is <0.2. To offer a quantitative and universal validity condition for diffusion approximation of synaptic dynamics, we also check the value of 
(i.e., the coefficient of variation of gs; Fig. 1C). In fact, it has been argued that the diffusion approximation is valid only when 
, meaning that the conductance mean 
should be much larger than the standard deviation 
(Richardson and Gerstner, 2005; Cai et al., 2012). It is however an impractical condition. Based on our observations of 
and 
, it can be conclusively said that the validity condition of diffusion approximation for the application of csPDM is 
, under which 
and 
are approximately <0.2. As a result, ΔVs should be <1.5 mV for the case of τs = 5 ms such that this criterion can be satisfied in the range of 
. In the following simulations, we choose ΔVs as 1 mV for the excitatory inputs.
Steady-state values of 
, 
and 
of a single population with respect to different νs and 
A, 
B, 
C, 
. Black dashed lines indicate the sets of (νs, ΔVs) where 
. Blank areas in A, B indicate the error ratio of >0.2.
Examples of steady-state marginal conductance density functions 
and marginal voltage density functions 
. A, D, 
and 
under 
and 
mV. B, E, 
and 
under 
and 
mV. C, F, 
and 
under 
and 
mV. As shown, the estimations of 
and 
by csPDM are similar with those by MCS if the input rate νs is increased. τs is set as 5 ms in all examples. Shaded area: MCS. Black thick line: csPDM.
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 ΔVs = 1 mV and Es = 0 mV. Figure 3 displays input-output curves computed by MSC, fdPDM, csPDM, and MMFM. Results show that csPDM gives an estimation of input-output 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 steady-state 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 demonstration that such overestimation or underestimation disappears when the synaptic time constant increases has not been reported.
Input-output curve comparison. The top panels in A, B show the input-output curves computed in four ways: MCS (black solid line with diamond markers), fdPDM (blue solid line with square markers), csPDM (orange dashed line with circle markers), and MMFM (pink dashed line with triangle markers). Synaptic time constant (τs) is 5 ms in A and is 100 ms in B. 
is the steady-state population firing rate as a function of a fixed input rate νs. The bottom panels in A, B show the errors of fdPDM, csPDM, and MMFM in 
from MCS. ΔVs = 1 and Es = 0 mV in this example.
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., 
. 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 
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, ΔVs, gl, κ and Vr, are chosen 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 ΔVs increases 
, meaning that the accuracy of csPDM is decreased. This is due to the fact that diffusion approximation becomes invalid when ΔVs is too large (as shown in Fig. 1). As expected, 
is less sensitive to the change of ΔVs than 
because ΔVs is arbitrary in the fdPDM. Surprisingly, although MMFM does not have limitations on ΔVs, in the case of τs = 5 ms, 
severely fluctuates over a range of 0.2 when changing ΔVs. Figure 5B–Dunravel that changing gl, κ and Vr almost does not affect the accuracy of 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 gl enlarges the error ratio 
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 gs-dimension in the master equation because, except for the requirement of more grid meshes along this dimension (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.
Dynamic population responses to a time-varying excitatory input. A, Time-varying input rate νs. B, top panel, Population responses r(t) in response to νs computed by MCS (shaded area), fdPDM (blue solid line with square markers), csPDM (orange dashed line with circle markers), and MMFM (pink dotted line with traiangle markers) under τs = 5 ms. Bottom panel, Corresponding errors of fdPDM, csPDM, and MMFM in r(t) from MCS. C, Similar to B except for τs = 100 ms. Average error ratios 
were computed over the interval from 600 and 1000 ms. 
, and 
are 0.058, 0.048, and 0.179 under τs = 5 ms, respectively. They are 0.042, 0.043, and 0.056 for τs = 100 ms. ΔVs = 1 and Es = 0 mV in this example.
Average error ratios of fdPDM, csPDM, and MMFM in population responses under different parameter sets. A, Subject to varying ΔVs. B, Subject to varying gl. C, Subject to varying κ. D, Subject to varying Vr. Blue solid line with square markers: fdPDM. Orange dashed line with circle markers: csPDM. Pink dotted line with traiangle markers: MMFM. 
means average error ratios corresponding to fdPDM, csPDM, and MMFM.
Comparison of computational time between fdPDM and csPDM. As shown, the computational speed of csPDM is almost 1000 times faster than fdPDM.
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, GABAA and GABAB (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 ΔVs = 1 mV, and inhibitory inputs mediated by GABAA-type receptors, which have τs = 10 ms and ΔVs = 0.25 mV, as well as GABAB-type receptors, which have τs = 100 ms and ΔVs = 0.25 mV. Each neuron is assumed to have 200 excitatory and 200 inhibitory synaptic connections.
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 GABAA and GABAB receptors share the common input rate in statistical 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 
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 
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 
are 0.065, 0.07, and 0.122, respectively. These results indicate that 
can be correctly captured by csPDM.
MCS and csPDM for a single population receiving excitatory and inhibitory inputs gated by three types of neurotransmitter receptors. A, The panels from top to bottom show (1) time-varying excitatory and inhibitory input rates (νe denotes the excitatory input rate, and νi denotes the inhibitory one); (2) raster plot of output spikes from 200 neurons in MCS; (3) the population firing rates r(t) computed by MCS (gray solid line) and csPDM (black dashed line); (4) the error in r(t) between MCS and csPMD; and (5) average membrane voltages across the population r(t) computed by MCS (gray thick solid line) and csPDM (black thick dashed line). Gray thin solid lines: voltage traces from 10 neurons in MSC. B, Three snapshots of the marginal density function 
at t = 1700 , 1800, and 1900 ms from left to right. Shaded area: MCS. Black thick line: 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, 
mv Es = 0 mV, τs,f = 700 ms, τs,r = 100, and Us = 0.05. They are the same except τs,f = 50 ms, τs,r = 300, and Us = 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 us and xs 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 
from 0.05 to ∼0.29. Meanwhile, the value of 
remains at a high level (
, 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 
roughly by 0.4, during which 
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 
between MSC and csPDM. In summary, the population response can inherit the STP property from synapses, and csPDM captures such property qualitatively.
STP population responses. A, Time-varying excitatory input rate. B, Population responses under short-term facilitation. The population response is facilitated during the sinusoidal input, resulting from the sharp increase in the running utilization parameter, 
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, 
. 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 
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, GABAA, GABAB, 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. 
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 (
) 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 
when gl 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, 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 103–104 (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^ calculus 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 large-eddy-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 
(τ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-fire-or-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 (Naud et al., 2008; 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 
, 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 
. This poor estimation is caused by a wrong assumption that us and xs are statistically independent for deriving mean-field equations for 
and 
(Tsodyks et al., 1998). In fact, us is independent of xs because xs does not appear in the Equation 6, but, by contrast, xs has a negative correlation with us. As a consequence, the mean-field equation for 
should include higher-order statistics of both us and xs. However, such an equation is unclosed. Before a closure method is available, the assumption of independence between us and xs is a compromised solution (Tsodyks et al., 1998). 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 ΔVs. If ΔVs 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 ΔVs 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 steady-state 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.
Footnotes
The authors declare no competing financial interests.
This work was supported by the National Cheng Kung University Hospital and funded by a grant (MOST 105-2314-B-006-080-MY3) from the Ministry of Science and Technology, Taiwan.
This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International license, which permits unrestricted use, distribution and reproduction in any medium provided that the original work is properly attributed.
References
In this issue
