Linear and non-linear integrate-and-fire neurons driven by synaptic shot noise with reversal potentials

The steady-state firing rate and firing-rate response of the leaky and exponential integrate-and-fire models receiving synaptic shot noise with excitatory and inhibitory reversal potentials is examined. For the particular case where the underlying synaptic conductances are exponentially distributed, it is shown that the master equation for a population of such model neurons can be reduced from an integro-differential form to a more tractable set of three differential equations. The system is nevertheless more challenging analytically than for current-based synapses: where possible analytical results are provided with an efficient numerical scheme and code provided for other quantities. The increased tractability of the framework developed supports an ongoing critical comparison between models in which synapses are treated with and without reversal potentials, such as recently in the context of networks with balanced excitatory and inhibitory conductances.

The steady-state firing rate and firing-rate response of the leaky and exponential integrate-and-fire models receiving synaptic shot noise with excitatory and inhibitory reversal potentials is examined.For the particular case where the underlying synaptic conductances are exponentially distributed, it is shown that the master equation for a population of such model neurons can be reduced from an integro-differential form to a more tractable set of three differential equations.The system is nevertheless more challenging analytically than for current-based synapses: where possible analytical results are provided with an efficient numerical scheme and code provided for other quantities.The increased tractability of the framework developed supports an ongoing critical comparison between models in which synapses are treated with and without reversal potentials, such as recently in the context of networks with balanced excitatory and inhibitory conductances.

I. INTRODUCTION
Neurons in active networks receive a barrage of competing excitatory and inhibitory synaptic input leading to a fluctuating membrane voltage and variability in the timing of outgoing spikes [1].Characterising how this incoming stochastic drive is integrated non-linearly and output spikes triggered is key to understanding singleneuron computation or how network states emerge and has been the subject of concerted theoretical effort for over half a century [2][3][4][5].Over the years, increasing biophysical details have been incorporated into an analytical framework that includes processes in the synaptic and membrane-response components of neuronal integration and spike generation.
At the synaptic level, a common approach has been to approximate synaptic amplitudes as small so that, after a Gaussian approximation is made, a Fokker-Planck approach can be used to examine the neuronal or coupled network responses [6][7][8].More recently, in part due to experimental evidence for long-tailed synaptic-amplitude distributions [9,10] and effects of presynaptic synchrony [11,12], there has been increasing interest in how finiteamplitude synaptic shot noise effects neuronal integration [13][14][15][16][17][18][19][20][21][22].Though the effect of synapses is often approximated as additive or current-based due to the reasonable desire for analytical tractability, they are more accurately implemented as conductances.Their stochastic activation therefore constitutes multiplicative noise due to the reversal-potential prefactor in the membrane current-balance equation.The aggregate conductance increase during strong presynaptic activity significantly affects the integrative properties of neurons [23], the Gaussianity of voltage fluctuations when coupled with finiteamplitude drive [24][25][26], the transmission of sensory signals [27] and, importantly, qualitatively changes the nature of the modelled balanced state between excitation and inhibition [28,29].
At the membrane level, the majority of results for the stochastic integration of synaptic drive have been derived for neurons with a linear and ohmic response, specifically for the leaky integrate-and-fire (LIF) model -see [30,31] for earlier reviews.However, from a theoretical reduction of Hodgkin-Huxley-type models, the current-voltage relationship was shown to be better captured by including an exponential non-linearity leading to the spike onset [32].It was subsequently shown experimentally that the resultant exponential integrate-and-fire (EIF) model provides a fairly accurate description of the integration properties of neocortical pyramidal cells [33] and fast-spiking interneurons [34].
Here, the framework for linear and non-linear integrate-and-fire neurons receiving exponentially distributed excitatory and inhibitory conductances is developed.The general framework previously introduced for population models [35][36][37] with linear subthreshold behaviour is first followed.However, the choice of exponential conductance distributions (extending the approach for the LIF [17] and EIF [21] models driven by exponentially distributed additive shot noise) is shown to reduce the integro-differential master equation to purely differential form.It also allows for a more direct numerical solution, avoiding the interpolation and integrative steps needed for the more general gamma-distributed amplitudes treated in earlier work on such systems [35].After introducing the model at the level of stochastic singleneuron dynamics in Section II, the reduction of the master equation to a differential form is demonstrated in Section III.In section IV and V the framework is applied to the LIF and EIF models to examine the steady-state firing rate and firing-rate response.Though analytical solutions to the master equation when both inhibition and excitation are present were not found, an efficient numerical scheme for directly obtaining the steady state and rate response is developed with Julia code [38] provided in the Supplemental Material [39].

II. DYNAMICS OF A NEURON
The neuron is isopotential with voltage v(t) and receives stochastic drive from excitatory and inhibitory presynaptic populations.The rate of charging of the membrane capacitance is proportional to the voltage derivative, the intrinsic current-voltage relation of the neuron is modelled as without history, potentially nonlinear and proportional to the function f (v), and the stochastic synaptic current s(v, t) is a function of both voltage and time as reversal potentials are included.Putting these terms together in the current-balance equation for the neuronal membrane yields a first-order nonlinear differential equation driven by multiplicative noise The spike is implemented via the integrate-and-fire mechanism: if the voltage passes a threshold v th it is directly reset to a lower value v re and a spike is registered.Two different current-voltage relationships are considered in this paper.The first is a leaky integrator where τ is the membrane time constant.The second forcing term considered is that of the EIF [32] which includes the non-linear effect of the spike onset As well as a stable resting voltage near v = 0, the model features an intrinsic spiking mechanism via an unstable fixed point.Above this value, the voltage diverges until it hits the threshold and is then reset.The quantity v T is the voltage at which the forcing term is at its minimum whereas δ T parameterises how the exponential nonlinearity grows.These two f (v) choices, corresponding to the LIF and EIF models, are plotted in Fig. 1a.

