The manner in which different distributions of synaptic weights onto cortical neurons shape their spiking activity remains open. To characterize a homogeneous neuronal population, we use the master equation for generalized leaky integrate-and-fire neurons with shot-noise synapses. We develop fast semi-analytic numerical methods to solve this equation for either current or conductance synapses, with and without synaptic depression. We show that its solutions match simulations of equivalent neuronal networks better than those of the Fokker-Planck equation and we compute bounds on the network response to non-instantaneous synapses. We apply these methods to study different synaptic weight distributions in feed-forward networks. We characterize the synaptic amplitude distributions using a set of measures, called tail weight numbers, designed to quantify the preponderance of very strong synapses. Even if synaptic amplitude distributions are equated for both the total current and average synaptic weight, distributions with sparse but strong synapses produce higher responses for small inputs, leading to a larger operating range. Furthermore, despite their small number, such synapses enable the network to respond faster and with more stability in the face of external fluctuations.
Experiments analyzing the distribution of synaptic weights impinging onto neurons typically observe low-amplitude peaks with only few, large-amplitude excitatory (EPSPs) or inhibitory post-synaptic potentials (IPSPs) [1]–[14]. These have been fitted by lognormal [1], truncated Gaussian [10], [11] or highly skewed non-Gaussian distributions [4], [6], [9], [13]. This raises the question of the functional role of such relatively rare but powerful synaptic inputs. The functional implications of such strong synapses can be very significant [15]. Circuits of the Mind [16] proposes a powerful computational brain architecture (the neuroidal model) to explain the brain' s remarkable flexibility to quickly memorize new events and associate them with previous stored ones. It is based on a very small fraction of powerful excitatory synapses. Furthermore, a recent study [17] of strong but sparse synapses, combined with weak and probabilistic synaptic amplitude distributions provided both computational justification as well as empirical support for the role of these rare yet powerful synaptic events in supporting low-frequency, spontaneous firing in neuronal networks at rest. This question has become more acute following several reports that individual cortical pyramidal neurons from human tissue recovered during surgery are sufficiently powerful to drive other neurons by themselves [18], [19], unlike equivalent cells in rodent cortex. The relation between probabilistic synaptic weight distribution and population dynamics can be studied using simulations [17]. However, the large parameter space to be explored and the need to repeat simulations many times makes this an impractical first method to apply, in particular when modeling the activity in large regions or even the entire mammalian brain. An alternative is the analytical Fokker-Planck method, based on a continuous stochastic process (a brief review of stochastic processes following [20] and [21]–[23] has been provided in Text (S1) ) that models the dynamics of homogeneous neuronal populations with a single partial differential equation (PDE). It is able to quickly explore parameter space and provide analytical insights [24] and characterizes any distribution of synaptic weights by just quantities-the drift and diffusion terms in the equation, corresponding to what the mean and variance of the membrane potential would be if the neuron did not have a threshold. But this method cannot reproduce dynamics in the presence of large synapses. Models based on jump stochastic processes [25], [26] treat synaptic input as being composed of pulses with finite amplitudes. The population dynamics in the presence of such synapses have been investigated using analytical and numerical techniques [27]–[38]. Richardson and Swarbrick [39] characterized analytically an integral formulation in the case of exponentially distributed stochastic jump processes. We here present a fast, semi-analytical approach to study the integral formulation of stochastic jump processes with arbitrary distributions. We apply this method to study the role of synaptic distributions on neuronal population dynamics.
The equation for the membrane potential relative to rest for a generalized leaky integrate and fire (gLIF) neuron with normalized capacitance is, (1) is a random variable characterizing the shot-noise synaptic current. It takes the value with probability and with probability. is the input synaptic event rate so that is the mean number of inputs in a time. represents the distribution of synaptic weights, with and and being the minimum and maximum synaptic weights respectively. When the synaptic input is sufficient to cause the membrane potential to exceed a threshold value, it is reset so that (2) This reset implementation is used to account for the shot-noise nature of the synaptic input. For nearly instantaneous synapses, represents the peak of the post-synaptic potential (either EPSPs or IPSPs) (for non-instantaneous synapses, see Methods: Non-instantaneous synapses). represents the sum of all non-synaptic currents, which can be voltage-dependent but not explicitly time-dependent. For the standard LIF neuron, , where is the membrane time constant. For an exponential integrate-and-fire neuron, . is the spike detection threshold and is the slope factor. The resulting equation for the probability of a neuron to have a voltage in at time -the master equation-is an integro-partial differential equation with displacement (DiPDE), (3) with is the unit step function and the threshold membrane potential. The first term on the RHS in equation (3) incorporates the drift due to non-synaptic currents. The second term removes the probability for neurons which previously were at potential and received a synaptic input. The third term adds the probability that a neuron away in potential receives a synaptic input such that its potential is changed to. The last term represents a probability current injection of the neuron which previously spiked. includes the effect of any excess synaptic input above the threshold at which spiking occurs, due to large super-threshold synaptic events. To account for the effect of the excess input with instantaneous synapses, the probability current is injected between the resting potential and, where represents the membrane potential that would be reached due to the excess input. This serves the purpose of ‘remembering’ the excess input, whose effect would have been held by the synaptic variables in the case of slow synapses, without losing it on resetting the membrane potential after spiking. The output firing rate is given by (4) In the Fokker-Planck formalism, is part of a boundary condition. If, then there is no probability current through threshold due to continuous processes for this system. The alternative case is discussed in (Methods: Boundary conditions). Since we include non-infinitesimal synaptic inputs in an infinitesimal time interval, it is possible for a neuron to cross the threshold by a large amount. In equation (3), this additional depolarization is accounted for in the term as outlined above. The case in which this additional depolarization is ignored is treated in (Methods: Boundary conditions). Simpifications obtained with -function distribution of synaptic weights are described in (Methods: Simplifications). It is also possible to include the effect of synaptic delays and distributions of synaptic delays within the DiPDE formalism, as discussed in (Methods: Non-instantaneous synapses). The stationary solution for equation (3) can be obtained as the solution to the following equation, (5) where is the stationary probability distribution for the membrane potential. Although the model can include both excitatory and inhibitory connections, all the results presented below, except for those in Figure S6 c) and d), are for feed-forward networks with excitatory connections. Intuitively, this model quickly diverges from the Fokker-Planck formalism when the synaptic strength is large. Consider the hypothetical case of a neuron starting at rest which has all the synapses equal and large. After a short time step, the Fokker-Planck formalism produces a narrow Gaussian distribution in membrane potential near rest, while the DiPDE formalism yields a membrane potential distribution which is a sum of two scaled delta functions: a large one at rest, and a small one at the synaptic weight. Over time, the Fokker-Planck equation converges to a single broader Gaussian, while the DiPDE formalism leads to a larger coefficient for the delta function at the synaptic weight value. A generalization to conductance-based synapses is presented in (Methods: Conductance-based synapses), while a generalization to the case of exponential integrate-and-fire neurons [40] is presented in (Methods: Exponential integrate-and-fire neurons). We solve equation (3) without the displacement terms using the method of characteristics. The characteristic equations can be solved to obtain a non-uniform discretization of the membrane potential. For the standard LIF neuron, the characteristic equations are solved analytically. We then numerically add the effects of the displacement terms at every time step (see Methods: Numerical Solutions). This semi-analytic technique has the advantage of reducing errors due to numerical diffusion at each time-step. The solution of the DiPDE equation (3) is in good agreement with simulations of 10,000 leaky integrate-and-fire neurons for both low frequency, large amplitude (Figure (1) ) and for high frequency, small amplitude distributions (Figure (S1) ). Differences between continuous and discontinuous stochastic processes can be seen in the transient behavior of the probability distribution of membrane potential in the top panel (Figure 1a, b, c) and are statistically significant. p-values for differences between the sub-threshold steady-state membrane potential distributions have been provided in Text (S4), Table (S2) and Table (S3). The steady state firing rate obtained from DiPDE is Hz and from simulations is Hz. The small discrepancy is caused in part by simulating synapses which are not instantaneous but rather have a time constant of ms (for more details see (Methods: Non-instantaneous synapses) ). Since the stochastic input is Poisson-distributed, expected 95% intervals for spike counts can be directly computed from DiPDE (see Methods: Expected 95% intervals for spike counts and Figure (1f) ). The transient time to firing, defined as the time taken to reach 10% of the equilibrium firing rate, is ms and provides a good estimate of how quickly the neuronal population responds to a given input. By comparison, the Fokker-Planck formalism results in a lower steady state firing rate of Hz and a slower transient time ms. The higher number of neurons closer to resting potential at equilibrium in Figure (1e) reflects the nature of the ‘jump’ stochastic process. For higher-frequency, low-amplitude inputs, all three methods converge to the same results (Text (S2) and Figure (S1) ). Thus, the formalism presented in equation (3) is in good agreement with simulations for instantaneous synaptic input and does not depend on the choice of time-step for the numerical solution (Figure (S2) ). For non-instantaneous synapses, upper and lower bounds on the steady state output firing rates can be obtained (Methods: Non-instantaneous synapses, Figure (S4) and Figure (S5) ). The DiPDE implementation also matches equivalent simulations for conductance-based synapses (Figure (1g), Figure (S3) ) and exponential integrate and fire neurons (Figure (1h) ). We used the DiPDE formalism to investigate the effect of generic synaptic weight distribution on the steady-state, subthreshold membrane potential distribution and output firing rates in feed-forward networks. Gaussian synaptic weight distributions with different mean weights whose input firing rates are adjusted to produce the same average synaptic current, result in different transient and steady-state firing rates as well as different equilibrium voltage distributions (top row in Figure (2), balanced excitation/inhibition Figure (S6) ). We then examined six different distributions tuned to have both the same mean (1 mV) synaptic weight and input rate (1,000 Hz), so that the average synaptic current was the same (see Methods: Matched synaptic distributions for the exact forms of the distributions and their respective variances). We find that the average synaptic current is not sufficient to accurately determine either the output firing rate or membrane voltage distributions (bottom row in Figure (2) ). We also tested distributions that were matched to have input synaptic current with the same mean and variance (Figure (S7) ). The undulating voltage distribution for -function input (Figure (2e) ), reflects the nature of the ‘jump’ stochastic process. Heavier-tailed distributions generate faster transient responses to inputs and monotonically reach steady-state output firing rates. For example, the power-law distribution leads to an equilibrium firing rate Hz and a transient time of ms, while the Gaussian distribution results in Hz and ms. Numerical results for all simulations are provided in Table (1). Even changing a small fraction of the synaptic weights can have a significant effect. For example, the -function distribution converges to the lowest steady state firing rate of Hz (for computation of 95% confidence-intervals, see Table (1) ). In contrast, the bimodal distribution which differs from the -function by only 3. 4% of synapses having a much higher weight, results in the highest steady state firing rate of Hz. The tail of the synaptic weight distribution has an even larger effect on the transient times starting from rest - it is 16. 2 ms for the -function, but only 2. 6 ms for the bimodal distribution. Since the network we study is a feed-forward network, the mean synaptic delay results in a simple time translation of the responses. However, the overshoot of steady state seen in Figure (2f) decreases as the variance in the distribution of synaptic delays increases (see Methods: Synaptic delays and Figure (S8) ) relative to the membrane time constant. To analyze what characteristic of the synaptic distribution is most important for the fast response observed for our heavy-tailed distributions, we generated more than 1000 random distributions matched to have the same mean synaptic weight and input rate (Methods: Tail weight numbers) and analyzed the responsiveness of a neuronal population (time to a fraction of the equilibrium firing rate) as a function of the moments of the synaptic distribution (Table (2) and top row of Figure (S10) ). Higher order moments explain part of the variance observed in the responsiveness. However, a larger fraction of the variance can be explained by introducing a set of measures specifically designed to quantify the number of strong synapses (Table (2) and bottom row of Figure (S10); see also Figure (S9) for our six chosen distributions). These are the tail weight numbers, the density above threshold of convolutions of the effective synaptic weight distribution: (6) where is related to the synaptic weight distribution by, (7) and represents a convolution. represents a Poisson process with mean and events occurring in a time-step. By definition, , , and so on. thus represents the average distribution of depolarization of a single neuron, when each neuron in the population receives excitatory inputs on average. then represents the fraction of neurons that spike in a neuronal population starting at rest when each neuron in the population receives excitatory inputs on average. Having established a measure of population activity when all neurons start at rest, we examined the dynamics resulting from an equilibrium different from rest (see Methods: Input-output curves). Mathematically, this amounts to analyzing the effect of synaptic weight distribution on the transient dynamics from one stationary solution to another stationary solution, when the input synaptic event rate in equation (3) is instantaneously changed by a constant from. The different synaptic distributions result in an overshoot of the population' s eventual equilibrium firing rate. The overshoot becomes smaller for heavier-tailed distributions. For the power-law, there is no overshoot present. The amount of overshoot is a measure of the stability of the neuronal population to sudden changes in input firing, given a distribution of synaptic weights. The weight of the tail in the synaptic distribution correlates with the stability of the system (Figure (3) and Figure (S12) ). The DiPDE formalism also enables insights into the influence of short-term synaptic depression (STSD), which is known to play a key role in neural network homeostasis and in the generation of multiple network states [41] via a Tsodyks-Markram mechanism [42]. With the inclusion of synaptic depression (see Methods: Implementation of short-term synaptic depression), the output firing rates for all six chosen synaptic weight distributions begin to saturate when the effective input synaptic events per second (integrated over all the synapses impinging onto a neuron) are around Hz (Figure (3f) ). Heavier-tailed synaptic weight distributions lead to a higher dynamic range (Table (S4) ). In the preceding analysis, we saw how perturbations of the input event rate in a population with different synaptic weight distributions affect the entire time-course of the population' s evolution from an initial to a final equlibrium. In contrast, we also analyzed how fluctuations in inputs (see Methods: Fluctuation Analysis) affect the instantaneous change in equilibrium firing rate of neuronal populations with different synaptic weight distributions (Figure (4), Figure (S13) and Figure (S14) ). To quantify the instantaneous change in the equilibrium firing rate, we make use of a closely-related variant of the tail weight numbers defined in equation (6). Instead of starting with an initial probability distribution for the membrane potential corresponding to all neurons at rest, we use the stationary probability distribution obtained from the solution to equation (5) with input rate to define, (8) where is as defined in equation (7). thus represents the fraction of neurons that spike in a neuronal population starting from the stationary distribution, when each neuron in the population receives additional inputs on average. This means that, effectively in a time, the total number of inputs becomes. Using, we can define the relative excitability of a neuronal population with a given synaptic weight distribution as (9) This is a measure of susceptibility of the neuronal population. If is the increase in input fluctuations for a population with additional inputs per neuron on average, is the corresponding increase in the output firing rate of this population. A relative excitability implies that the additional instantaneous response to an external input is independent of other inputs received at the same time. Starting from the equilibrium of Figure (2), for all synaptic weight distributions, the relative excitability initially rises (Figure (4a) ). In this case, it can be partially compensated by synaptic depression (Figure (4b) ). Table (S5) lists the maximum relative excitabilites for all the distributions, without and with synaptic depression. Results for an equilibrium obtained with a higher input firing rate have been presented in Figure (S14). Heavier-tailed distributions of synaptic weights cause smaller changes in relative excitability of the neuronal population. A few strong synapses can exert an undue large influence on the mean response of the population. Consider the unimodal -function synaptic distribution, whose weight is set at 1 mV and the bimodal distribution, where 96. 6% of the synapses have 0. 5 mV weight, while 3. 4% of its synapses are much larger at 15 mV (an even more extreme case with the two distributions differing in only 0. 1% of their synapses is presented in Figure (S11) ). Even though both the average synaptic input and the total current are the same for both distributions, their steady-state firing rates are 19. 6 and 28. 7 Hz, their transient times are 16. 2 and 2. 6 ms and their dynamic ranges (see Table (S4) ) are 10. 22 and 31. 83 respectively. More dramatic is the manner in which these two distributions react to synaptic fluctuations. When each neuron in the population receives additional inputs on average the relative excitability for the -function distribution is, while that for the bimodal distribution is (see Figure (4a) ).
Our study provides an advance on two fronts, computationally and neurobiologically. First, we developed and validated a semi-analytic method to model the sub-threshold membrane potential probability distribution and the firing rate of homogeneous neuronal populations with finite synaptic inputs. Second, we apply this method to explore the effect of varying synaptic weight distributions on equilibrium and transient population characteristics. From a methodological standpoint, the DiPDE formalism reproduces population behavior from aggregate simulations of identical point neurons, without the need to run the large-scale simulations themselves with the attendant computational costs (see Text (S3) and Table (S1) ). The Fokker-Planck equation works best in the regime where the synaptic distribution and input firing rates approach a continuous process. It does not accurately model the response to even rare high-amplitude (‘jump’) synaptic inputs. By contrast, the ability of the DiPDE formalism to model both jump and continuous processes makes it a powerful framework for modeling single as well as multiple interacting neuronal populations. This population statistics method makes possible the modeling and characterization of the membrane potential distribution and the spiking statistics of very large networks of distinct but homogeneous populations of hundreds of neuronal cell types in numerous brain regions. This would be relevant, for instance, when modeling the resting state activity at the cellular level throughout the awake or the sleeping brain. Second, in order to understand the biological significance of powerful but rare EPSPs [7], [18], [19], we used the DiPDE formalism to compare distributions which contain some very strong synapses against distributions with a preponderance of small EPSPs (but that are matched in their mean and variances). We will refer to the distributions as heavy-tailed, but a more quantitative description is presented in Methods: Tail weight numbers. We showed that power law and related heavy-tail distributions lead to faster transient behavior than non-heavy-tailed distributions (Figure (2f) and Figure (3f) ). Second, heavy-tailed synaptic distributions lead to a higher dynamical range (Figure (3f) ) than non-heavy-tailed distributions. Third, heavy-tailed distributions are much less sensitive to random fluctuations in synaptic activity (Figure (4) ). All three properties associated with such heavy-tailed distributions can be functionally advantageous compared to matched synaptic input with no such large synapses. For example, faster transient responses would be desirable when a potential prey needs to detect the presence of a predator and escape. Higher dynamical ranges would be useful when sensory systems need to respond external stimuli over a large range, without the need to have multiple networks designed to respond over smaller ranges. The lesser sensitivity of heavy-tailed distributions to random fluctuations in synaptic activity can be useful when the response to a stimulus needs to be independent of its context. Teramae and colleagues [17] recently published a joint computational and physiological investigation into the role of strong but sparse excitatory EPSPs, superimposed onto a very large pool of weak EPSPs and discovered the critical role that the former play in generating and sustaining long-term, low-frequency spontaneous firing activity in mixed excitatory-inhibitory neuronal networks in the absence of sustained external input or NMDA synapses. Such states allow the neurons to often reside one synaptic input away from the threshold, thereby leading to high correlation between the pre-synaptic and post-synaptic neurons of the strong synapses. Such investigations highlight the need to focus the attention of electrophysiologists studying synaptic transmission in vivo onto the critical role that such rare, “Black Swan”-like, events can play in the day-to-day life of the brain. Given the large biological and instrumental noise present in synaptic measurements, in particular under in vivo conditions, distinguishing between these distributions in practice would not be easy as it would require recordings of long duration to collect the relevant statistics. For experimental purposes, the numbers of independent recordings of sub-threshold steady state membrane potential required to differentiate between different synaptic weight distributions have been provided in Table (S2) and Table (S3). Yet as shown here, large but rare excitatory synaptic inputs can exert undue influence on population dynamics and robustness. These conclusions are very important for the emerging field of connectomics: a weighted graph description for neuronal networks in which each node represents a homogeneous population and each connection is characterized by the mean postsynaptic current is not sufficient. The network has to be either refined to individual neurons or expanded to include knowledge of the entire distribution of postsynaptic currents. The method further developed here facilitates the modeling and characterization of the membrane potential distribution and the spiking statistics of very large networks of distinct but homogeneous populations of hundreds of neuronal cell types in numerous brain regions.
A simplified case can be obtained if all EPSPs are assumed to have the same value; that is, if. Equation (3) becomes: (10) This is a first order partial differential equation with displacement (DiPDE). For numerical solutions of equation (10), the matrix obtained by discretizing is as sparse as the matrix needed to solve the Fokker-Planck equation. If, then equation (5) reduces to a delay differential equation (DDE) (11) In equation (3), if the additional depolarization after the neuron crossed threshold is neglected, the equation becomes (12) where is the two-dimensional Dirac -function and corresponds to the resting potential. If, then there is an additional contribution from the continuous processes to the probability flux through threshold. The output firing rate is now given by, (13) Numerical simulations, against which the Fokker-Planck (for current-based synapses) and DiPDE formalisms were compared, were performed by invoking the NEST simulator after writing the code in PyNN. For current-based synapses, simulations were performed for a population of independent and identical leaky integrate-and-fire (LIF) neurons with decaying exponential post-synaptic current. With the resting and reset potential being equal and denoted by, the neuron parameters were chosen to be mV, mV, membrane time constant ms, membrane capacity nF, refractory period ms, input resistance M and decay time ms for excitatory synapses. All neurons in the population were supplied with a Hz Poisson excitatory input of amplitude nA (although this is a huge current, with our choice of parameters, the effective charge deposited on the neuron is given by pC) and the membrane potential was recorded over a simulation time of ms with a time-step of ms. The recorded voltages were then analyzed to obtain the probability distributions shown in Figure (1). The recorded spike-times were binned with an interval of ms and the resultant output firing rate was obtained as shown in Figure (1d). For simulation results presented in Figure (S1), we used Hz Poisson excitatory input of amplitude nA. For conductance-based synapses, simulations were performed for a population of independent and identical leaky integrate-and-fire (LIF) neurons with decaying exponential post-synaptic conductance. The neuron parameters were chosen to be mV, mV, excitatory reversal potential mV, membrane time constant ms, membrane capacity nF, refractory period ms and decay time ms for excitatory synapses. All neurons in the population were excited by external Poisson input with input firing rate Hz and peak synaptic conductance µS and their membrane potential was recorded over a simulation time of ms with a time-step of ms. The recorded voltages were then analyzed to obtain the probability distributions shown in Figure (S3). The recorded spike-times were binned with an interval of ms and the resultant output firing rate was obtained as shown in Figure (1g) and Figure (S3c). Simulations were also performed for a population of independent and identical exponential integrate-and-fire (EIF) neurons with conductance-based synapses without adaptation. The neuron parameters were chosen to be mV, mV, mV, slope-factor mV, excitatory reversal potential mV, membrane time constant ms, membrane capacity nF, refractory period ms and decay time ms for excitatory synapses. The adaptation parameters were chosen to be ms, and, such that adaptation was absent. All neurons in the population were excited by external Poisson input with input firing rate Hz and peak synaptic conductance µS and their membrane potential was recorded over a simulation time of ms with a time-step of ms. The recorded spike-times were binned with an interval of ms and the resultant output firing rate was obtained as shown in Figure (1h). To solve the evolution equation (3) for the probability density, we first solve the advection (leak) portion of the equation, and then include the synaptic terms to calculate the overall time derivative. We then march the solution for forward in time using an explicit first-order-in-time scheme with constant time step. Based on the general explicit scheme criterion, we use a sufficiently small time step to ensure stability. This scheme is conservative; the integral of the probability distribution is not affected by numerical errors. For the standard LIF neuron with membrane time-constant, the leak term is linear with. The ODE has an analytical solution and the method of characteristics provides the full solution for the advection term. In order to ensure an exact implementation of the leak term, we use a geometric binning scheme for the membrane potential, with the bin ratio determined by the product of leak and time-step. The bin-edges are determined, starting from threshold membrane potential, using (14) A total of bin-edges are generated until the first bin between the resting potential and is at least as small or smaller than the first bin generated from. Mathematically, the lower bound on the number of bin-edges is given by the condition (15) since we have chosen and the second equality follows from equation (14). Increasing the numbers of bins beyond this lower bound increases the accuracy, but also increases the computational cost involved (see Text (S3) ). For more general forms of the leak term (for e. g, the exponential integrate-and-fire (EIF) neuron; see Methods: Exponential integrate-and-fire neurons) the ODE can be solved numerically starting from the threshold membrane potential in order to generate the bins. At each time-step, the probability distribution for membrane potential (which is a vector) is first evolved with the leak term and then the synaptic input. With our non-uniform binning scheme, the evolution of the probability distribution for the membrane potential due to the leak reduces to a single-index shift towards the resting potential and reduces the error due to numerical diffusion at each time step. To implement the effect of instantaneous synaptic input, including the effect of any excess input above the threshold the input distribution of synaptic weights is first used to construct a transition matrix. The additional columns per row in are used to keep track of the effect of excess synaptic input above. If this additional depolarization due to super-threshold inputs is ignored and the membrane potential is reset to zero upon exceeding, then is just a matrix. This is then used to construct an effective transition matrix (16) where (17) represents the probability for the neuronal population to receive inputs in a time from an external homogeneous Poisson process with mean input rate. The Poisson process is truncated at a sufficiently high value of such that. For inhomogeneous Poisson inputs, an effective transition matrix would have to computed at each time step, depending on the input rate. Each off-diagonal element in the effective transition matrix represents the proportion of neurons in the population receiving a specific synaptic input. The probability distribution for the membrane potential at each time-step is then multiplied by this transition matrix to generate the new probability distribution after synaptic input. The case for non-instantaneous synapses has been discussed below in Methods: Non-instantaneous synapses. In order to compare results with simulations, we need to match the total synaptic input to the neuronal population with that from simulations. The DiPDE formulation is exact when implementing instantaneous synapses. The synaptic weight is equated with the maximum depolarization of membrane potential (peak EPSP) achieved after input from an exponentially decaying current-based synapse, (18) This value is reached after a time (19) For our choice of simulation parameters outlined above mV provides good agreement with simulations as shown in Figure (1). The time-step was chosen to be ms to match with simulations. Another constraint imposed by the nature of the numerical solution is that the discretization of the membrane potential means that it is impossible to have the initial probability distribution to be a sharp -function at mV. The initial probability distribution is therefore spread uniformly over the width of the first voltage bin in both DiPDE and simulations. The underlying stochastic process in equation (3) assumes Poisson-distributed inputs. If an equivalent numerical simulation with input firing rate involves a finite population of neurons and the resultant spikes are binned in time intervals, then there are inputs per bin, with and. The relative variance in the number of inputs per bin is. Expected 95% intervals for output spike counts can then be obtained as follows. At each time step, the solution obtained from equation (3) for a given distribution of synaptic weights with a given input rate is additionally subjected to input rates and the corresponding output rates are calculated. These give the expected 95% intervals to the mean output firing rate. The 95% intervals for the expected spike counts per time-bin for a population neurons are in the interval between. Figure (1f) shows the output firing rates corresponding to the expected 95% intervals for results presented in Figure (1d). For a population of neurons and a bin size of ms, the 95% interval for expected output spike counts per bin obtained from DiPDE is (90. 3,99. 6) with corresponding output firing rates of Hz and Hz respectively. For simulations, the 95% interval for spike counts corresponds to (90. 5,104. 5). For non-instantaneous synapses, upper and lower bounds on the output firing rate can be obtained. The exponentially decaying current-based synapses used here in simulations (see (Materials and Methods: Current-based synapses) ) result in an exact EPSP given by for synaptic input at time. is the unit step function. If is the total charge deposited by an instantaneous, excitatory current-based synapse on a neuron with normalized capacitance at time, then the corresponding EPSP is given by, is the synaptic weight used in equation (3). Setting corresponds to equating the total charge, while setting corresponds to equating the maximum depolarization. Figure (S4) shows the output firing rates obtained from DiPDE after equating the total charge, the maximum depolarization, an intermediate estimate obtained by setting and that obtained by setting, along with the corresponding numerical simulation. A lower bound on the output firing rate can obtained by equating the maximum depolarization, with corresponding to. Equating the total charge with corresponding to provides an estimate of the output firing rate. An upper bound can be obtained by taking with corresponding to. Figure (S5) shows a plot of the different EPSPs obtained with mV, ms and ms for synaptic input at time. For, a formal proof for the upper and lower bounds can be provided as follows. Let, and represent the EPSPs due to synaptic input at time. The maximum depolarization attained due to inputs can be represented as: (20) For any, . Therefore, for all times after the synaptic input, so that. Hence, using in equation (3) an upper bound for the output firing rate is obtained. For the lower bound, consider the difference between and. For all, since, . For, . Using this, it is straightforward to show that when. Thus for all times and hence. The time translation by does not affect the result since' s have been defined to be the maximum depolarization over all times. For current-based synapses, the change in membrane potential resulting from synaptic input is independent of the initial membrane potential. For a conductance-based synapse however, this change is proportional to the difference between the initial membrane potential and the reversal potential for the channels present in a particular synapse. Equation (3) becomes, (21) Because of the additional dependence on the membrane voltage and synaptic reversal potential, modeling an instantaneous change in voltage through a conductance-based synapse is not as straightforward as with a current-based synapse. To model this type of synapse efficiently, we introduce the non-dimensional quantity as in equation (21) above, which represents the instantaneous voltage change due to synaptic activation as a fraction of the charge needed to reach the reversal potential from a given membrane potential. For example, indicates that a single synaptic event would shift the neuron from its current membrane potential to halfway between the current voltage and the synaptic reversal potential. Mathematically, if synaptic events have already occurred, then the change in membrane potential induced by the -th synaptic input can be represented by, (22) For such synaptic events, it can be verified that this reduces to: (23) This non-dimensionalization of synaptic weight for conductance-based synapses thus automatically rescales all weights to lie between zero and 1, and provides a simple way to model these synapses as generating instantaneous changes in membrane voltage. For our choice of simulation parameters, the maximum depolarization achieved by a neuron starting from rest due to synaptic input with and mV is mV. This is obtained from the numerical solution of the equation for a LIF neuron with exponentially decaying conductance-based synapses (with time-step chosen to be 0. 1 ms to match with simulations). Figure (S3) shows that numerical simulations (see Methods: Simulations for parameters used in simulations) are in good agreement with numerical solutions of DiPDE for conductance-based synapses. The equilibrium output firing rate is Hz and the transient time to firing is ms. For neurons and a bin size of ms, the 95% confidence interval for the spike counts per bin is (43. 0,47. 6) corresponding to output firing rates of Hz and Hz. For multiple neuronal populations connected to each other, one can generalize equation (21) to a system of equations, (24) while defining two synaptic weight matrices, The DiPDE formalism (equation 3) can be generalized to other types of neurons such as for e. g - the exponential integrate-and-fire (EIF) neuron. The membrane potential dynamics for the EIF neuron with excitatory conductance-based synapses is governed by (25) where the first two terms on the RHS contribute to the leak and the third term corresponds to the synaptic input. represents the leak reversal potential, is the spike detection threshold, is the slope-factor, is the membrane capacitance, is the membrane time-constant and is the excitatory reversal potential. In the absence of synaptic input, the leak portion can be solved numerically to obtain the discretized non-uniform bins for the membrane potential. For the exponential integrate-and-fire (EIF) neuron, the numerical solution of the equation with exponentially decaying conductance-based synapses gives. With mV, this means that the maximum depolarization achieved by a neuron starting from rest is mV (with time-step chosen to be 0. 05 ms to match with simulations). Figure (1h) shows that the output firing rate obtained from numerical simulations (see Methods: Simulations for parameters used in simulations) is in good agreement with that obtained from the numerical solution of DiPDE. The equilibrium output firing rate is Hz and the transient time to firing is ms. For neurons and a bin size of ms, the 95% confidence interval for the spike counts per bin is (48. 5,53. 1) corresponding to output firing rates of Hz and Hz. The Fokker-Planck equation for current-based synapses was solved using an explicit Forward-Time Centered-Space (FTCS) scheme. As in the DiPDE, the leak term was solved analytically, although the voltage discretization was kept uniform in order to be compatible with the standard Centered-Space discretization used for the diffusion term. To ensure that the discretization was sufficiently fine (especially since FTCS is only first-order accurate), we compared the equilibrium firing rate generated by this scheme against the analytical solution [43] given by (26) where is the analytical firing rate, is the refractory period and is the membrane time constant. and are parameters related to the drift and diffusion coefficients in the Fokker-Planck equation, and are given by (27) where and are excitatory and inhibitory synaptic weights respectively, while and are the corresponding input firing rates. The effect of various matched synaptic weight distributions on the population dynamics in feed-forward networks was investigated using the DiPDE formalism. The top panel of Figure (2) showed that two self-similar Gaussian distributions matched to produce the same average synaptic current resulted in different steady-state output firing rates and sub-threshold membrane potential distributions. In the regime when the firing is primarily driven by drift, the difference in output firing rates is small as shown in the top panel of Figure (S6). The higher-mean Gaussians result in an steady state output firing rate of Hz and Hz, while the lower mean Gaussians lead to Hz and Hz. If the firing is driven exclusively by variations in input (Gaussian distributions with balanced excitation and inhibition), the differences in output firing rates are large as shown in the bottom panel of Figure (S6). For the Gaussians with 6 mV mean and either 1 mV or 2 mV standard deviations, the steady state output firing rates are Hz and Hz respectively. The steady state output firing rates for the Gaussians with 3 mV mean and either 0. 5 mV or 1 mV standard deviations are Hz and Hz respectively. These results imply that population response is determined not only by the total current, but also by the mean synaptic weights. For the bottom row of Figure (2), six different synaptic weight distributions were tuned to have the same mean (1 mV). Using the same input rate (1,000 Hz) for the all the distributions, the average synaptic current was the same. The corresponding expressions for these distributions normalized between and are: The' s represent the normalization constants. The corresponding zero-centered second moments are: Delta (1. 0 mV2), Gaussian (1. 3 mV2), Exponential (1. 8 mV2), Lognormal (2. 2 mV2), Bimodal (8. 0 mV2) and Power-law (8. 7 mV2). The firing rates and transient times corresponding to distributions for which the mean input current was matched (Figure (2) is provided in Table (1). An alternate way to match synaptic weight distributions is to match the drift and diffusion, corresponding to what the mean and variance of the membrane potential would be if the neuron did not have a threshold, while allowing the input rates to vary. In the Fokker-Planck formalism, such distributions matched for drift and diffusion would have led to the same results. Figure (S7) shows the differences in output firing rates and equilibrium membrane potential distributions obtained from eq. (3) for synaptic weight distributions matched for drift and diffusion. The input rates are: Delta (555. 5 Hz), Gaussian (745. 98 Hz), Exponential (952. 38 Hz), Lognormal (1,062. 15 Hz), inverse power-law (4,166. 72 Hz) and Bimodal (1,825. 33 Hz). The corresponding expressions for these distributions normalized between and are: The' s represent the normalization constants. The heavier-tailed distributions still lead to faster transients, however the steady-state firing rates now decrease with increasing heaviness of the tail in the distributions. This is in contrast to the case when the input currents were matched (Figure (2). In order of increasing heaviness of the tail distributions, the equilibrium output firing rates and the transient firing times for the different distributions are provided in Table (1). We also implement synaptic delays within the DiPDE formalism. This is done by a using a queue to store the output firing rate, which is then accessed and updated depending on the distribution of synaptic delays. Figure (S8) shows the effect of synaptic delays on the output firing rate of a feed-forward network with different synaptic weight distributions in the presence of -function and Gaussian distributions of synaptic delays with the same mean delay. Since the network is feed-forward, the mean delay simply translates the responses in time. Higher variance in the distribution of delays relative to the membrane time constant leads to lesser overshoot of the steady state. The tail-heaviness of the synaptic weight distribution can be characterized by either constructing the moments of the distribution or a set of tail weight numbers. These numbers were defined in equations (6) and (7) and represent the fraction of neurons that spike in a neuronal population starting at rest when each neuron in the population receives excitatory inputs on average. While higher order moments of the synaptic weight distributions explain part of the variance observed in the responsiveness, a larger fraction of the variance can be explained by the tail weight numbers. Figure (S9) shows the values of the tail weight numbers as a function of (the average synaptic weight) for distributions matched for mean input current (top left) and for drift and diffusion (top right). These numbers start increasing much quicker for heavier-tailed distributions as is evident from the lower panel of Figure (S9). To test which key characteristic of synaptic weight distributions better describes the transient times, we generated 1222 random distributions between 0 and, each with a mean of 2 mV. With an input firing rate of 500 Hz, we computed the output firing rates obtained from DiPDE and the corresponding transient times taken to reach different fractions of the steady state firing rate (time taken to reach % of equilibrium firing rate denoted as here) for all these matched distributions. Figure (S10) shows a scatter plot of as a function of the first few moments and response numbers respectively for all these distributions along with the best-fit exponential curves (shown in red). Scatter plots as a function of tail weight numbers provide a better fit to the transient times than the moments. This can be seen from Table (2), which lists the sum of squared errors of the residuals for each of the best fit exponentials for different' s as a function of the first few moments and tail weight numbers respectively. Plots for the other look similar to the plots for shown in Figure (S10). Input-output curves in Figure (3) and Figure (S12) were generated using the following protocol: for a given fraction and base input rate, the DiPDE was solved first with an input firing rate. After 500 ms (at which time all runs had reached an equilibrium voltage distribution), the input rate was instantaneously changed to. The resulting output firing rate was calculated as the maximum output firing rate achieved after the instantaneous input rate change. We quantified the spread in the input-output relationship due to different values of the fraction for a given base input rate by measuring the slope, defined as follows: (28) The most comprehensive treatment of synaptic utility would involve a single 3-D integro-differential equation where the probability distribution depends additionally on the synaptic utility. However, since this formalism is computationally more expensive and given that the two separate 2-D integro-partial differential equations provide adequate results in cases where the full joint probability can be expressed as a product of the marginals over the voltage and synaptic utility, we opted to model synaptic utility as a separate differential equation via a Tsodyks-Markram mechanism [42]. This additional equation describes the probability of a synapse to have a particular utility. The synaptic utility ranges from 0 to 1, with 1 indicating maximum availability of synaptic vesicles. Its recovery is governed by (29) In the pre-synaptic population, for all probability at the membrane potential threshold, the synaptic utility is decremented using, (30) where thus represents the fraction of synaptic vesicles available after a spike. For our analysis, we use and ms respectively. Given the exponential form of equation (29), the synaptic utility axis (ranging from 0 to 1) is also discretized geometrically, similar to the voltage discretization for linear leak (equation (14) ). Thus, at each time step, the pre-synaptic probability matrix undergoes a single-index shift towards full synaptic utility (). When a non-zero proportion of the probability for a given pre-synaptic neuronal population is at threshold (i. e. there is spiking in the population), the output synaptic weight distribution from that population to the target population is convolved with the synaptic utility distribution, (31) and the utilities are subsequently modified according to equation (30). To analyze the effect of fluctuations generated by different synaptic distributions, the DiPDE numerical solution was evolved with a given synaptic distribution for 300 ms, long enough for the voltage distribution to reach equilibrium. At this point, to simulate fluctuations which occur on timescales, extra synaptic inputs are added per neuron on average to the population in a single time step. This is done by convolving the stationary probability distribution obtained from a given external input rate with the corresponding synaptic weight distribution and the instantaneous change from the equilibrium firing rate is quantified using equation (9). This same analysis was then conducted with synaptic depression during the fluctuation stage, starting with the stationary distribution obtained without synaptic depression. Unlike in the no-depression case, the synaptic distribution was scaled by the synaptic utility (using the method outlined in (Methods: Implementation of short-term synaptic depression) ) after each additional synaptic event. The utility did not recover (). The utility fraction, which represents the decrement in pre-synaptic utility after an input (equation (30) ), was calculated for each synaptic distribution separately. Denoting the peak of the relative excitability curve generated from the no-depression case by after additional inputs per neuron on average, we chose the decrement to be: (32) Intuitively, with this choice, after additional inputs per neuron on average, the relative excitability in the presence of synaptic depression gets closer to unity. All the results for DiPDE simulations were obtained using MATLAB. The entire code base is available for download at http: //download. alleninstitute. org/publications/the_influence_of_synaptic_weight_distribution_on_neuronal_population_dynamics/.
Neurons communicate via action potentials. Typically, depolarizations caused by presynaptic firing are small, such that many synaptic inputs are necessary to exceed the firing threshold. This is the assumption made by standard mathematical approaches such as the Fokker-Planck formalism. However, in some cases the synaptic weight can be large. On occasion, a single input is capable of exceeding threshold. Although this phenomenon can be studied with computational simulations, these can be impractical for large scale brain simulations or suffer from the problem of insufficient knowledge of the relevant parameters. Improving upon the standard Fokker-Planck approach, we develop a hybrid approach combining semi-analytical with computational methods into an efficient technique for analyzing the effect that rare and large synaptic weights can have on neural network activity. Our method has both neurobiological as well as methodological implications. Sparse but powerful synapses provide networks with response celerity, enhanced bandwidth and stability, even when the networks are matched for average input. We introduce a measure characterizing this response. Furthermore, our method can characterize the sub-threshold membrane potential distribution and spiking statistics of very large networks of distinct but homogeneous populations of 10s to 100s of distinct neuronal cell types throughout the brain.
The Influence of Synaptic Weight Distribution on Neuronal Population Dynamics