A Diffusion Approximation and Numerical Methods for Adaptive Neuron Models with Stochastic Inputs

Characterizing the spiking statistics of neurons receiving noisy synaptic input is a central problem in computational neuroscience. Monte Carlo approaches to this problem are computationally expensive and often fail to provide mechanistic insight. Thus, the field has seen the development of mathematical and numerical approaches, often relying on a Fokker-Planck formalism. These approaches force a compromise between biological realism, accuracy and computational efficiency. In this article we develop an extension of existing diffusion approximations to more accurately approximate the response of neurons with adaptation currents and noisy synaptic currents. The implementation refines existing numerical schemes for solving the associated Fokker-Planck equations to improve computationally efficiency and accuracy. Computer code implementing the developed algorithms is made available to the public.


INTRODUCTION
Linking neuronal membrane dynamics, synaptic input statistics and spike train statistics is a central problem in computational neuroscience. The forward problem of mapping membrane dynamics and input statistics to spiking statistics, such as firing rates, is important for understanding how biological neurons transmit sensory and motor signals. The inverse problem of mapping spike train statistics to input statistics and membrane dynamics is critical for inferring synaptic connectivity and neuron properties from extracellular recordings (Pillow et al., 2008;Pernice and Rotter, 2013). Efficient and accurate solutions to both the forward and backward problem rely on computational models of neurons that balance simplicity and computational tractability with biological realism. The adaptive exponential integrate-and-fire (AdEx) neuron model has been proposed as a balance between these needs: It is simple enough to be amenable to mathematical analysis, computationally efficient to simulate and its parameters can be fit to accurately capture the responses of diverse types of biological neurons (Fourcaud-Trocme et al., 2003;Richardson et al., 2003;Brette and Gerstner, 2005;Jolivet et al., 2008;Naud et al., 2008;Touboul and Brette, 2008;Richardson, 2009;Hertäg et al., 2014;Ocker and Doiron, 2014).
Noise is pervasive in the brain. Unreliable neurotransmitter release, ion channel noise and irregular arrival of synaptic inputs combine so that spiking activity is highly irregular with mean spike counts approximately equal to their variance (Softky and Koch, 1993;Newsome, 1994, 1998;Faisal et al., 2008). While the deterministic dynamics of the AdEx model are well-characterized Touboul and Brette, 2008), its response to the strong noise observed in cortical circuits is more difficult to analyze mathematically. Spike train statistics can be approximated from Monte-Carlo simulations, but it is computationally expensive to obtain accurate approximations. Fokker-Planck techniques offer an alternative method for approximating spiking and membrane potential statistics. The threshold integration method (Richardson, 2007(Richardson, , 2008; Richardson and Swarbrick, 2010) provides efficient numerical schemes for solving Fokker-Planck equations associated with integrate-and-fire neuron models, but is only applicable when the Fokker-Planck equation is one-dimensional. Temporally correlated inputs and adaptation currents introduce additional dimensions that prevent the direct application of threshold integration. This shortcoming is partially overcome by approximations that reduce the dimensionality of the Fokker-Planck equations: A quasi-static approximation can be used to capture the effects of slow and weak adaptation with one dimension (Richardson, 2009;Ocker and Doiron, 2014). Matched variance approximations capture the effects of temporally correlated synaptic inputs in one dimension by scaling the noise coefficient to capture passive membrane variability (Alijani and Richardson, 2011;Hertäg et al., 2014).
In this article, we present a suite of numerical algorithms for solving Fokker-Planck equations associated with the AdEx model driven by stochastic synaptic input. We extend the matched variance approach (Alijani and Richardson, 2011;Hertäg et al., 2014) to account for the membrane potential fluctuations introduced by a voltage-based adaptation current and to account for a wider class of temporally correlated input. This approximation more accurately captures the effects of strong or fast voltage-based adaptation on neuronal response statistics than previous approximations for adaptive neuron models (Richardson, 2009;Ocker and Doiron, 2014). Moreover, we developed a higher order solver for the numerical implementation of the threshold-integration scheme (Richardson, 2007(Richardson, , 2008 to solve the associated Fokker-Planck equations more efficiently. All of these numerical schemes are combined into user-friendly computer code that is publicly available at http://www.mathworks.com/matlabcentral/ fileexchange/55337-adex-threshold-integrator.

Model Description
We consider integrate-and-fire model neurons with linear adaptation, defined by Brette and Gerstner (2005) V ≥ V th ⇒ V → V re and w → w + b where Equation (3) indicates that each time V(t) reaches a threshold at V th , it is reset to V re and w(t) is incremented by a fixed amount, b. In this expression, C m is the membrane capacitance, V(t) is the membrane potential, g L is a leak conductance, E L is the leak reversal potential, w(t) is an adaptation variable with time constant τ w , b captures the spikedependence of w(t) and a captures the voltage-dependence of w(t). Our numerical methods are derived for general ψ(V), but all simulations use the adaptive exponential integrate-andfire (AdEx) model (Fourcaud-Trocme et al., 2003;Brette and Gerstner, 2005), for which Our methods are also easily adapted to the closely related Izhikevich model, which contain a quadratic, instead of exponential, non-linearity (Izhikevich, 2003). We develop numerical approximations that can be applied to any stationary stochastic input current, I(t), but all examples use a model of the form Here, J e > 0 and J e < 0 are excitatory and inhibitory synaptic weights, α x (t) is a postsynaptic current waveform (EPSC and IPSC), and t x k are spike times for x = e, i. Without loss of generality, assume that α x (t)dt = 1. We consider two instantiations of this model, one is temporally uncorrelated and the other temporally correlated.
The temporally uncorrelated input model is achieved by letting excitatory and inhibitory spike arrive as homogeneous Poisson processes with rates r e and r i , and by taking the synaptic currents to be Dirac delta functions, α x (t) = δ(t), x = e, i.
The temporally correlated input model is defined by introducing temporal correlations to the spike times and to the synaptic kinetics. Excitatory and inhibitory spike times arrive as inhomogeneous Poisson processes. The time-dependent firing rates, ν x (t), obey Ornstein-Uhlenbeck dynamics, τ ν dν x dt = −ν x + r x + σ ν 2τ ν η x (t), for x = e, i where τ ν sets the timescale of firing rate fluctuations, r x is the stationary mean rate, σ ν is the stationary standard deviation and η x (t) is standard Gaussian white noise. As long as σ ν ≪ r x , firing rates are positive with overwhelmingly large probability (Merkel and Lindner, 2010;Rosenbaum et al., 2012b). Synaptic kinetics for the temporally correlated input model are captured by setting (Dayan and Abbott, 2001)   timescales of IPSCs respectively. Synaptic time scales were chosen to mimic the kinetics of AMPA and GABA B mediated kinetics (see Table 1 and Dayan and Abbott, 2001). We consider the neuron response statistics and the accuracy of our approximations under a variety of different parameter values. The "default" parameter values, used in all simulations except where explicitly stated otherwise in figure captions and axes labels, are given in Table 1.
For all Monte Carlo simulations, Equations (1-3) were solved using a forward Euler scheme with a time bin size of dt = 0.1 ms, except where specified otherwise. Each Monte Carlo simulation was of length 65 s. The first 5 s of all simulations were not used in computing statistics, so exactly 60 s were used. Statistics were computed by averaging over 10 such Monte Carlo simulations, except where specified otherwise. All simulations and numerical computations were performed on a MacBook Pro running OS X 10.11.2 with a 2.3 GHz Intel Core i7 processor and 16 GB of 1600 MHz DDR3 RAM. Simulations were run in Matlab R2015b (Mathworks) using the mex environment to compile C code that can be run from the Matlab command line.

A Review of Spectral and Statistical Measures of Stationary Stochastic Processes
We are interested in the steady-state statistics of membrane potential and spiking activity when the input current, I(t), is modeled as a stationary stochastic process. We first briefly review the statistical measures of stationary stochastic processes used in this study. A more in depth treatment can be found elsewhere (Yaglom, 2004;Tetzlaff et al., 2008). We model the spike train of a neuron as a sum of Dirac delta functions, where t k is the k-th threshold crossing of V(t). The steady-state firing rate is given by A common measure of covariability between two processes is the cross-covariance function, defined by for stationary stochastic processes, X(t) and Y(t). Analytical computations are often simplified by transitioning to the Fourier domain, defining the cross-spectral density, and the power spectral density, C XX (f ). The analytical and numerical computation of cross-spectral densities is simplified by the relationship (Yaglom, 2004), where * denotes the complex conjugate, is the normalized, finite-time Fourier transform of X(t) with its steady-state mean, µ X , subtracted and similarly for Y T (t). This relationship is particularly useful because it allows convolutions to be treated easily. In particular, suppose that U(t) = (K * X)(t) where * denotes convolution and K(t) is a deterministic kernel with finite L 2 norm. Then for sufficiently large T where is the Fourier transform of K(t), and therefore Finally, the steady-state variance can be computed from the integral of the power-spectral density, These relationships will be helpful in computing several statistics used below.

A Review of the Quasi-Static Approximation for Stochastic Neuron Models with Adaptation Currents
We first consider the AdEx model from Equations (1-3) with temporally uncorrelated inputs. In this case, a classic diffusion approximation (Gluss, 1967;Capocelli and Ricciardi, 1971;Ricciardi and Sacerdote, 1979) replaces I(t) from Equations (4) with where µ I = J e r e + J i r i is the mean input bias, is the effective diffusion coefficient and η(t) is standard Gaussian white noise. This choice of coefficients, µ I and D I , assures that the steady-state mean of the input current is captured exactly by the input model and also that the first two infinitesimal moments of the membrane potential are the same for the diffusion approximation as they are for the temporally uncorrelated input model (Gluss, 1967;Capocelli and Ricciardi, 1971;Ricciardi and Sacerdote, 1979). Under the diffusion approximation in Equations (9), the AdEx model from Equations (1-3) represents a coupled, two-dimensional system of stochastic differential equations. Therefore, the joint probability density of V and w obeys a two-dimensional Fokker-Planck equation (see Hertäg et al., 2014 for full formulation of the two-dimensional Fokker-Planck equation). In principle, this equation could be solved numerically and the resulting bivariate probability density could be used to compute the firing rate of the neuron. However, the Fokker-Planck equation is computationally expensive and difficult to solve, in part because it has two spatial dimensions and in part because of the non-local boundary conditions associated with the threshold-reset condition in Equation (3) (Richardson, 2009;Hertäg et al., 2014;Ocker and Doiron, 2014).
Previous studies of stochastic neuron models with voltageactivated adaptation currents (Richardson, 2009;Hertäg et al., 2014;Ocker and Doiron, 2014) resolve the difficulty inherent in solving a two-dimensional equation by utilizing the fact that the adaptation time constant, τ w , is typically much larger than the membrane time constant, τ w ≫ τ m = C m /g L . This separation of timescales justifies a quasi-static approximation in which w in Equation (1) is replaced by its steady-state mean value, µ w . Under this approximation, the AdEx model can be viewed as a nonadaptive exponential integrate-and-fire model (EIF, obtained by setting a = b = 0) with an extra bias term accounting for the mean adaptation current. The only difficulty is that the steady-state mean adaptation current depends on the steady-state firing rate and steady state mean membrane potential, and vice versa. In Richardson (2009) and Ocker and Doiron (2014), this difficulty was overcome by numerically computing the fixed point of their dependence as described next.
The linearity of Equation (2) allows the steady-state mean value of w(t) to be computed exactly in terms of the firing rate and the steady state mean of V(t). Specifically, taking expectations in Equation (2) and including the spike-based perturbations from Equation (3) gives, where w is the time-dependent mean of w(t), V is the time-dependent mean membrane potential and r(t) is the instantaneous firing rate of the neuron. The last term in this equation captures the fact that w(t) is incremented by b each time the neuron spikes, as indicated in Equation (3). The steady-state mean is given by finding the unique fixed point of this differential equation, where µ V and r 0 are the steady-state mean membrane potential and firing rate respectively. Once w(t) is replaced by its steady-state mean value and I(t) is replaced by Gaussian white noise plus a bias term, µ I , Equation (1) is a one-dimensional stochastic differential equation, Thus, the steady-state probability density of V and the steady state firing rate can be computed by solving a time-independent Fokker-Planck equation, parameterized by µ w (see Appendix A for precise formulation of the one-dimensional Fokker-Planck equation). This Fokker-Planck equation can be solved numerically efficiently using the threshold integration method developed in Richardson (2007Richardson ( , 2008. We improved computational efficiency of previous implementations of the threshold integration method in several ways. First, we implemented the solver in the C programming language. To maintain usability, we used the Matlab "mex" environment so that the compiled C code can be called from the Matlab command line. Additionally, solver allows for a variable mesh size so that the mesh can be refined where the derivative of the solution is expected to be large in magnitude. Finally, we improved the numerical scheme used to implement the threshold-integration method. Previous implementations used a "modified Euler scheme" that exploits the linearity of the Fokker-Planck equation to write the solution in terms of integrals over mesh elements (Richardson, 2008). Whereas these previous implementations used the midpoint rule to approximate the resulting integrals, we used the more accurate Simpson's rule. A detailed description of the numerical solver is given in Appendix A. The performance of our implementation is compared to previous implementations below. In summary, the mean membrane potential and firing rate, µ V and r 0 , can be computed as a function of the mean adaptation current, µ w , by numerically solving the Fokker-Planck equation. Similarly, the steady-state mean of the adaptation variable can be computed in terms of the steady-state mean membrane potential and firing rate using Equation (12). Combining these two strategies, the steady state statistics can be approximated using fixed point iteration (see Appendix B). We hereafter refer to the approximation obtained this way as the "quasi-static approximation" to steady-state statistics. This approach was used previously in Richardson (2009) and Ocker and Doiron (2014).
The quasi-static approximation is expected to be accurate when adaptation is slow and weak (Richardson, 2009;Ocker and Doiron, 2014). We tested the accuracy of the quasi-static approximation with stronger and fast adaptation (see Methods) as a function of the the rate, r e , of excitatory presynaptic spikes (Figure 1). While the quasi-static approximation captured the overall shape of the dependence of steady state firing rate on r e ( Figure 1A, compare red circles and gray curve), the errors were substantial in magnitude and percentage (Figures 1B,C, gray curve).
We next show that these errors are due largely to the fact that the quasi-static approximation does not account for subthreshold, passive membrane potential variability elicited by sub-threshold, voltage-based adaptation. We then propose an improved approximation that accounts for this variability.

A Matched Variance Approximation to Sub-Threshold Adaptation Dynamics
We conjectured that the errors made by the quasi-static approximation were largely due sub-threshold fluctuations in w(t), which are not accounted for by the quasi-static approximation. Our goal in this section is to derive a diffusion approximation that accurately captures the steady-state free membrane potential variance when sub-threshold adaptation is taken into account.
To demonstrate the source of the errors made by the quasistatic approximation, first consider the free membrane potential under the quasi-static approximation, defined by removing active currents and spiking from Equation (13) to obtain This represents an Ornstein-Uhlenbeck process with steady-state variance (Gardiner, 1985) var where τ m = C m /g L . Thus, the variance of the free membrane potential under the quasi-static approximation does not depend on adaptation dynamics at all (since Equation 14 does not depend on µ w , τ w , a or b).
For the full AdEx model, the sub-threshold membrane potential is affected by fluctuations in w(t) through voltage-based adaptation. This can be seen by considering the free membrane potential and adaptation dynamics defined by omitting spikeactivation and spiking dynamics from Equations (1-3) to obtain Here, U(t) represents the free membrane potential and W(t) the "free adaptation" current. Since these equations are coupled, fluctuations in W(t) affect fluctuations in V(t) and vice versa.
In Appendix C, we show that the steady-state variance of U(t) under Equations (15) for the temporally uncorrelated input model is given by where D I = (J 2 e r e + J 2 i r i )/2 as above. This reveals a major source of error in the quasi-static approximation. As noted above, the free membrane potential variance for the FIGURE 2 | Steady-state free membrane potential variance and the error in estimating it with a quasi-static approximation. (A) Steady-state free membrane potential variance when sub-threshold adaptation is accounted for, from Equation (16), as a function of voltage-based adaptation strength, a, and adaptation timescale, τ w . (B) Percent error of quasi-static approximation to the steady-state membrane potential variance, from Equation (14), vs. the variance from (A). Adaptation strengths, a, are reported in units kHz·C m . Voltage variance has units mV 2 . All parameters other than a and τ w are as in Table 1.

quasi-static approximation (Equation 14
) does not depend on adaptation at all. On the other hand, the free membrane potential variance of the full model is affected by sub-threshold adaptation currents, since Equation (16) depends on a and τ w . Comparing Equations (14) and (16) reveals that the steady-state variance under the quasi-static approximation is least accurate when adaptation is strong (a larger relative to g L ) and/or fast (τ w smaller relative to τ m ). This difference is demonstrated by comparing the approximated and actual free membrane potential variance for different values of a and τ w (Figure 2). The quasi-static approximation does especially poorly in capturing the free membrane potential variability when a is large or τ w is small ( Figure 2B). These errors in approximating the variability of sub-threshold membrane potential variability introduce errors in the resulting firing rate approximation. Notably, the percent error made in approximating the firing rate depends similarly on the adaptation parameters as the percent error made in approximating the membrane potential variance (Figures 3C,D, compare to Figure 2B).
We now derive a one-dimensional "matched variance" diffusion approximation that accurately captures the steady-state free membrane potential variance of the two-dimensional model in Equation (15). This diffusion approximation is defined by where η(t) is standard Gaussian white noise and the diffusion coefficient, D eff , is scaled to account for the variability of V(t) in the presence of sub-threshold adaptation. The steady-state free membrane potential variance from Equation (17) is Gardiner (1985) var We want to choose D eff so that this variance is equal to the variance in Equation (16) from the two-dimensional model. This is achieved by setting Equation (16) equal to Equation (18) and solving for D eff to obtain where D I = (J 2 e r e + J 2 i r i )/2. This choice of D eff defines the "matched variance" approximation under the temporally uncorrelated input model (see Methods).
In summary, the matched variance diffusion approximation, given by using Equation (17) with D eff from Equation (19), achieves the same steady state free membrane potential variance as the two-dimensional model in Equation (15) under the temporally uncorrelated input model (see Methods). This approach is similar to previous matched variance approximations that account for the effects of temporally correlated noise, but not sub-threshold adaptation, on the free membrane potential variance (Alijani and Richardson, 2011;Hertäg et al., 2014). Below, we consider a combination of these approaches that accounts for both voltage-based adaptation and temporally correlated inputs.
Once D eff is computed, the matched variance approximation to the firing rates can be computed using the same Fokker-Planck solver and iterative methods used for the quasi-static approximation reviewed above and in Appendices A and B. The only difference between the quasi-static and matched variance approximations is the choice of diffusion coefficient, D I or D eff .
The matched variance diffusion approximation corrects a substantial portion of the error made by the quasi-static approximation, especially when the neuron is in the fluctuation dominated regime (Figure 1, compare blue and gray). This is because the fluctuations introduced by the adaptation current alter firing rates the most in the fluctuation driven regime.
Moreover, the matched variance approximation provides a substantially improved firing rate approximation whenever voltage-based adaptation is strong and/or fast (a large, τ w small; Figures 3E,F). Notably, the absolute error in firing rate is less than 1 Hz and the relative error less than 10% for the entire set of parameters we considered (Figures 3E,F). It was greater than 1 Hz and 10% for a considerable portion of the parameters using the quasi-static approximation (Figures 3C,D). Nevertheless, the matched variance approximation replaces the temporally correlated adaptation current with white noise which introduces some error, especially when adaptation is strong and fast ( Figure 3F).
The matched variance approximation offers essentially the same computational efficiency as the quasi-static approximation since they only differ by the value of the diffusion coefficient used. The average CPU time used to compute each firing rate data point in Figure 1 was 1.4 ms for the quasi-static approximation and 1.5 ms for the matched variance approximation.
A comparison of numerical schemes for solving the Fokker-Planck equation shows that the modified Euler scheme with Simpson's rule substantially outperforms the previously proposed modified Euler scheme with the midpoint rule (Figure 4, compare blue and red). Interestingly, a standard Euler scheme slightly outperforms the modified Euler scheme with the The temporally uncorrelated input model was used. All parameters other than a and τ w are as in Table 1. midpoint rule (Figure 4, compare green and red), but only gives a reasonable approximation for very fine meshes, due to the large exponential nonlinearity in ψ(V) for V near V T (Richardson, 2008).
Despite the fact that the matched variance approximation captures some of the variability introduced by sub-threshold voltage-based adaptation, errors are introduced by several assumptions made by the model. We next explore these other sources of error by modulating various parameters of the model.
Both the quasi-static and matched variance approximations approximate spike-based inputs with Gaussian white noise. This diffusion approximation is only mathematically valid in the limit of small synaptic weights, J e , J i ≪ 1, and large presynaptic firing rates, r e , r i ≫ 1 (Ricciardi and Sacerdote, 1979). To test the dependence of our numerical approximations on the synaptic weights, we computed the error in approximating postsynaptic firing rates as the synaptic weights were scaled and firing rates were scaled inversely, so that the mean input, µ I = J e r e + J i r i remained fixed. Specifically, we set J x → cJ x and r x → r x /c for x = e, i and for a range of scalings, c. We found that the matched variance approximation accurately predicted the firing rate over a large range of scalings (Figure 5 blue curve, absolute error < 5% for all data points) and substantially outperforms the quasi-static approximation (Figure 5, compare blue and gray).
While the matched variance approximation corrects for some variability introduced by sub-threshold, passive adaptation, it does not correct for fluctuations in the adaptation current evoked by action potentials. Some of this variability is introduced the reset of the membrane potential after a spike, V → V re . This reset rule affects the adaptation current, w, through the FIGURE 4 | Convergence of numerical schemes for computing steady state firing rate. (A) Relative numerical error (defined as the relative deviation of the computed firing rate from the rate computed with a fine uniform mesh of size dv = 10 −4 mV which gives 7 × 10 5 mesh points) plotted as a function of the number of mesh points. The mesh was uniform except near V re and V th where it was refined (see Appendix A). (B) Same error plotted as a function of CPU time to compute the firing rate. Firing rates were computed using the threshold integration method with a standard Euler solver (green), the modified Euler method from Richardson (2007Richardson ( , 2008 using a midpoint rule for integration (red) and the modified Euler method using Simpson's rule (blue). For the Euler method, the error is not plotted for coarser meshes because they yielded a predicted firing rate of zero. Circle indicates the mesh size used in all other figures. All parameters are as in Table 1. voltage-based adaptation. If V re is far from the steady-state mean membrane potential, the voltage-reset rule has a larger impact on the adaptation current, so we expect the matched variance approximation to perform poorly since it ignores this impact altogether. Indeed, we found that the matched variance approximation performs especially poorly for V re < −80 mV Presynaptic firing rates, r e and r i , were scaled inversely so that µ I = J e r e + J i r i remained constant. The temporally uncorrelated input model was used. All parameters other than J e , J i , r e and r i are as in Table 1. FIGURE 6 | Firing rate and error as a function of reset potential. (A) Steady state firing rate from simulations (red circles), quasi-static approximation (gray) and matched variance approximation (blue) as a function of the reset potential, V re . (B) Percent errors from (A). The temporally uncorrelated input model was used. All parameters other than V re are as in Table 1.
Surprisingly, we found that a lower reset potential (V re more negative) produced a larger firing rate, and this effect was captured by simulations and both numerical approximations ( Figure 6A). This is a consequence of strong voltage-based adaptation. Resetting the membrane potential far below E L also reduces the adaptation current, w, through Equation (2). The membrane potential recovers from its reset before w (since τ m < τ w ), during which time the neuron has an increased likelihood of spiking.
While the matched variance approximation captures some sub-threshold fluctuations introduced by voltage-based adaptation, it uses the same approximation to the effects of spike-based adaptation that the quasi-static approximation uses. Thus, we expect both approximations to perform poorly when spike-based adaptation is too strong or fast. We confirmed this expectation by comparing both approximations to Monte-Carlo simulations over a range of spike-based adaptation strengths, b, and adaptation timescales, τ w (Figure 7). As expected, the matched variance approximation performed most poorly when spike-based adaptation was both fast and strong (Figures 7E,F). The quasi-static approximation performed poorly for faster adaptation, but interestingly performed better for stronger spikebased adaptation (Figures 7C,D). Surprisingly, the quasi-static approximation actually performed better than the matched variance approximation for strong and slow spike-based adaptation.

Numerical Methods to Compute the Linear Response to a Modulation of Input Current
So far, we considered numerical methods to compute the steady state firing rate of an adaptive model neuron with stationary inputs. We now consider the response of the neuron to inputs whose statistics exhibit an explicit time-dependence. In general, this would require solving a time-dependent formulation of the Fokker-Planck equation that is not amenable to the thresholdintegration methods used for the steady state problem. However, if the temporal modulation of the statistics are weak, thresholdintegration methods can be used to solve the time-dependent Fokker-Planck equation to first order in the magnitude of the modulation (Richardson, 2007(Richardson, , 2008. This can be combined with linear response theory to quantify neuronal signal transfer, interneuronal correlations and recurrent network stability (Lindner et al., 2005;de la Rocha et al., 2007;Richardson, 2008;Shea-Brown et al., 2008;Ostojic et al., 2009;Ledoux and Brunel, 2011;Trousdale et al., 2012). We follow previously developed methods for computing the response of an adaptive neuron model to a modulation (Richardson, 2009;Ocker and Doiron, 2014), with some refinements. Details of the implementation are given in Appendix D and we review the overall approach here.
In short, we consider a weak perturbation to the synaptic input of the form where I 0 (t) is the un-perturbed input considered above and u(t) is some stationary perturbation. Linear response theory can be used to show that the finite-time Fourier transform of the firing rate can be written as for sufficiently large T where S r (f ) is the the susceptibility function of the neuron's firing rate (Risken, 1996;Fourcaud and Brunel, 2002).
The susceptibility function can be computed from a onedimensional Fokker-Planck equation, linearly coupled with an equation for the modulation of the adaptation current. This can in turn be solved using threshold-integration methods. This approach was first developed in Richardson (2009) and also used in Ocker and Doiron (2014). We review the approach and our implementation of it in Appendix D. Similar to our computation of steady-state statistics, we improved the efficiency and accuracy of the susceptibility computation by implementing the solver in C, using Simpson's method to compute integrals and using the matched variance approximation for the effective diffusion coefficient.  Table 1.
It is often necessary to take an inverse Fourier transform of susceptibility functions (see the examples below, for instance). Using Fast Fourier Transform methods for this inverse requires a uniform grid of frequencies. This grid must be fine enough to capture the low-frequencies introduced by adaptation, but must also extend to sufficiently high frequencies to capture faster spiking and membrane potential timescales. As a result, the susceptibility functions must be evaluated at a large number of frequency values. This can be computationally expensive since threshold-integration must be performed separately for each value of f (see Appendix D).
To overcome this high computational cost, we note that susceptibility functions of neuron models are typically smooth and slowly varying functions of f , especially at higher frequencies and especially when viewed on a log-log scale (see Figure 8, for example). Instead of evaluating them at each value of f , we first evaluate them at a coarser, potentially non-uniform mesh of frequencies. Since the functions change more sharply at lower frequencies, we use a mesh that is uniform on a log scale. Once the susceptibility functions are computed on this coarser mesh, interpolation is used to estimate their values on a finer, uniform mesh. The inverse Fast Fourier Transform can then be applied to the interpolated values, since the mesh is uniform and sufficiently fine. We specifically found that using piecewise-cubic Hermite interpolating polynomials (Fritsch and Carlson, 1980) via Matlab's built-in pchip command yields accurate results when the coarser mesh has only 30 points (spaced evenly on a log scale). We next verify this numerical approach by approximating spike triggered averages and inter-neuronal spike train correlations, which require the application of an inverse transform to the susceptibility functions.

Computing Spike-Triggered Average and Cross-Spectral Density between the Spike Train and an Input Current Perturbation
To test the accuracy of our approximation of the response to a modulation, we performed Monte-Carlo simulations with a weak Gaussian white noise term was added to I(t), where u(t) = 2D n η n (t) and η n (t) is Gaussian white noise. The excitatory and inhibitory conductances were left un-perturbed for this example.
We considered two ways to quantify the response to the white noise perturbation. We first considered the crossspectral density, C us (f ), between the white noise input and the postsynaptic spike train (see Methods for definition). A direct application of linear response theory shows that the crossspectral density is given to first order in D n by C us (f ) = 2D n S r (f ). Thus, the cross-spectral density can be numerically approximated by computing the susceptibility function. We also considered the spike-triggered average noise, defined by STA(τ ) = avg j u(τ − t j ) In all panels, results from Monte Carlo simulations are in red and numerical results from the matched variance approximation are in blue. The filled blue circles in indicate the coarse frequency values at which threshold-integration was applied. The smooth blue curve indicates the interpolated values used to compute the inverse Fourier transforms for the STAs. The temporally uncorrelated input model was used with the additional white noise input. All parameters except r e are as in Table 1. where t j are the spike times of the neuron. These STA can be computed from the cross-spectral density by first noting that (Dayan and Abbott, 2001) where C us (τ ) = cov(u(t), s(t + τ )) is the cross-covariance function between the white noise input and the postsynaptic spike train, which is the inverse Fourier transform of C us (f ) (see Methods and Tetzlaff et al., 2008). Putting this all together gives where F −1 denotes the inverse Fourier transform and * the complex conjugate. Thus, the spike-triggered average can be approximated directly from the susceptibility through an inverse Fourier transform of the susceptibility function. A comparison of the numerically computed values of STA(τ ) and C us (f ) to their estimates from Monte-Carlo simulations of the full AdEx model shows good agreement (Figure 8). Note that, since the cross-spectral density is proportional to the susceptibility, Figure 8 also reveals the shape of the susceptibility function and demonstrates the interpolation used (compare coarse mesh indicated by blue circles to interpolated fine mesh indicated by solid blue curve). The peak in the susceptibility function arises from resonance inherited by the spike-and voltage-based adaptation rules Richardson et al., 2003;Richardson, 2009;Ocker and Doiron, 2014). This resonance is captured by the numerical approximation because the fixed point computation of the susceptibility function takes into account the linear filtering of the adaptation currents (see Appendix D) as in Richardson (2009)  (2014). Computing the numerical approximation to the crossspectral density and spike triggered average over a larger range of presynaptic excitatory rates, r e shows that this peak is a robust feature of the adaptive model at various firing rates (Figure 9).

and Ocker and Doiron
Numerical computation of the linear response function, crossspectrum and spike triggered average was extremely efficient, with all three computed in 21 ms on average for each value of r e used in Figure 9. As noted above, this efficiency was obtained by first computing S r (f ) at a coarse and non-uniform mesh of frequencies, then using polynomial interpolation to approximate the response at a fine, uniform grid of frequencies so that the inverse Fast Fourier Transform can be used to calculate STA(τ ). Using this approach, the threshold-integration scheme was only called 30 times (one for each coarse frequency mesh point), but capturing the fast and slow timescales in STA(τ ) required interpolating these 30 values to a fine mesh of 20,001 frequency values (0 to 2 kHz with 0.1 Hz steps). To estimate the efficiency added by our interpolation scheme, we tried computing the susceptibility directly on the finer mesh (with r e = 10 kHz). This required 20,001 calls to the threshold-integration scheme and took 2.5 s. Thus, our interpolation approach is about 100 times faster than computing the susceptibility at each frequency directly.

Computing the Cross-Covariance Function between the Spiking Response of Two Model Neurons
We now consider two identical versions of the pyramidal neuron model, each identical to the individual models considered above, but with the additional assumption that a proportion c = 0.1 of the excitatory presynaptic spike times are shared between the neurons. This scenario models two neurons with overlapping presynaptic pools. These shared presynaptic spikes introduce correlations between the inputs to the two neurons, which in turn introduces correlations between their postsynaptic spike trains, where t k,j is the jth spike of neuron k = 1, 2.
Linear response theory can be used to show that the crossspectral density is given to first order in c by Lindner et al. (2005),

FIGURE 10 | Cross-covariance between two identical neurons receiving shared excitatory synaptic input. (A) Cross-covariance and (B)
cross-spectral density between the spike trains produced by two identical AdEx model neurons that share 10% of their excitatory presynaptic spike times, which occur as a homogeneous Poisson process. Computed from Monte Carlo simulations (red) and from the matched variance approximation (blue). The temporally uncorrelated input model was used. All parameters are as in Table 1. where S r (f ) is the susceptibility function of both neurons and C I 1 I 2 (f ) = c J 2 e r e is the cross-spectral density between the neurons' input currents. The cross-covariance can then be computed through an inverse Fourier transform of the cross-spectral density (see Methods).
Comparing the cross-covariance and cross-spectral density approximated numerically using this method to those computed from Monte-Carlo simulations shows reasonable agreement (Figure 10). Moreover, the numerical approximations were computed efficiently, taking less than 50 ms.
An analogous approach can be applied when the neurons share inhibitory presynaptic inputs too. Moreover, an extension of this method can be applied to compute the correlation between synaptically coupled neurons and the entire matrix of crosscorrelations in large, sparsely coupled networks (Ostojic et al., 2009;Trousdale et al., 2012). Thus, our numerical methods can be used to efficiently approximate the correlation structure of large populations of synaptically coupled adaptive model neurons.

A Matched Variance Approximation with Temporally Correlated Input Currents
So far we have considered numerical approximations obtained for a model of temporally uncorrelated synaptic input currents (see Methods). Synaptic currents in vivo are temporally correlated due to synaptic kinetics and temporal correlations in presynaptic spike timing. This shortcoming has been partially overcome by approximations that scale the diffusion coefficient to capture the steady-state free membrane potential variance of the full model (Alijani and Richardson, 2011;Hertäg et al., 2014), analogous to the matched variance approximation to adaptation dynamics described above. In these previous approaches, the approximation was only applied to synaptic input currents with an exponentially-shaped auto-correlation function, which can be used to model Poisson presynaptic spiking with exponential post-synaptic current waveforms. Moreover, these previous approaches were either applied to non-adaptive neuron models (Alijani and Richardson, 2011) or did not account for the effects of sub-threshold adaptation on the free membrane potential variance (Hertäg et al., 2014). We next extend this approach to account for sub-threshold adaptation and arbitrary stationary stochastic input.
First note that the linearity of Equations (15) allows the steadystate variance of U(t) to be computed for any stationary input I(t). Specifically, note that each line in Equation (15) implements a linear filter so that, for sufficiently large t, where c 1 (t) and c 2 (t) are asymptotically constant as t → ∞ and * denotes convolution. The kernels of the membrane potential and adaptation filters are given by where τ m = C m /g L and (t) is the Heaviside step function.
Transitioning to the Fourier domain and using Equation (6), we can write Equations (20) as and K w (f ) = 1 1 + 2πif τ w are the Fourier transforms of K v (t) and K w (t). Here, U T (f ), I T (f ) and w T (f ) are defined by Equation (5). This is a linear system of algebraic equations that is easily solved at each frequency, f , to obtain The power spectral density of U(t) is therefore given by (see Equation 7) where C II (f ) is the power spectral density of I(t). The steady-state variance of the free membrane potential is then given by This expression is valid for any stationary input, I(t), for which C II (f ) is defined and well-behaved (Yaglom, 2004).
Under the diffusion approximation in Equation (17), the free membrane potential variance is given by Equation (18). Combining these expressions shows that the diffusion approximation accurately captures the free membrane potential variance by taking which defines the matched variance approximation to temporally correlated inputs. In the absence of voltage-based adaptation (a = 0) and when I(t) has an exponential-shaped autocovariance function, this approximation is identical to the ones in Alijani and Richardson (2011) and Hertäg et al. (2014).
We consider an input model where temporal correlations are introduced by synaptic kinetics and by temporal correlations in the timing of presynaptic spikes (see Methods). Under this input model, the power spectral density of the synaptic input current is given by Rosenbaum et al. (2012b) where α e (f ) = 1 1 + 4π 2 τ 2 e f 2 is the Fourier transform of the excitatory postsynaptic current waveform,  Table 1.
is the Fourier transform of the excitatory postsynaptic current waveform and C νν (f ) = 2σ 2 ν τ ν 1 + 4π 2 τ 2 ν f 2 is the power spectral density of the firing rate fluctuations, ν(t). Under this model, the integral in Equation (21) can be approximated numerically (we used the integral command in Matlab) to obtain the corrected diffusion coefficient that defines the matched variance approximation.
Testing this matched variance approximation against Monte-Carlo simulations shows that it captures the overall shape of the dependence of firing rate on presynaptic excitatory rate ( Figure 11A, compare blue and red), but introduces large errors (Figures 11B,C blue). Nevertheless, it substantially out-performs the quasi-static approximation obtained using the diffusion coefficient in Equation (10) (Figures 11B,C gray).
The timescale of temporal input correlations is determined by the timescale of firing rate fluctuations, τ ν , as well as the timescales of synaptic kinetics, τ e , τ r,i and τ d,i . In the limit that all of these timescales approach zero, the input is approximated accurately by the temporally uncorrelated model. Thus, we expect the matched variance approximation for temporally correlated inputs to be approximately as accurate as for temporally uncorrelated inputs when these timescales are small.
To test this prediction, we compared the matched variance (and quasi-static) approximation to Monte-Carlo simulations as the input timescales were scaled (Figure 12). As expected, the matched variance approximation performs well when input timescales are fast, but not when they are slow (Figure 12, blue curves). The quasi-static approximation, obtained using the diffusion coefficient in Equation (10), does not capture any dependence of the firing rate on the input timescales and performs worse than the matched variance approximation (Figure 12, compare blue and gray curves).
These examples demonstrate that the matched variance approximation offers an improvement over approximations that do not correct for synaptic timescales at all, but there is substantial room for improvement in the approximation of firing rates of neuron models drive by temporally correlated inputs.

DISCUSSION
We derived and tested a suite of accurate and efficient numerical methods for approximating the response statistics of the adaptive integrate-and-fire neuron models with noisy synaptic input. We introduced an extension of the matched variance approximation (Alijani and Richardson, 2011;Hertäg et al., 2014) to capture the effects of sub-threshold, voltage-based adaptation currents on passive membrane variability. This approximation is more accurate than previously posed quasi-static approximations to the effects of adaptation (Richardson, 2009;Ocker and Doiron, 2014), especially when voltage-based adaptation currents are strong and/or fast.
Our approximation is also applicable in the case of arbitrary, temporally correlated synaptic inputs. This extends previous approximations, which are applicable only when inputs have exponentially-shaped auto-correlations and which do not account for the effects of voltage-based adaptation currents on sub-threshold membrane potential variability (Alijani and Richardson, 2011;Hertäg et al., 2014).
In the examples we considered, however, the matched variance approximation performed relatively poorly when strong adaptation was combined with temporally correlated inputs (Figures 11, 12).
The methods presented here rely heavily on previously developed techniques. The threshold-integration method was developed previously (Richardson, 2007(Richardson, , 2008 and we only improved the computational efficiency. Similarly, the fixed point iteration method for approximating the effect of an adaptation current was considered previously (Richardson, 2009;Ocker and Doiron, 2014) and we derived a re-scaled diffusion coefficient to account for variability in sub-threshold adaptation currents. A matched variance approximation has been developed to account for temporally correlated synaptic inputs and we extended this method to account for sub-threshold adaptation currents.
The response of model neurons to temporally correlated inputs has also been approximated using timescale separation when synaptic kinetics are much slower or faster than membrane dynamics (Moreno-Bote and Parga, 2006). However, with the FIGURE 12 | Firing rates and errors as a function of temporally correlated input timescale. (A) Steady state firing rate from simulations (red circles), quasi-static approximation (gray) and matched variance approximation (blue) as the timescales (τ e , τ r,i , τ d,i and τ ν ) are scaled by different factors. (B,C) Errors and percent errors from (A). The excitatory presynaptic rate was increased to r e = 11 kHz and the time bin width was decreased to dt = 0.05 for simulations. The temporally correlated input model was used. All parameters other than r e , τ e , τ r,i , τ d,i and τ ν are as in Table 1. exception of NMDA synaptic kinetics, synaptic and membrane timescales in cortex are not so separated, especially when one accounts for faster "effective" membrane time constant arising from the barrage of conductance-based inputs received by neurons in vivo (Destexhe and Paré, 1999;Destexhe et al., 2003;Kuhn et al., 2004).
Computational efficiency is often ignored when reporting on numerical methods for stochastic neuron models. We carefully designed our methods to be computationally efficient and tested this efficiency. As a consequence, our methods allow the computation of membrane potential response properties orders of magnitude faster than some previous methods for which timing has been reported (Hertäg et al., 2014), though our methods are potentially less accurate. While the difference between milliseconds and seconds of computation time is unimportant when computing the firing rate of a single neuron model for one set of parameters, it could be critical for applications where response properties need to be computed many times. For example, computing rates and correlations in a large recurrent network requires solving a fixed point point problem in which stationary rates and susceptibility functions are computed many times per neuron (Amit and Brunel, 1997;Brunel and Hakim, 1999;Trousdale et al., 2012). Additionally, fitting a network model to recorded multicellular data can require the repeated evaluation of rates and susceptibility functions for several parameter values and several neurons (Pernice and Rotter, 2013).
Our methods compromise some accuracy for biological realism. Specifically, the inclusion of temporally correlated noise and adaptation currents forced us to use approximations that introduced errors in the computation of response statistics. Real neurons receive temporally correlated input and have adaptation currents. Thus, the inclusion of these complicating factors is a necessity for realistic modeling. We hope that future work can refine the approximations discussed here to reduce their errors further.
One shortcoming of our approach is that we used a highly simplified synaptic input model. We only considered synaptic timescales modeled after AMPA and GABA B mediated kinetics.
NMDA mediated kinetics are much slower and are best captured by a conductance-based synapse model that depends non-linearly on voltage (Dayan and Abbott, 2001). Thus, we expect our approximations, which rely on a current-based synapse model with fast kinetics, to perform poorly for NMDA mediated synaptic kinetics. The matched variance approximation be extended to conductance-based synapse models, using approximations for the free membrane potential variance for such a models Richardson, 2004;Richardson and Gerstner, 2005;Lindner and Longtin, 2006).
Another shortcoming of our approach is that we were unable to accurately approximate the power spectral density, autocorrelation or Fano factor of the postsynaptic spike trains. The inability to compute temporal spike train statistics is primarily due to the non-renewal properties introduced by the strong and slow adaptation current as well as temporally correlated synaptic input. Methods have been developed for approximating the spiking statistics of integrate-and-fire models with adaptation currents and temporally correlated inputs (Middleton et al., 2003;Lindner, 2004;Moreno-Bote and Parga, 2006;Ly and Tranchina, 2009;Richardson and Swarbrick, 2010;Dummer et al., 2014;Shiau et al., 2015). It would be interesting to investigate whether some of these methods could be integrated with the methods developed here.
We used a linear model of adaptation in the sense that the dependence of V ′ (t) on w in Equation (1) and the dependence of w ′ (t) on V in Equation (2) is linear. This linearity was used in the matched variance approximation because it allowed the variance of the free membrane potential to be computed in closed form in the presence of sub-threshold adaptation. The quasi-static approximation has been applied to models in which the adaptation variable depends nonlinearly on V (Richardson, 2009). Extending the matched variance approach to nonlinear adaptation would be a challenge since it would require the computation of the free membrane potential variance in a nonlinear system of stochastic differential equations.
Another approach to accounting for adaptation currents would be to solve the full two-dimensional Fokker-Planck equation for the bivariate density of w(t) and V(t). Analytical approaches have been developed for spike-based adaptation when noise is weak and spiking is mean-driven . Numerical solutions to the general problem are difficult because threshold-integration cannot be applied directly to the two-dimensional equation. Instead, one could use a numerical scheme such as finite difference, finite volume or finite element methods (see a similar approach for a two-neuron model in Rosenbaum et al., 2012a). Our future work will consider these approaches.

AUTHOR CONTRIBUTIONS
The author confirms being the sole contributor of this work and approved it for publication.

A. Steady State Fokker-Planck Equation and Numerical Scheme for Its Solution
The stochastic differential equation in Equation (13) can be re-written in the form where η(t) is Gaussian white noise, is the drift coefficient and is the diffusion coefficient. The matched variance approximation is identical except the diffusion coefficient from the input current, D I , is replaced by an effective diffusion coefficient, D eff , that accounts for the fluctuating adaptation current. The Fokker-Planck equation for the stationary probability density, P 0 (v), of the diffusion process, V(t), can be written in the form of a continuity equation (Risken, 1996;Richardson, 2007Richardson, , 2008, is the stationary probability flux and r 0 = J 0 (V th ) is the stationary firing rate of the neuron. The absorbing boundary at threshold imposes the Dirichlet boundary condition, P 0 (V th ) = 0. An additional condition is given by noting that P 0 (v) is a probability density so that We now describe a numerical method for computing the steady state statistics of the AdEx model with the adaptation current, w(t), replaced by its steady-state mean, w . This is identical to solving a non-adaptive EIF model with an adjusted input current, I(t) → I(t) − w . An identical method was used in Richardson (2007Richardson ( , 2008 except that a midpoint rule was used to approximate integrals whereas we use Simpson's rule. The steady state probability flux is defined by which satisfies the trivial differential equation with boundary condition J 0 (V th ) = r 0 where r 0 is the stationary firing rate. Recall also that P 0 (V th ) = 0. Now define p 0 (v) = P 0 (v)/r 0 and j 0 (v) = J 0 (v)/r 0 . We can then re-write the Fokker-Planck equation as a system, with boundary conditions j 0 (V th ) = 1 and p 0 (V th ) = 0. The solution for j 0 (v) is given trivially as and the solution for p 0 (v) can be found numerically by integrating backwards from v = V th . Specifically, we re-write the equation for p 0 (v) as where G(v) = −A 0 (v)/D 0 (v) and H(v) = j 0 (v)/D 0 . This linear differential equation can be solved numerically for the scaled density p 0 (v). Once p 0 is found, the firing rate and the true density are given by To solve Equation (A2) numerically, we choose a lower boundary, V lb , sufficiently below E L and V re that the density is approximately zero there, then discretize the state space [V lb , V th ] into a possibly non-uniform mesh {v k } n k=1 . We denote by p k 0 the numerical approximation to p 0 (v k ). Since the last element represents threshold, v n = V th , the boundary conditions give us p n 0 = 0 and j n 0 = 1. Equation (A2) can then be integrated backwards from threshold using a numerical integration scheme.
The most straightforward integration scheme for Equation (A2) is a standard Euler scheme defined by setting where G k and H k represent the functions G(v) and H(v) evaluated at the midpoint (v k + v k+1 )/2 and where v k = v k+1 − v k . This method is only accurate for very fine meshes, due in part to the large advection term for v > V T (see Figure 4, green curve). The "modified Euler" scheme proposed in Richardson (2007Richardson ( , 2008 takes advantage of the linearity of Equation (A2) to first write an analytical solution Gdv du.
(A3) Using a midpoint rule to approximate the integrals in this expression gives an increment defined by For mesh points where G k is approximately zero, numerical inaccuracies are circumvented by replacing the fraction in this expression with its limit, as G k → 0. This scheme gives better convergence results than the Euler scheme for coarse meshes, but not for finer meshes (Figure 4).
To improve the numerical accuracy over all mesh sizes, we propose a new scheme that uses the same analytical expression (A3), but approximates the integrals using Simpson's method instead of a midpoint method. This is achieved by taking the increment defined by is the Simpson's rule approximation to the first integral in Equation (A3). Approximating the second integral is more difficult because it is a double integral. Also, H(v) = j 0 (v)/D 0 has a jump discontinuity inherited from j 0 (v). Since Simpson's method does not converge for integrals of discontinuous functions, we must evaluate H(v) at the midpoints only. Thus, the outer integral is approximated by where I k 3 is the Simpson's rule approximation to the integral of G from v k to (v k + v k+1 )/2, which is given by While this numerical scheme requires several more arithmetic operations per increment, it converges much more quickly than the other two methods, giving much more accurate results at coarser meshes using less CPU time (Figure 4). For all figures, except where indicated otherwise in figure caption, we set V lb = −100 mV and used a mesh that was uniform with step size v = 0.1 mV, except for the two mesh points surrounding V re and all mesh points above V T where the mesh size was refined to v/2 = 0.05 mV. This yielded 654 mesh points. For Figure 4, an analogous mesh was used, with v scaled by different factors.

B. Iterative Method for Computing Fixed Point of Steady-State Statistics
The steady-state mean adaptation current, µ w , can be computed from the steady-state rate, r 0 , and steady-state mean membrane potential, µ V , in closed form using Equation (12). Conversely, the steady state rate and mean membrane potential can be computed from the steady-state adaptation current by numerically solving the associated Fokker-Planck equation as described in Appendix A. To find the steady-state statistics, we used an iterative method to compute the fixed point of these mappings, similar to the approach in Richardson (2009) and Ocker and Doiron (2014).
Specifically, define to be the function determined by Equation (12) that maps mean membrane potential and firing rate to mean adaptation current. Similarly define h(w) = (v, r) to be the mapping from mean adaptation current to mean membrane potential and firing rate determined by numerically solving the Fokker-Planck equation.
Choose an initial value w 0 = 0, then iteratively calculate and for some small relative error tolerance, ǫ. We used ǫ = 10 −4 . Here, 0 < q ≤ 1 controls the size of the iterative step. Larger values of q can potentially yield faster convergence, but also make the iterations more susceptible to oscillations that do not converge. Smaller values of q converge slowly even in the best case, but are more robust. To obtain a robust and fast iterative method, we implemented a scheme that first searches for a fixed point using q = 0.25. If convergence is not realized within 50 iterations, q is successively decreased until convergence is obtained.

C. Steady-State Free Membrane Potential Variance for Temporally Uncorrelated Inputs
We now derive Equation (16) for the steady-state variance of the free membrane potential under the two-dimensional dynamics of Equation (15) and under the temporally uncorrelated input model. First note that, when I(t) is temporally uncorrelated, the first two moments of U(t) and W(t) are identical to the first two moments under the substitution I(t) → µ I + √ 2D I η(t), which gives the two-dimensional stochastic differential equation where η(t) is white noise, and C is a constant 2 × 1 vector. Thus, X(t) is a two-dimensional Ornstein-Uhlenbeck process and the steady-state covariance matrix has a known closed form (Gardiner, 1985) = det(A)BB T + QBB T Q T 2Tr(A) det (A) where Q = A − Tr(A)I 2 and I 2 is the 2 × 2 identity matrix. Substituting A and B into this expression, the steady state variance of U(t) is given by the top-right element of the covariance matrix, var(U) = (1, 1), which is easily simplified to give Equation (16).

D. Numerical Methods for Computing the Susceptibility and Linear Response Functions
We first consider the response of the neuron to a periodic perturbation, which can be combined with linear response theory to calculate the response to more general forms of modulations (Risken, 1996;Fourcaud and Brunel, 2002). Consider a complex, oscillatory perturbation of the synaptic input of the form where I 0 (t) is the un-perturbed input with stationary statistics considered above. Here, f is the frequency of the modulations and the small parameter ǫ scales the amplitude. Using a complex modulation allows the amplitude and phase shift of the response to be calculated simultaneously (Richardson, 2007(Richardson, , 2008. When ǫ = 0, we recover the steady state model from Appendices A and B. To first order in ǫ, the time-dependent part of the statistics are periodic with frequency f after an aperiodic transient (Risken, 1996;Fourcaud and Brunel, 2002). Therefore, for sufficiently large t and to first order in ǫ, we can write where P 0 (v), r 0 , µ V and µ w are the steady state statistics from above obtained by setting ǫ = 0 and computed using the methods outline in Appendices A and B.
The terms, P 1 (v), S r (f ), S V (f ) and S w (f ) are the susceptibility functions of the probability density, firing rate, mean membrane potential and mean adaptation current. We next derive approximations of these quantities by adapting previously used methods for models with voltage-gated adaptation currents (Richardson, 2009;Ocker and Doiron, 2014).
First note that, due to the linearity of Equation (2), the value of S w can be computed exactly in terms of S V and S r . This is accomplished by making the substitutions from Equation (A4) into Equation (11), collecting all time-dependent terms and solving the resulting equation to obtain S w = aS V + τ w bS r 2πif τ w + 1 .
We next note that substituting w(t) from Equation (A4) into the equation for V(t) is equivalent to a non-adaptive EIF model with input current I(t) = I 0 (t) − µ w + ǫ(1 − S w )e 2π ift + O(ǫ 2 ). With this substitution, the oscillation in of the firing rate and mean membrane potential are to first order in ǫ by where S EIF r (f ) is the susceptibility of the firing rate and S EIF V (f ) of the mean membrane potential for the non-adaptive EIF model (obtained by setting a = b = 0 in Equations 1-3).
The EIF susceptibility functions, S EIF r and S EIF V , can be computed at a given frequency, f , using the threshold integration developed in Richardson (2007Richardson ( , 2008. We again improved previous implementations of the threshold integration method by implementing the algorithm in C and using Simpson's rule to evaluate the resulting integrals, as we describe below. Note that the susceptibility functions must be evaluated with a bias term shifted by the steady state mean adaptation current, µ I → µ I − µ w , which is computed using the steady state solver described above. In summary, S V and S r can each be computed as a function of S w using Equations (A6) and S w can be computed as a function of S V and S r using Equation (A5). Similar to the steady state statistics, this yields a fixed point problem for the modulated statistics, S v , S r and S w . This fixed point problem could be solved iteratively, as done for the fixed point statistics above, but a simpler and more efficient approach is given by noting that the relationships between S r , S V and S w are linear, unlike the nonlinear dependence of r 0 and µ V on µ w for the steady state problem. Therefore, the fixed point values of S r , S V and S w are easily found in closed form by solving the linear system determined by Equations (A5) and (A6).
In summary, the susceptibility functions, S r (f ), S V (f ) and S w (f ) are computed by first solving the steady-state problem for r 0 and µ V on µ w as described above, then computing the EIF susceptibility functions, S EIF r (f ) and S EIF V (f ), using threshold integration, then solving the linear system determined by Equations (A5) and (A6).
We now describe our numerical scheme for computing the EIF susceptibility functions, S EIF r (f ) and S EIF V (f ), that are used to compute the susceptibility functions for the full adaptive model as described above.
We adapted the numerical "threshold-integration" scheme introduced in Richardson (2007Richardson ( , 2008) except that we use Simpson's rule in place of the midpoint rule to approximate integrals, as described below. We review this approach now. A more detailed description can be found in Richardson (2007Richardson ( , 2008. The modulated membrane potential, V(t), obeys a stochastic differential equation of the form The time-independent drift and diffusion coefficients, A 0 (V) and D 0 , are the same as for the steady state equation considered above and the coefficient for the time-dependent drift is A 1 = C −1 m . The Fokker-Planck equation for the time-dependent probability density, P(v, t), can be written as a continuity equation If we make the substitution for P(v, t) from Equations (A4) and expand to linear order in ǫ, we recover the steady state Fokker-Planck Equation (A1) for P 0 (v) and derive a new continuity equation for P 1 (v), given by is the linear order perturbation of the probability flux defined by J(v, t) = J 0 (v) + ǫJ 1 (v)e 2π ift + O(ǫ 2 ).
Boundary conditions are P 1 (V th ) = 0 and The modulation in the firing rate is given by S r = J 1 (V th ) and the modulation in the mean membrane potential is given by integrating vP 1 (v). Note that P 0 (v) appears in the equation for P 1 (v), so the steady state statistics must be computed first. The continuity equation can be rewritten as as system of ordinary differential equations As for the steady state equation, we reformulate the problem by defining rescaled density and flux functions. First decompose P 1 (v) and J 1 (v) as P 1 (v) = S r p r (v) + p 1 (v) and J 1 (v) = S r j r (v) + j 1 (v). The components satisfy the two systems of equations with boundary conditions p r (V th ) = p 1 (V th ) = j 1 (V th ) = 0 and j r (V th ) = 1. As for the steady state equation for p 0 (v) considered above, these equations for p 1 (v) and p r (v) can be written as −p ′ (v) = G(v)p(v) + H(v). To integrate these equations, we use the same modified Euler method with Simpson's rule that we used for the steady state density (see Appendix A). As above, we use a forward Euler step to integrate the equations for j 1 and j r . Once p 1 , j 1 , p r and j r have been found, the modulation in the firing rate can be found by noting that J 1 (V lb ) ≈ 0 assuming that V lb was chosen sufficiently below E L and E i . Thus, we have S r = − j 1 (V lb ) j r (V lb ) .
The modulation in V(t) can be found by noting that