Shot noise with reversal potentials
The second term on the right-hand-side of Eq. ( 1) is now examined in detail.In this paper, a synaptic input is approximated as being unfiltered and implemented as a conductance impulse having the effect of immediately increasing (excitation) or decreasing (inhibition) the membrane potential by a voltage-dependent amplitude.In the context of active presynaptic populations, the barrage of excitatory and inhibitory input is modelled by timedependent Poissonian rates R e (t) and R i (t) such that where ϵ e and ϵ i are the reversal potentials and {t e k }, {t i k } the set of presynaptic-pulse arrival times.The quantities h e k and h i k are unitless (scaled conductance) impulses drawn from a biophysically plausible distribution.In the context of the fast-synapse limit used here, these quantities can be thought of as being equal to the timeintegral of a fast synaptic-conductance waveform divided by the membrane capacitance and therefore directly proportional to the distribution of synaptic conductance strengths Care is needed in the interpretation of the effect of s(v, t) on the voltage dynamics in Eq. (1) due to the delta functions having a voltage-dependent prefactor.Using an isolated excitatory pulse as an example both sides are first divided by (ϵ e −v) and then integrated over a short time window that includes the delta pulse.The term including f (v) will vanish with the size of this time window leaving the solution where w was the voltage before and v > w is the voltage after the pulse.This equation can be re-arranged to give the voltage jump from its initial value w as where b = 1−e −h is a convenient measure of the synaptic amplitude [35] and is bounded between 0 and 1.The amplitudes of the synaptic impulses in Eq. ( 4) are stochastic and drawn, using excitation as an example, from a distribution H e (h).Of interest, in the population-level description to be considered later, will be the fraction of synaptic impulses that bring the voltage above some value v from a lower value w.This tail distribution can be written as where the lower bound h vw is the value given in Eq. ( 6) and ensures that only jumps taking the voltage above v are included.

Exponentially distributed synaptic conductances
A biophysically plausible exponential amplitude distribution of conductance impulses is chosen here.As will be seen later, this special case of the more general gamma distribution used in reference [35] allows for a considerable simplification of both the analytical description of the master equation and its numerical solution at the population level.Using excitation as an example, the distribution is written where h e is the mean conductance amplitude.As introduced earlier in Eq. ( 7), the voltage increase can be conveniently expressed in a transformed variable b = 1−e −h where b = 0 corresponds to a negligibly small input and b = 1 an input sufficiently strong to place the voltage right at the excitatory reversal potential ϵ e .An exponential distribution for h implies a special case of the beta distribution for b where β e = 1/h e and the mean of b is b e = 1/(β e +1).The distribution parameters h e or β e can therefore be related to the mean synaptic amplitude from rest a e via Eq.( 7) which is This gives the typical jump size from a voltage w as a e (1 − w/ϵ e ).Note that with this definition the limit ϵ e → ∞ with a e held constant recovers a current-based implementation of the synaptic drive (see also Eq. 53 of the Appendix).Examples of these results, and the equivalent for inhibition, are presented in Fig. 1b

III. POPULATION DYNAMICS
The voltage trajectories of a population of neurons, each obeying Eq. ( 1) with an uncorrelated but statistically identical realisation of the stochastic drive in Eq. (4), can be described by a probability density P (v, t).The stochastic single neuronal dynamics allow for the construction of a master equation that describes the deterministic dynamics of the ensemble at the population level.As well as the probability density, it is convenient to consider the probability flux J(v, t).This describes the flow-rate of trajectories passing a particular voltage.Note that the flux at threshold J(v th ) is equal to the instantaneous spike-rate r(t) of the population with the flow then reinserted at the reset v re .These quantities are connected by a continuity equation which is statement of conservation of total density.The flux J can be resolved into a deterministic contribution equal to the forcing term multiplied by the density f P and two terms J e > 0 and J i < 0 coming from the stochastic excitatory and inhibitory synaptic events in the s(v, t) term of Eq. ( 1).This gives the flux equation The synaptic fluxes J e , J i for conductance-based synaptic drive s(v, t) can be straightforwardly derived from the amplitude distributions considered in the previous section.For example, the excitatory flux across a voltage v is the excitatory presynaptic rate times the fraction of amplitudes that bring the neuron from any lower voltage w to a voltage greater than v. Hence, where T e (v, w) is the tail distribution given in Eq. ( 8).A similar form is derivable for inhibition but will be negative.Equations (12)(13)(14) constitute the master equation for fast synaptic shot noise implemented with reversal potentials and having a general distribution for the synaptic-conductance amplitudes.Such coupled integrodifferential equations sets are generally difficult to treat analytically.In the next section it is shown that the description simplifies considerably when the underlying distribution of synaptic conductances is exponential.

Differential form of the master equation
The exponential form for H e given above (Eq.9) is now substituted into the tail equation for the distribution T e (v, w) and integrated to give Substituting in for h vw given in Eq. ( 6) yields where β e = 1/h e .It can be noted that this is a special case of the gamma-distributed conductance amplitudes considered previously [35].Following the approach in reference [17], the derivative with respect to voltage is taken and the resulting integral in one of the terms identified as being proportional to J e .This results in differential equations for the excitatory and inhibitory synaptic fluxes that take similar forms These two synaptic flux equations together with the continuity (12) and flux (13) equations describe the dynamics of an ensemble of neurons subject to exponentially distributed conductance-based shot noise.

Steady-state rate and firing-rate response
The master equation describes the full dynamics with arbitrarily strong modulations of the incoming synaptic rates R e (t) and R i (t).A full solution of the system appears difficult to obtain, given the complexity of the timevoltage partial-differential equation set and thresholdreset boundary conditions.However, simpler quantities such as the steady-state rate and firing-rate responsethe response to weak modulations of the incoming rates -nevertheless provide important information on the behaviour at the neuronal-population level and are central quantities needed for network stability and emergent oscillations [6].Using modulated excitation as an example, the incoming rate can be written in complex form as where Re is the steady-state value, Re the (potentially complex) amplitude of the modulation and ω its angular frequency.The modulatory amplitude is considered sufficiently weak such that evoked modulations in all downstream variables take (using the excitatory flux as an example) the form These expansions can be substituted into the master equation defined through Eqs.(12,13,17) with the resulting steady-state quantities providing a self-contained set of equations and the modulations providing a second set of equations.The master equations in the steady-state is with the steady-state flux equation J = f P + Je + Ji .Note that the first equation implies J = rθ(v − v re ) for v < v th .Next, at the level of a weak excitatory and inhibitory modulation of the master equation a similarlooking set of differential equations is found where now Ĵ = f P + Ĵe + Ĵi .For these quantities, the master equation has therefore been reduced to a set of coupled ordinary differential equations for which there is some hope of analytical solution or at least an efficient numerical scheme.Before treating the equations, a general statement about the limit of their behaviour at high frequency modulation is first provided.

Synaptic fluxes at high frequency
For the analysis of the high-frequency asymptotics, it is useful to consider the behaviour of the synaptic fluxes in that regime.First, it can be noted that because of the prefactor iω in the continuity equation ( 23), the modulated probability density P will decay with frequency.This allows for the following observation, again using excitation as an example: expanding the integral (Eq.16) in terms of the modulated components gives Ĵe = Re As P → 0 at high frequency, the first term becomes increasingly dominant.It can be noted that this term is simply the steady-state excitatory synaptic flux multiplied by the ratio Re / Re .This argument follows for both excitatory and inhibitory modulation and so, for increasing high modulation frequencies, the modulated synaptic fluxes can be written in terms of their steady-state values These results allow for the high-frequency firing-rate response to be related to steady-state quantities significantly simplifying the asymptotic analysis.

IV. LEAKY INTEGRATE & FIRE MODEL
For the LIF model, the current-voltage relation f (v) is linear (Eq.2, Fig. 1a) and the threshold for action potential generation v th is at the beginning of the spike with an instantaneous reset to v re .The threshold and reset are both considered to be above the resting potential, with the former criterion ensuring that the neuron is in the fluctuation-driven firing regime.
Boundary conditions.In the fluctuation-driven regime, the threshold can only be crossed by an excitatory synaptic event so the firing rate is identical to the excitatory flux at threshold: The second result -zero probability density at thresholdcomes from the flux equation ( 13) and the equality of the total and excitatory steady-state fluxes combined with zero inhibitory flux at threshold (there are no neurons with v > v th ).This requires P (v th ) = 0 given that the forcing term f (v th ) is non-zero at threshold.

LIF steady-state rate
For the LIF driven by current-based shot noise [17] it was possible to solve the master equation using a bilateral Laplace transform.This approach does not appear practical for combined excitatory and inhibitory shot noise with reversal potentials due to the additional voltage dependencies in the flux equations ( 17) compared to those for current-based drive (see Appendix Eq. 56).However, an efficient numerical scheme developed for additive Gaussian [40] and additive shot-noise [17] can be extended to account for the reversal potentials.The method works by integrating Eqs.(20)(21)(22) in the direction of convergence in two domains: The only inhomogeneous term in the equation set is proportional to r (in Eq. 20).It must be, therefore, that all quantities are proportional to r and so this factor can be scaled out.Using lower-case letters for the scaled versions, the flux is trivially given as ȷ = θ(v − v re ) from the solution of Eq. ( 20).This leaves the remaining pair of synaptic fluxes (ȷ e , ȷi ) to describe the system where the scaled probability density can be written in terms of these quantities by re-arranging the flux equation p = (ȷ − ȷe − ȷi )/f .For domain (ii) the equations (21,22) with threshold conditions (1, 0) are integrated downwards to v = 0.For domain (i) the equations are integrated up to v = 0 starting from just above the inhibitory reversal potential with initial conditions (0, −1).The excitatory flux is then matched on either side of the origin by scaling the solutions in domain (i).Finally, the unknown steady-state rate r is recovered by the normalisation condition on the probability density P (v) via r p(v)dv = 1.The steady-state rate for excitatory and inhibitory conductance-based synaptic shot noise is given for the LIF model in Fig. 2a as a function of the subthreshold mean voltage.The inset features a semilog plot on the y-axis shows that the rate is close to exponential with the mean voltage.The behaviour is also compared to that for additive shot noise (see Appendix B).In Fig. 2b the steady-state probability density and synaptic fluxes are shown for the case where r = 5Hz.Note that the distribution of the current-based case is curtailed above the reversal potential for inhibition ϵ i .

LIF firing-rate response
It was not obvious how an analytical solution to Eqs. (23)(24)(25) for the modulated rate can be obtained for combined excitatory and inhibitory conductance-based shot noise; however, it is again straightforward to generalise the numerical approach.At the modulated level, there are two inhomogeneous terms, proportional to r and either Re or Ri in the equation set (taking either excitatory or inhibitory modulation separately).The method is therefore more involved than for the steady-state case and can be found in full in the Appendix.Examples of the firing-rate response are provided in Fig. 2c for modulated excitation and in Fig. 2d for modulated inhibition.They are compared with previously derived results for current-based shot noise and shown to be broadly similar, notably in their high-frequency asymptotics.Though it was not possible to solve the equations analytically, the asymptotics for the firing-rate response for the LIF model can be derived.
Excitatory modulation at high frequencies.For the LIF, the modulated firing rate is given by the modulated excitatory flux at threshold Ĵe (v th ) as stated in Eq. ( 28).This quantity in turn becomes increasingly proportional to the steady-state excitatory flux as the modulating frequency increases (see Eq. 27).Applying this gives the result r ∼ Re Re r (30) and so the response to modulation of the presynaptic excitatory rate tends to a constant in the limit of high frequency with zero phase lag.An example of this asymptotic behaviour is provided in Fig. 2c and compared to the current-based case that has an identical form.
Inhibitory modulation at high frequencies.Though it is the incoming inhibitory rate that is modulated, the outgoing firing rate is still given by the modulated excitatory flux but with the first term in the expansion Eq. ( 26) absent as Re = 0.The modulated firing rate therefore reduces to just the second term and so r = Re where it remains only to approximate P (u) at high frequency.The argument is as follows: in Eq. ( 25) there is an inhomogeneous term that is proportional to the steady-state density P that does not decay with increasing frequency.Driven by this term, inhibitory flux is therefore dominant over the f P and Ĵe terms in the modulated flux equation, so that Ĵ ≃ Ĵi at high frequencies.
Inserting this into equation (23) gives P = −(d Ĵi /dv)/iω and then, using the result from Eq. ( 27) that Ĵi ≃ Ji Ri / Ri at high frequencies allows for substitution into Eq.(31).Following a final integration-by-parts step gives the modulated rate in terms of the steady-state inhibitory flux r ∼ Ri Re Ri It was not obvious how to further reduce this equation.Nevertheless, it shows that the high-frequency asymptotics decay with the reciprocal of frequency and with a phase of 90 • by virtue of the 1/i term and that the steady-state inhibitory flux Ji is negative.An example of this result is provided in Fig. 2d and compared to the result for current-based inhibitory shot noise which shares the same asymptotic behaviour.

V. EXPONENTIAL INTEGRATE & FIRE MODEL
The EIF current-voltage term f (v) given by Eq. ( 3) has two zeros determined by the parameters v T and δ T : a stable fixed v s just above v = 0, similar to the LIF model, and an unstable fixed point v u above v T .These fixed points can be found in terms of Lambert W functions [41] but to leading order approximation are roughly In the absence of synaptic drive, a starting voltage above v u will take the voltage to infinity in finite time which can be considered the ultimate threshold.However, physiologically, the range of validity of the exponential term is up to ∼ 10mV above the unstable fixed point (see inset to Fig. 2a in reference [33]) and so it is reasonable, as well as convenient numerically, to set a finite threshold at some value v th that is above v u but below ϵ e and at which point a spike is registered and the voltage is reset to v re .For the EIF model, therefore, the ultimate threshold v th is set at a higher value than for the LIF model.
Boundary conditions.At threshold the firing rate is given by the flux J(v th , t).However, unlike for the LIF model, in this voltage range both the (positive) deterministic component of the flux and the excitatory flux contribute to crossing the threshold where the second result is due to there being no neuronal trajectories above v th due to the instantaneous reset condition.Unlike the LIF model, there is no equality between the total flux and excitatory flux at threshold so the probability density at threshold does not vanish.Additionally, the flux equation ( 13) reduces to J = J e +J i at the two fixed points (when f (v) = 0).

EIF steady-state rate
Analytical solutions to the steady-state master equation (Eqs.[20][21][22] driven by excitatory and inhibitory shot noise were not forthcoming and so again a numerical scheme was developed.However, some care is needed due to the presence of the fixed points at v s and v u where f (v) vanishes.This requires that the steady-state master equation is integrated in three domains with matching across the domain interfaces v s and v u as well as the boundary conditions at ϵ i and v th respected.
It is convenient to first solve in domain (iii) from the unstable fixed point up to the threshold, then in domain (ii) from the unstable fixed point down to the stable fixed point and finally in domain (i) from the inhibitory reversal potential to the stable fixed point.Like for the LIF model, the only inhomogeneous term is in Eq. ( 20) and is proportional to r.All downstream quantities will therefore be proportional to this unknown quantity and it can be scaled out (lower-case letters are used for these quantites).The scaled flux is again ȷ = θ(v −v re ).It remains only to solve for the two synaptic fluxes (ȷ e , ȷi ) with the probability density written in terms of these quantities using the flux equation in the steady state, as was done for the LIF.Numerical stability requires that integration of the equation set is in the direction away from the unstable fixed point and towards the stable fixed point.This adds the complication that the ratio of excitation to inhibition is unknown at v u .As stated above, the method is to first solve the equations in region (iii) and to do so twice.For the first solution (ȷ A e , ȷA i ) initial conditions at v u of (1, 0) are chosen using pA = (1 − ȷA e − ȷA i )/f .For the second solution (ȷ B e , ȷB i ) the initial conditions at v u are (1, −1) and pB = −(ȷ B e + ȷA i )/f is used.The requirement of a vanishing inhibitory flux at threshold is then used to get the correct mix of these solutions fixing the constant κ.Once achieved, the correct combination of A and B solutions can be integrated in region (ii) from v u down to the stable fixed point v s .Matching at v s is achieved by rescaling the solution from region (i), which is found in the same way as for the LIF model.Finally, the steady-state rate r is recovered from the normalisation condition on P = r p. Examples of the steady-state rate (Fig. 3a), density and flux distributions (Fig. 3b) are provided for the conductance-based case with comparison to the case of current-based additive shot noise.

EIF firing-rate response
The master equation at the level of modulation, given by Eqs.(23)(24)(25) again appears difficult to solve analytically.However, a numerical approach similar to that used for the steady-state rate can be developed.As for the LIF case, at the level of modulation, the numerical solution is relatively involved to construct due to the multiple inhomogeneous terms and difficulty in accounting for boundary conditions and the fixed points of f (v) for the EIF model: the method is therefore described in detail in the Appendix.Example results for the numerical solutions are provided in Fig. 3c, 3d and compared to the matched current-based shot-noise model.
Modulation at high frequencies.For the case of excitatory current-based shot noise, an analytical approach was developed that allowed the rate-response at high frequencies to be related to steady-state synaptic fluxes [21].A similar framework can be developed for the case of combined excitatory and inhibitory flux and is now presented.
First, the notion of the escape time T is introduced.In the absence of the synaptic input term s(v, t) in Eq. ( 1), the time taken for the voltage to diverge to infinity from a voltage above the unstable fixed point v > v u is where Here v = v u corresponds to T = ∞ and v = ∞ corresponds to T = 0.Separately, in the same domain of v > v u , the modulated continuity equation (Eq.23) can be rewritten for P and substituted into the modulated flux equation Ĵ = f P + Ĵe + Ĵi to give These two results can now be combined by using the derivative form in Eq. ( 38) to convert the differential equation in v to one in T (v).Making use of the boundary conditions that Ĵ → r in the limit v → ∞, T → 0 and At high frequencies, using Eq. ( 27), the asymptotics can be written in terms of the steady state fluxes allowing for the high-frequency asymptotics -a parameterisation of the dynamics -to be derived if an analytically convenient form for the steady-state synaptic fluxes can be found in terms of the escape time T (v).This is in the nature of a fluctuation-dissipation relation (though see reference [42] for a complete relation for neuronal systems).Unfortunately, this does not appear to be as straightforward as for the current-based case [21] in which the asymptotics had a dependency on the relation between the ratio of the mean excitatory amplitude and the spike sharpness δ T .However, the broad behaviour remains as can be seen in Fig. 3c, 3d in which the rateresponse of the models with conductance or and currentbased drive have similar forms over moderate frequencies.

VI. DISCUSSION
The population-level framework for linear and nonlinear neurons driven by exponentially distributed excitatory and inhibitory synaptic condutances was examined.For this biophysically reasonable choice of conductance distribution it was demonstrated that the generic integro-differential master equation (see reference [35]) can be reduced to a set of more tractable differential equations.This description of the master equation allowed for the development of an efficient scheme for the numerical derivation of the steady-state density distributions and rate as well as those at the level of the firingrate response to weakly modulated presynaptic rates.The analyses of these equations presented here extend treatments of the effects of synaptic conductance on the LIF and EIF firing rates previously derived in the smallamplitude appproximation [3,40,43].They also extend previous results for finite-amplitude synaptic drive for the LIF [17] and EIF [21] models where synapses were treated as current-based and the noise was therefore additive rather than conductance-based and multiplicative.
For the case of additive shot-noise, it was possible to derive analytical forms for the steady-state and modulated rates [17] for the LIF with excitatory and inhibitory input.These results were later extended to non-exponentially distributed inhibition [18,19].Unfortunately, that analytical approach, which used bilateral Laplace transforms, does not appear to extend straightforwardly to the case of reversal potentials: analytical solutions for these quantities remain a topic for future research.Interestingly, the steady-state firing rate as shown in the inset to Fig. 2a appears close to exponential, suggesting that there may nevertheless be a simple approximation, if not a full solution, to the steady-state rate with combined excitation and inhibition An analytical solution for the rate-response for the LIF model was similarly illusive, though the high-frequency asymptotics could be obtained.The difficulty in finding an analytical solution of course also extends to the more detailed EIF model, though again an intriguing nearexponential rise in the steady-state firing rate can be seen in the inset to Fig. 3a This suggests that a simple solution (or, at the least, an accurate approximation) to the steady-state rate might be possible.
It was previously shown [21] that, at the level of the firing-rate response for the EIF model with current-based drive, there is an interplay between the mean amplitude a e and the sharpness of the spike δ T in setting the asymptotic high-frequency exponent.Though it was again possible to show that the asymptotics of the high-frequency response can be related to the steady-state excitatory flux -a form of fluctuation-dissipation relation -the simplification of the excitatory flux for the case of reversal potentials was not so straightforward to derive as for the current-based case because the mean synaptic amplitude is conditional on voltage.However, the dependency on the rapidity of the response as a ratio of the amplitude around threshold and the spike sharpness broadly holds over moderate frequencies (Fig. 3c).
These results represent an initial foray into models of linear and non-linear integrate-and-fire models with exponentially distributed conductances.While it is clear from the structure of the master equation that the steady-state voltage distribution can be trivially derived in integral form for both the LIF and EIF in the case of excitation only, in the context of synaptic reversal potentials it is the inclusion of inhibition that is of relevance.The solution for the steady-state rate for combined ex-citation and inhibition as well as the firing-rate response remain therefore an important theoretical goal.
A further case of interest would be to relate a changing amplitude distribution to the role of short-term synaptic plasticity [44,45].This has been analysed in the Gaussian approximation [46,47] but would be interesting to examine as a driver for how the amplitude distribution changes as a function of presynaptic activity, the history dependence of inputs at particular fibres and the effect of non-Poissonain statistics on the population equations [48,49].The high-frequency asymptotics of the rate response are sensitive to the interplay of the spikesharpness and mean synaptic amplitude [21].Including a mechanism like short-term plasticity provides a motivation for analysing how the rapidity of the neuronal response can be modulated by the level of presynaptic network activity and is a physiologically valid question for future research.

APPENDIX A: NUMERICS FOR THE FIRING-RATE RESPONSE
The numerical solution for the steady-state master equation (20)(21)(22) for the conductance-based shot-noise LIF and EIF models were described in the main text.The approach is an extension of the Threshold Integration methods developed for Gaussian-white noise [40], LIF with additive shot noise [17] and EIF with additive excitatory shot noise [21].In this section, the approach for modulated excitatory or inhibitory conductance-based shot noise is outlined.The general components common to the LIF and EIF are first described and then the boundary conditions specific to the two models explained in detail.All Julia code [38] is provided in the Supplemental Material [39].
The method now presented is similar to that used for the steady state, though in the modulation case there are three inhomogeneous components proportional to r, Re and Ri in the region (ii) for the LIF and (ii) and (iii) for the EIF (see Eqs. 29 and 36).The system is linear so it suffices to solve the problem with either excitatory modulation or inhibitory modulation separately, with effects of combinations of the remaining two modulations simply adding at the population level in the weak modulation approximation.Throughout the following, the example of excitatory modulation will be used without loss of generality.There are therefore two inhomogeneous terms that are considered together: that of r, Re allowing the solutions for the various fluxes to be separated into a sum of two subsolutions, like for the total flux: Equations (23)(24)(25) can be resolved into two sets of equations with variables (ȷ r , ȷr e , ȷr i ) and (ȷ e , ȷe e , ȷe i ).The first set has r = 1 and Re = Ri = 0 and the second set has r = 0 and Re = 1 with Ri = 0 still because here only excitatory modulation is considered.The solution approach is relatively straightforward for the LIF in region (ii) because the boundary conditions can be imposed at the threshold, which is the starting point for the integration back to the stable fixed point.For the EIF, however, there is an added complication because it is necessary to integrate from the unstable fixed point (where boundary conditions are not directly specified) up to the threshold where they are -the strategy for handling this complication is explained in detail later.
In domain (i), from the reversal potential to the stable fixed point for the LIF or EIF, it is numerically convenient to run the integration from a lower-bound v lb that is slightly above the inhibitory reversal potential ϵ i .A zero-flux condition is imposed at v lb .The excitatory flux is zero here (as there are no neuronal trajectories below) so the inhibitory flux balances the deterministic component J i (v lb ) = −f P .Calling the magnitude of this inhibitory flux q, the equations in region (i) for both the LIF and EIF can be resolved into two components.So, again using the total flux as an example: Ĵ = qȷ q + Re ȷe .
Once each subsolution has been obtained with the appropriate boundary conditions at v lb (see later) the equation in domain (i) can be scaled to match with that in (ii) at the stable fixed point (v = 0 for the LIF, v = v s for the EIF) where the ± means either just below or above the stable fixed point.These two equations can be trivially solved to finally provide q and the desired modulatory rate r in response to presynaptic excitatory modulation.The case of inhibitory modulation is derived using an identical approach but with Re = 0 and Ri = 1.The specific boundary conditions for the LIF or EIF model are now described in more detail.

LIF boundary conditions
For the LIF there are two domains to integrate over (see Eq. 29).Starting with domain (ii), first integrate 0 ← v th with initial conditions (1, 1, 0) for the (ȷ r , ȷr e , ȷr i ) solution and (0, 0, 0) for the (ȷ e , ȷe e , ȷe i ) solution.Then integrate in domain (i) up from v lb → 0 with initial conditions (0, 0, −1) for the (ȷ q , ȷq e , ȷq i ) solution and (0, 0, 0) for the (ȷ e , ȷe e , ȷe i ) solution.The matching criterion then gives r as required from the solution of Eqs.(44,45).The approach for inhibitory modulation is analagous.

EIF boundary conditions
The method is similar to that described above for the LIF except that, like for the steady-state case, the so-lution in the domain (iii) from v u → v th needs to be properly constructed before the solutions in the other domains (ii) and (i) are found.Again, for simplicity of exposition, a case is considered where there is excitatory modulation only.
Region (iii).Stability requires integration to be in the direction v u → v th .However, the boundary conditions, that Ĵ = r and Ĵi = 0, are imposed at v th .Hence, linearly independent solutions with different initial conditions at v u need to be appropriately combined so that the desired conditions at v th are met.First, a solution (ȷ E , ȷE e , ȷE i ) to Eqs. (23)(24)(25) with Re = 1, Ri = 0 are integrated from v u → v th with initial conditions (0, 0, 0).Two other solutions, (ȷ A , ȷA e , ȷA i ) and (ȷ B , ȷB e , ȷB i ) are then integrated with Re , Ri = 0, 0, with initial conditions (1, 1, 0) and (0, 1, −1) respectively.The solution that deals with the inhomogeneous term proportional to Re can now be constructed as a linear combination of these three solutions such that, for example, the total flux is written ȷe = aȷ A + bȷ B + ȷE . ( Given freedom to vary a and b (the two homogeneous solutions), a scenario is chosen where both the total and inhibitory flux vanish at threshold for this component This fixes the values of a and b in terms of the integrated solutions at v th .The r component (ȷ r , ȷr e , ȷr i ) can also be constructed from the A and B subsolutions, so for example for the total flux ȷr = cȷ A + dȷ B . ( Given the way the Re component was constructed to have vanishing total flux at threshold, the r component now has to satisfy the boundary conditions ȷ = 1 as well as ȷi = 0. Hence, the requirements that together fix c and d.These provide the correct combination of solutions for the two inhomogeneous solutions proportional to Re and r in domain (iii).

Region (ii).
Here the integration is downwards (v s ← v u ) with the initial conditions from the previous case being It should be remembered that the initial conditions for the solution (ȷ E , ȷE e , ȷE i ) were (0, 0, 0) and so don't contribute to the first of the conditions above.Both sets of equations, for the r and Re inhomogeneous solutions, are integrated down to v s .Note that because v re is in this domain, the integration for the r solution includes the Dirac-delta function in Eq. 23.Region (i).As explained above, a lower bound v lb is imposed just above the inhibitory flux.At this point, all fluxes will be zero except for the inhibitory one which is set as being proportional to a quantity q (see Eq. 43).The initial conditions for the solution (ȷ q , ȷq e ) are (0, 0, −1) and (0, 0, 0) for (ȷ e , ȷe e , ȷe i ).Finally, r and q are determined from the linear equations (44,45) given earlier thereby completing the numerical solution.

APPENDIX B: ADDITIVE SHOT NOISE
The framework for additive, current-based shot noise with exponentially distributed amplitudes was developed in reference [17] for the LIF driven by excitatory and inhibitory input and later for the EIF driven by excitatory input only [21].In this section, the analysis in [21] is extended to combined excitation and inhibition.For current-based drive, the amplitude distribution is independent of voltage with the drive term in Eq. ( 1) written as where a e k is the amplitude and t e k the time of the kth excitatory input, with a similar definition for inhibition.The voltage amplitudes are drawn from exponential distributions A e (a) with mean a e > 0, using excitation as an example, so that The tail distribution T e (a) is the probability that an amplitude is greater than a and, unlike for the conductancebased case, is independent of voltage.The inhibitory jumps are also exponential distributed, though with a i < 0. Example distributions are provided in Fig. 1b lower panel for comparison to the case with reversal potentials.

Synaptic flux equations for additive shot noise
The excitatory flux across a voltage v is equal to the rate that jumps from all values of the voltage w < v cross it.For current-based shot noise, this is just a convolution of the probability density and tail distribution, which for exponentially distributed input can be written with a similar form derivable for inhibition.This integral form can be recognized as the solution of a first-order linear differential equation for J e with the density P acting as an inhomogeneous term, with a similar result for inhibition ∂J e ∂v + J e a e = R e P and dJ i dv Together with the continuity and flux equation (Eqs.12-13) these differential equations constitute the master equation that fully describes the dynamics of an ensemble of neurons subject to current-based shot noise with exponentially distributed amplitudes.

LIF with additive shot noise
The master equation for this model can be solved to find both the steady-state rate and firing-rate response through a bilateral Laplace transform of the voltage [17].Both rates can be written in terms of the voltage moment-generating function Z(s) This is the shot-noise generalization of the simplified form [6] of the Ricciardi formula [5].The corresponding firingrate response for weak modulation of either the excitatory or inhibitory presynaptic rate is written rκ = Rκ τ r where κ = e, i for excitation or inhibition.Though these analytical forms exist for the LIF with additive shot noise, it is often more convenient to generate the solutions numerically using a similar method to that described in Appendix A for conductance-based shot noise.The principal difference is that the lower bound is not constrained to be above any reversal potential for inhibition but rather chosen to be sufficiently low that it has little effect on the results.Finally, in the limit of high frequencies, the asymptotics can be shown to be re ≃ r Re R e and ri ≃ r Ri iω a i a e − a i (61) for excitatory and inhibitory modulation.

EIF with additive shot noise
Unfortunately, the Laplace-transform solution used for the LIF does not transfer easily to the EIF with additive shot noise.The solutions are therefore obtained numerically, in the same way as for the conductance-based case but, again, with the lower bound v lb no-longer constrained by any inhibitory reversal potential ϵ i .Though a numerical approach is required for the full solution, the firing-rate response asymptotics can nevertheless be obtained, as was shown for the case of the EIF driven by only excitatory current-based shot noise [21].In the following section these results are updated for the case of excitatory and inhibitory input.
Excitatory modulation at high frequencies.For the excitatory-only drive it was shown [21] that the asymptotics for excitatory modulation depend non-trivially on the ratio of the excitatory synaptic amplitude a e to the spike sharpness δ T .The method used in reference [21] is identical when background steady-state inhibition is included with the results taking the same form where it should be remembered that a i < 0. Well above the unstable fixed point the voltage varies as P (w) ≃ (rτ /δ T )e −(w−vT)/δT .Substituting in this form and performing the integral gives an approximation valid at large voltages: Ji ≃ Ri a i δ T − a i rτ e −(v−vT)/δT ≃ Ri a i δ T − a i rT (v) (67) In the second form, the escape time T (v) has been substututed in by noting that f ∼ δ T e (v−vT)/δT /τ so that T ∼ τ e −(v−vT)/δT at large voltages.Performing the integral transform (see Eq. 41) gives the high-frequency asymptotic r ≃ Ri rτ iωτ for modulation of the inhibitory rate modulation.These results are presented in Fig. 3.

Numerics for the EIF with additive shot noise
In cases where the analytical solution cannot be found in convenient closed form or even for reasons of numerical convenience, the solutions can be found for the steady-state rate or firing-rate modulation using the same method as for the conductance-based case described in Appendix B. The two differences are: first, the forms of the flux equations so that Eqs.(56) are used instead of Eqs.(17); and second, the position of the lower bound v lb .Because there is no inhibitory reversal potential, the lower bound v lb can be placed at some sufficiently low value such that it does not have a material effect on the results (i.e. in a sufficiently hyperpolarised region where the probability density is very low).Julia code [38] for the EIF with additive shot noise is also provided in the Supplemental Material [39].

APPENDIX C: SIMULATIONS
Simulations of the voltage dynamics (Eq. 1) with synaptic-amplitude distributions drawn from Eqs. (9) or (54) were performed using a forward Euler scheme to provide a comparison to analytical or numerically exact results.All simulations were run at a time step of 0.01ms and averaged such that results were of the order of the symbol size in figure panels.For the oscillatory responses the amplitudes of the modulated rates for the LIF model were Re = 0.025kHz and Ri = 0.075kHz for both conductance-based and current-based modulations.For the EIF model these parameters were 0.075kHz in all cases.

FIG. 1 :
FIG. 1: Comparison of the Leaky and Exponential Integrateand-Fire (LIF and EIF) models driven by exponentially distributed conductance or current-based synaptic shot noise.(a) Current-voltage relationship f (v) normalised by the capacitance for the LIF and EIF models.(b) Excitatory (green) and inhibitory (red) synaptic-amplitude distributions as a function of three initial voltages (−5, 0, 5mV) for exponentially distributed synaptic conductances (upper panel, Eq. 10) or currents (lower panel, Eq. 54) with distributions having mean amplitudes matched at v = 0. (c) Voltage dynamics for the LIF (upper panel) and EIF (lower panel).Implementations with conductance-based (Eq.4, black lines) and current-based (Eq.53, grey lines) drive for the distributions in panel (b) exhibit clear differences in membrane voltage and firing times.Parameters: τ = 20ms, δT = 1mV, vT = 10mV, vre = 5mV, v th = 10mV for the LIF or v th = 20mV for the EIF; reversal potentials ϵe = 60mV, ϵi = −10mV; mean synaptic amplitudes as marked; steady-state presynaptic rates for panel (c) were Re = 0.2kHz and Ri = 0.02kHz.Code used to generate the figure is provided in the Supplemental Material [39].

FIG. 2 :
FIG.2: Leaky I&F neuron driven by conductance or current-based shot noise with excitatory and inhibitory reversal potentials.(a) Steady-state outgoing firing rate with the incoming presynaptic rates ( Re, Ri) the same for the conductance and currentbased cases.The presynaptic rates were parameterised, for convenience, from the current-based equations (see Eq. 58) for the voltage mean ⟨v⟩ with the standard deviation fixed at σ = 5mV.Note the near-exponential rise for the conductance case seen in the inset.(b) Probability density together with excitatory and inhibitory fluxes for neurons firing at 5Hz (intersections with dotted line in panel a).The presynaptic rates ( Re, Ri) were (0.393, 0.650)kHz and (0.365, 0.762)kHz for conductance and current-based cases, respectively.(c) Amplitude and phase of the firing-rate response to modulation of the presynaptic excitatory rate with asymptotics shown.(d) Same, but for the case of modulation of the inhibitory presynaptic rate.For panels (c-d) the presynaptic rates from panel (b) were used (corresponding to r = 5Hz) with all other parameters for the LIF provided in caption to Fig.1.Simulational results are provided by symbols in panels (a,c,d) and histograms in upper panel (b).All code used to generate the figure is provided in the Supplemental Material[39].

FIG. 3 :
FIG.3: Exponential I&F neuron driven by conductance or current-based shot noise with excitatory and inhibitory reversal potentials.(a) Steady-state outgoing firing rate with the incoming presynaptic rates ( Re, Ri) the same for the conductance and current-based cases.The presynaptic rates were parameterised, for convenience, from the current-based equations (see Eq. 58) for the (subthreshold) voltage mean ⟨v⟩ with the standard deviation fixed at σ = 5mV.As for the LIF, again note the near-exponential rise of r for the conductance case seen in the inset.(b) Probability density together with excitatory and inhibitory fluxes for neurons firing at 5Hz (intersections with dotted line in panel a).The presynaptic rates ( Re, Ri) were (0.446, 0.440)kHz and (0.397, 0.636)kHz for conductance and current-based cases, respectively.Note the relatively curtailed distribution near the inhibitory reversal potential ϵi = −10mV for the conductance-based case.(c) Amplitude and phase of the firing-rate response to modulation of the presynaptic excitatory rate with asymptotics shown.(d) Same, but for the case of modulation of the inhibitory presynaptic rate.For panels (c-d) the presynaptic rates from panel (b) were used (corresponding to r = 5Hz) with all other parameters for the EIF provided in caption to Fig. 1.Simulational results are provided by symbols in panels (a,c,d) and histograms in upper panel (b).Code used to generate the figure is provided in the Supplemental Material [39].

A
e (a) = θ(a) e −a/ae a e and T e (a) = θ(a)e −a/ae .(54) 57) giving the steady-state voltage mean and variance as ⟨v⟩ = a e τ Re + a i τ Ri and Var(v) = a 2 e τ Re + a 2 i τ Ri .(58) Fixing these two quantities specifies Re and Ri which is used as the basis of the steady-state rate in Figs 2a, 2b.The steady-state firing rate can be written in terms of an integral over e s −e svre .
−e svre e iωt (60) r ≃ Re rτ iωτ a e δ T − a e for a e < δ T (62) r ≃ Re rτ iωτ log(iωτ ) for a e = δ T (63) r ≃ Re rτ I e (iωτ ) δT/ae Γ δ T a + 1 for a e > δ T .(64)Theintegral I e is a function of the steady-state densityI e = ∞ −∞ due (u−v T )/ae P (u)/rτ (65)and therefore includes the presence of both excitatory and inhibitory drive.Inhibitory modulation at high frequencies.As inhibitory modulation was not considered in reference[21] a little more detail is provided.The steady-state inhibitory flux can be written Ji = − Ri ∞ v dwP (w)e −(v−w)/ai (66)