The Morris–Lecar neuron model embeds a leaky integrate-and-fire model

We show that the stochastic Morris–Lecar neuron, in a neighborhood of its stable point, can be approximated by a two-dimensional Ornstein–Uhlenbeck (OU) modulation of a constant circular motion. The associated radial OU process is an example of a leaky integrate-and-fire (LIF) model prior to firing. A new model constructed from a radial OU process together with a simple firing mechanism based on detailed Morris–Lecar firing statistics reproduces the Morris–Lecar Interspike Interval (ISI) distribution, and has the computational advantages of a LIF. The result justifies the large amount of attention paid to the LIF models.


Introduction
Much effort has been made to create a realistic but still easily computed stochastic neuron model, primarily by combining subthreshold dynamics with firing rules.The result has been a variety of, usually one dimensional, leaky integrate-and-fire (LIF) descriptions with a fixed membrane potential firing threshold [4,11,18,19], or with a rate of firing depending more sensitively on membrane potential [15,21].These models are useful both for obtaining analytical results and for ease of simulation.
By contrast, the two-dimensional stochastic Morris-Lecar (ML) neuron model, a simple cousin to the more detailed Hodgkin-Huxley (HH) model, describes the dynamics of firing in a way more closely motivated by the biology.It has been better respected by biologists than the LIF class of models, but has received little attention owing to the difficulty of mathematical analysis of this rather complicated stochastic dynamical system.
In Section 4 of this paper we show that in fact a LIF model is embedded in the ML model as an integral part of it, closely approximating the subthreshold fluctuations of the ML dynamics.This result suggests that perhaps the firing pattern of a stochastic ML can be recreated using the embedded LIF together with a ML stochastic firing mechanism.We construct such a model in Section 5 and 6, and show in Section 7 that its Interspike Interval (ISI) distribution is similar to that of the ML.Our model, while of the type described in our first paragraph, combines the realism of the ML with the ease of analysis and computation of a one dimensional LIF-type model.The work invested in LIF models is further justified by this new model.
Before we set up our stochastic ML model and write analytical details, let us have an informal look at how it works.The principal dynamics of the ML, in the central range of the input current, consist of a stable limit cycle (Fig. 1A) corresponding to firing, which encloses a stable fixed point.In between there loops an unstable limit cycle.The path of the stochastic model has two quasi-stable patterns (Fig. 1B).One is succesive firings, where the dynamics makes "large" noisy circuits around the stable limit cycle, the other is membrane fluctuations between spikes, where the dynamics makes "small" noisy circuits around the fixed point inside the unstable limit cycle.The system would continue forever in one of these two patterns were it not for the noise which causes switching from firing to subthreshold fluctuations and back again at random times when the dynamics cross the unstable limit cycle.Our analysis will show that the dynamics between spikes, of random cycling inside the unstable limit cycle followed by crossing to the stable limit cycle outside it, can be identified with the sample path behavior of a two-dimensional Ornstein-Uhlenbeck (OU) process times a rotation.
A main ingredient in our result is the stochastic dynamical phenomenon that oscillations which damp to a fixed point in a deterministic system will be sustained by the stochasticity in a corresponding stochastic system.Damped oscillations in a two-dimensional system are signalled by a local linear structure defined by a matrix having a pair of conjugate complex eigenvalues with negative real part.A corresponding stochastic system will not damp, being prevented by the noise.Instead, a quasi-stationary stochastic process is set up, which cycles in a random pattern around the fixed point.Using recent results of [2] we are able to identify, approximately, this stochastic process which is part of the subthreshold dynamics of the ML.Up to a fixed linear transformation, the approximating process is the product of a steady fast rotation with a two-dimensional OU process.The identification allows us to cement in place the correspondance, for a particular set of model parameters, a particular LIF model as the appropriate subthreshold phase between ML firings.

The Morris-Lecar model
There exists a large variety of modeling approaches to the generation of spike trains in neurons (see e.g.[6,11,14]).Most famous is the Hodgkin-Huxley (HH) model [13] consisting of four coupled differential equations, one for the membrane voltage, and three equations describing the gating variables that model the voltage-dependent sodium and potassium channels.A large amount of research effort is currently directed towards understanding how neural coding carries information through nervous systems.Basic to the subject is how single neurons transmit information.As in any modeling effort, we must ignore or summarize details and focus on what, we hope, are a few essential aspects.The ML model [20] has often been used as a good, qualitatively quite accurate, two-dimensional model of neuronal spiking.It is a conductance-based model like the HH model, introduced to explain the dynamics of the barnacle muscle fiber.The original ML model was three-dimensional, including a fast responding voltage-sensitive Ca 2+ conductance, and a delayed voltage-dependent K + conductance for recovery.To justify the two-dimensional version, one uses that the Ca 2+ activation moves on a much faster time scale than the other variables, and can conveniently be treated as an instantaneous variable, by replacing it by its steady-state value given the other variables.
The parameter values in our computations were chosen from [22,23], and are given in Table 1 together with the interpretation of variables and parameters.The variable V t represents the membrane potential of the neuron at time t, and W t represents the normalized conductance of the K + current.This is a variable between 0 and 1, and could be interpreted as the probability with the auxiliary functions given by Equation ( 1) describing the dynamics of V t contains four terms, corresponding to Ca 2+ current, K + current, a general leak current, and the input current I.The functions α(•) and β(•) model the rates of opening and closing, respectively, of the K + ion channels.The function m ∞ (•) represents the equilibrium value of the normalized Ca 2+ conductance for a given value of the membrane potential.In Fig. 1A the phase-state of the model is plotted.The system has two stable attractors; a stable fixed point corresponding to quiescence of the neuron, and a stable limit cycle corresponding to repetitive firing.In between the two attractors is an unstable limit cycle, which splits the state space into two parts from either of which the deterministic process cannot escape, once trapped there.membrane voltage V t normalized conductance W t q FIGURE 1: Phase-state plots of the normalized conductance Wt against membrane voltage Vt.The full drawn magenta curve is a stable limit cycle, the dashed magenta curve is an unstable limit cycle, and the magenta point is a stable fixed point.Black curves are sample trajectories.Panel A: model without noise, (1)- (5).If the process is started between the stable and the unstable limit cycle, or outside the stable limit cycle, the solution is seen to spiral out, respectively in, towards the stable limit cycle, corresponding to repetitive firing of the neuron.If the process is started inside the unstable limit cycle, the solution spirals into the stable fixed point, corresponding to subthreshold fluctuations of the neuron.Note that three trajectories are plotted.Panel B: model with noise, (1), (3)-( 5) and (8), σ * = 0.05.Only one trajectory is plotted, and the solution is seen to switch between periods of firing and quiescence.model, by potential dependent averages.However, it has become apparent that the stochastic nature of ion channels must be explicitly modeled if we are to capture essential features of neuron dynamics.Changes in the states of channels cannot be tracked explicitly because of their vast number.Hence, it is useful to model the role of ion channels as a stochastic process, W t , the proportion of channels open at time t.We therefore add channel noise by changing the ordinary differential equation system (1) -( 5), to a stochastic differential equation system, replacing the conductance equation ( 2) by where B t is a standard Wiener process, and the function h(•) has to be chosen.The diffusion coefficient h(•) in ( 6) should be based on the drift coefficient which gives the rate of change of fraction of open ion channels due to random openings and closings.A natural choice of the function h(•), following the diffusion approximation of [16], would be the square root of the sum of the two rates in the drift coefficient, times a factor 1/

√
N where N is the number of ion channels involved.However, this choice has the problem that it is not zero when all the channels are closed, and the resulting (6) would produce negative solutions with positive probability.To avoid this difficulty, for fixed V t we let W t be a Jacobi diffusion.In fact, in the class of Pearson diffusions [9], i.e. one-dimensional diffusions with linear drift, and with h 2 (•) a polynomial of at most degree two, this is the only bounded diffusion.Living on (0, 1), it has the form where θ > 0 and µ ∈ (0, 1).It is named for the eigenfunctions of the generator, which are the Jacobi polynomials.It is ergodic provided that γ 2 ≤ min(µ, (1 − µ)), and its stationary distribution is the Beta distribution with shape parameters µ/γ 2 and (1 − µ)/γ 2 .It has mean µ and variance γ 2 µ(1 − µ)/(1 + γ 2 ).In our case, because the diffusion coefficient in ( 7) should be of the same order as the one given by the Kurtz approximation [16], γ is proportional to 1/ √ N .By equating the drift terms in ( 6) and ( 7), we have (6) will stay bounded in (0, 1).Since α(V t ) and β(V t ) are strictly positive, we can put 2 , with σ * ∈ (0, 1], and specify the conductance equation ( 6) as In the next Section we compute the equilibrium point (V eq , W eq ) of the system (1)-( 5) for the chosen parameters.By equating the diffusion coefficient as it would occur in the diffusion approximation of [16] with the one in ( 8) at (V eq , W eq ) we will obtain σ * in terms of 1/ √ N , where N is the number of channels involved.
It can be shown by a coupling argument that also for varying V t will W t given by ( 8) stay bounded in (0, 1), since V t is bounded once it is started inside some interval [7].
In Fig. 2, the model defined by ( 1), ( 3)-( 5) and ( 8) is simulated for different values of σ * , where these can be thought of as corresponding to different total numbers of ion channels.

The linear approximation of the stochastic Morris-Lecar during quiescence
To identify the process of subthreshold oscillations, i.e. the dynamics close to the stable fixed point between firings, we analyze the linearized system around this point.Consider the system where the functions f (•), g(•) and h(•) are given by ( 1), ( 3)-( 5) and (8).
For the chosen parameter values given in Table 1, the deterministic system, obtained for h(•) = 0, has a unique locally stable equilibrium point (V eq , W eq ) given by and V eq is the solution to the equation f (V eq , W eq (V eq )) = 0, which cannot be solved analytically, but can be found numerically.The input current value I = 90µA/cm 2 is a typical value well inside the range of I where the deterministic dynamics has a stable limit point inside an unstable limit cycle as shown in Fig. 1A.The equilibrium point for I = 90µA/cm 2 is (V eq , W eq ) = (−26.6mV , 0.129).
In terms of the centered variables the system becomes dX t + V eq ), (X t + V eq ), (X t + W eq ) dt + h (X t + V eq ), (X t + W eq ) dB t .
We write t ) T , where T denotes transposition.Note that B (1) t does not enter the dynamics, but is introduced to ease the matrix notation, as will be clear in the following.When the noise is small and the process X t is started near the equilibrium point, x = (0, 0), we expect the dynamics to concentrate around the equilibrium point.A local approximation is obtained by linearizing ( 9)-( 10) around (0, 0).The diffusion term is approximated by setting X (1) t = X (2) t = 0 in the diffusion coefficients.The linearized system is where where σ = 0.034σ * .By evaluating the diffusion approximation of [16] at (V eq , W eq ) and equating to the above we obtain σ * = 1/ W eq (1 − W eq )N ≈ 3/ √ N .In the Appendix the matrix M is detailed.Solutions of (11) with G = 0 are given in terms of the eigenvalues of M which are complex conjugates and given by −λ ± ωi = −0.0094± 0.0803i Thus, near the equilibrium point the solution of (11), with σ = 0, is where C contains the initial conditions In Fig. 3 the solution of the deterministic model, (1)-( 5) with σ = 0, is compared to the linear approximation (13).

Identification of the stochastic process of quiescence
In this Section we identify the stochastic process defined by the linearized system (11) in the limit of small λ, i.e. under the condition λ ω.The deterministic system (13) has decaying oscillations, whereas for the stochastic system (11), the noise will prevent the decay of the oscillations.Can we describe the resulting process specifically?The answer is that, after a linear change of variables, this process can be approximated in distribution by a fixed matrix times a deterministic circular motion modulated by an OU process.
We follow the development in [2], where a first step is to transform the matrix M into a form which reveals the slow decay towards the equilibrium point and the fast oscillatory structure of the deterministic dynamics.Let Q be a 2 × 2 matrix such that  1)-( 5) with σ = 0 (black full drawn curves) is compared to the linear approximation (13)  A possible choice for Q is where C = Q −1 G.A further change of variables moves the rotation to form part of the diffusion coefficient of the linear stochastic system.We define is the counterclockwise rotation of angle s.Then by Ito's formula The infinitesimal covariance matrix in ( 14) is

Now define
where we have used that (m 11 + λ) 2 + ω 2 = −m 12 m 21 .Finally, we rescale Xt so that we can compare with a standardized two-dimensional OU process.Let Relation (15) becomes where Bt = √ λB t/λ is another standard two-dimensional Brownian motion.The following Theorem from [2] allows us to approximate the process U t given by (17), by a two-dimensional OU process with independent coordinates.Theorem 4.1.For each fixed t * > 0 and x ∈ IR 2 the distribution of {U t : 0 ≤ t ≤ t * } given by (17) with U 0 = x converges as λ/ω → 0 to the distribution of the standardized two-dimensional OU process {S t : 0 ≤ t ≤ t * } generated by Here S t follows a normal distribution, S t ∼ N S 0 e −t , 1  2 (1 − e −2t )I , where I is the 2 × 2 identity matrix.The proof of this Theorem uses a martingale problem convergence argument and involves the notion of stochastic averaging, where fast oscillations integrate out revealing the remaining structure determined by slower oscillations.Another result of this type obtained by a different method, called multiscale analysis, is in [17].
Thus, the process U t is approximated by S t if λ ω.In our case λ is one order of magnitude smaller than ω.
Putting together the transformations and the final approximation we have, in the sense of stochastic process distributions, Let us denote by X a t the stochastic process on the right hand side of (18), i.e.
To get a sense of how closely the process X a t approximates the dynamics of the ML process in a neighborhood of (V eq , W eq ) we compare their power spectral densities, as well as that of the solution of the linearized system (11).The spectral density of X a t and that ofX t satisfying (11) can be calculated explicitly using the power spectrum formula of [10] for linear diffusions of the form (11).In fact X a t is such a diffusion: the effect of the stochastic averaging can be seen as replacing C from ( 14) by a multiple of the identity in the system (14), so the approximation to X satisfies d Xa t = A Xa t dt + τ dB t , where τ is given by ( 16).If we transform this equation by X a t = Q Xa t , we see that X a t satisfies The spectral density of the first coordinate of X a is whereas the spectral density of the first coordinate of the linearized system, (11), is In Fig. 4 the theoretical spectral densities for the two approximations are plotted, together with the estimated spectral density of the quiescent process from simulations of the stochastic ML model ( 1),( 3)-( 5) and ( 8).The spectral density is estimated by averaging over at least 20 estimates from paths started at 0 of at least 450 ms of subthreshold fluctuations, and scaled to have the same maximum as the theoretical spectral density from (20).The averaging is done to reduce the large variance connected with spectral density estimation, avoiding any smoothing.Thus, the estimator is approximately unbiased, see also [8] where this approach is treated.The estimation is done for σ * = 0.03, 0.05 and 0.1.For higher noise, the lengths of subthreshold fluctuations between spikes are too short to reliably estimate the spectral density.Moreover, σ * = 0.1 corresponds to a number of ion channels N ≈ 900, which can be considered a minimum acceptable number for the diffusion approximation to be relevant.The value σ * = 0.03 corresponds to N ≈ 10, 000.Remember that σ = 0.034σ * , see (12).The approximations are only acceptable for small noise, which is expected, since larger noise brings the process to areas further away from the fixed point, where non-linearities become increasingly important.

Reconstructing the stochastic ML firing mechanism
In this Section we construct a firing mechanism matching that of the stochastic ML neuron.In Section 6 we will define a new LIF-type process by combining this firing mechanism with the radial OU process.This new model will, for small σ, have an ISI distribution similar to that of the ML.
Firing in model ( 1), ( 3)-( 5) and ( 8) occurs when the stochastic dynamics shifts from a path circulating the stable equilibrium, modulated by an OU, to a noisy circuiting of the stable limit cycle.This shift happens, roughly, when the orbit passes from the inside to the outside of the unstable limit cycle.When the orbit comes close to the unstable limit cycle, it will follow this limit cycle for a short time, and then escape either to the inside, i.e. continue its subthreshold oscillations, or to the outside and a spike will occur.This understanding is not accurate enough to be implemented as a firing scheme for the radial OU process (27), as we discuss further in Section 7. Hence, we embed the process X a defined by (19) in the stochastic ML model by constructing a firing mechanism mimicking that of the ML itself.It is clear that in the ML model, starting inside the unstable limit cycle, a spike will occur with increasing probability, the further away the process is from the fixed point.In order to construct a firing mechanism matching that of ML, we will estimate, from simulations, the conditional probability that the ML fires, given that the trajectory of the ML crosses the line L = {(v, w) : v = V eq , w < W eq }.We computed estimates from simulated data using crossings of the line L as follows.
For a given value of σ * and distance l from the fixed point, a short trajectory starting in (V eq , W eq − l) was simulated from model ( 1), ( 3)-( 5) and ( 8), and it was registered whether firing occurred in the first cycle of the stochastic path around (V eq , W eq ).Firing was defined by the path crossing the line v = 0, which is well above the largest level inside the unstable limit cycle, see Fig. 1B.This was repeated 1000 times, and estimates of the conditional probability of spiking, p(l, σ * ), were computed as the frequency of the trajectories where firing occurred.The procedure was repeated for l = l i = iδ, i = 1, . . ., 25, where δ is the distance to the stable limit cycle divided by 20.In this way a grid of possible l values was covered, starting from l = 0 at the fixed point, where the probability of firing is close to zero, to a point on L below the stable limit cycle, where the probability of firing is close to one.The estimation was, furthermore, repeated for σ * = 0.01 to 0.08 in steps of 0.01.
For each fixed σ * , the estimates of the conditional probability appear to depend in a sigmoidal way on the distance from the fixed point.We assumed the conditional firing probability to be of the form The parameters α and β were estimated using non-linear regression of the 25 estimates of p(l i ; σ * ) on l.In Fig. 5A these parametric estimates are plotted, as well as the individual nonparametric estimates p for σ * = 0.02, 0.05 and 0.08.We see that the family of estimates, p, fits the hypothetic curve quite well for each value of σ * .Regression estimates are reported in Table 2.Note that α is the distance along L from W eq at which the conditional probability of firing equals one half.For all values of σ * , the estimate of α is close to the distance along L between W eq and the unstable limit cycle, which equals 0.0172.In other words, the probability  q q q q q q q q q q q q q q q q q q q q q q q q q 0.005 0.015 0.025   21).The dashed line indicates where the unstable limit cycle crosses L, the full drawn line where the stable limit cycle crosses L. (B) The fitted curves in the transformed space for σ * = 0.02, 0.03, 0.04, 0.05, 0.06, 0.07 and 0.08 (right to left), as a function of the distance from the fixed point in the transformed coordinates.The crosses and boxed crosses indicate the crossing of the unstable and stable limit cycles of L, respectively, which depend on σ = 0.034σ * . of firing, if the path starts at the intersection of L with the unstable limit cycle, is about 1/2.The parameter β indicates the width of a band around α where the conditional probability essentially changes.For instance, if l ∈ α ± β then p(l) ∈ (0.27, 0.73), if l ∈ α ± 2β then p(l) ∈ (0.12, 0.88).As expected, the estimate of β increases with increasing σ * , and for small noise the conditional probability approaches a step function since the process is mostly dominated by the drift.A step function would correspond to the firing being represented by a first-passage time of a fixed threshold.Note though that β is approximately proportional to σ * , and thus, as we said earlier and will see in the following, a fixed threshold at the crossing of the unstable limit cycle does not reproduce the desired spiking characteristics.
In order to simplify the construction in Section 6 of a LIF model which, together with a firing rule, behaves like the stochastic ML, we will change coordinates as follows.Observe that ( 19) can be written so for fixed t, √ λQ −1 X a t /τ is the clockwise rotation by angle ωt of the orthogonal pair (S λt ).We define a transformation of the space (v, w) by centering at (V eq , W eq ) and normalizing as in (22).Let be the coordinates of the transformed space.In the new coordinates our process is simplified to a rotation modulated by a standard two-dimensional OU process with independent components.
The transformation depends on σ = 0.034σ * , namely, the transformed unstable limit cycle becomes smaller with increasing noise, through the value of τ given in (16).This is exactly what is causing a higher firing probability for larger σ * .The line L will in the transformed space be for l ≥ 0. A distance l will thus transform to a distance r = ( √ 2λ/σ)l, and the conditional probability of firing (21) transforms to where α * = α √ 2λ/σ and β * = β √ 2λ/σ.The fitted curves of (24) for σ * = 0.02 − 0.08, as a function of the distance from the fixed point in the transformed coordinates are given in Fig. 5B, with indication of the crossings of the unstable and stable limit cycles, respectively, which now depend on σ.Note that in the transformed space, the width of the band where the conditional probability is essentially different from 0 or 1 is nearly constant, see Table 2. From here on we use the coordinates defined by (23).

Construction of a leaky-integrate-and-fire model with ML firing statistics
The simpler stochastic LIF models sacrifice realism for mathematical tractability [4,11].In these models, a neuron is characterized by a single stochastic differential equation describing the evolution of neuronal membrane potential depending on time, where X t corresponds to V t in the ML model, together with a threshold firing rule, In this Section we define a LIF model which does not make this compromise, using the result of Section 4 and the firing mechanism defined in Section 5.
The distance of the approximate process √ λQ −1 X a t /τ of ( 22) from the point (0, 0) at time t is given by the modulus of the two-dimensional standardized OU process S λt .The modulus of S λt at time t is given by the process which is a standard radial OU process with two degrees of freedom.It has state space (0, ∞), and solves the stochastic differential equation see e.g.[3].We define a new LIF process by ( 27), and firing mechanism derived from (24).
After each firing, we will reset the time to 0 and assume the process reset to 0, i.e.R 0 = 0, corresponding to S 0 = (0, 0) and (V 0 , W 0 ) = (V eq , W eq ).By Ito's formula, the process Y u = R 2 u satisfies the stochastic differential equation and is thus a square-root process, see e.g.[5], also called a Feller or a Cox-Ingersoll-Ross process.This process is ergodic, and its stationary distribution is the exponential distribution with mean one.It follows that the stationary distribution of R u has density f (r) = 2re −r 2 on (0, ∞), i.e. it follows a Rayleigh distribution.The transition density of Y u starting at y 0 at time 0, is a non-central χ 2 -distribution with two degrees of freedom and non-centrality parameter δ(u, y 0 ) = 2y 0 e −2u /(1 − e −2u ).Then 2Y u /(1 − e −2u ) follows the standard non-central χ 2distribution F χ 2 (2y/(1 − e −2u ), 2, δ(u, y 0 )).It is particularly simple because of the integer degrees of freedom.Transforming to the radial OU we obtain the transition density of R u starting at s at time 0 where I 0 (x) = 1 π π 0 e x cos θ dθ is the modified Bessel function of the first kind of index 0. Writing the two-dimensional process S u in polar coordinates, R u and θ u , where θ u is the angle at time u to the positive part of the first coordinate, we find that the modulus and the angle are independent, and that θ u is uniformly distributed on (0, 2π).This can e.g.be seen from the fact that S Let T denote the firing time random variable.We want to compute the density of the distribution of T , and for this we find it convenient to express this density in terms of the conditional hazard rate, This function is the density of the conditional probability, given the position on L is r at time t, of a spike occuring in the next small time interval, given that it has not yet occurred.
From standard results from survival analysis, see e.g.[1], we obtain The unconditional distribution is then given by where E(•) denotes expectation with respect to the distribution of R. The density is thus The firing is defined to be initiated from L, and on average the process crosses L every 2π/ω = 78.2time units.Using (24), the estimated conditional probability of firing given the position on L is r, which by definition does not depend on t, we estimate the hazard rate as Note that it is bounded.This is not realistic, since a very large value of r should cause immediate firing.
In [21] a firing rule with unbounded hazard rate was proposed, and in [15] it was shown to fit well to experimental data.Therefore, we will also see how our model performs if we use in the firing mechanism a hazard rate of the form for α, β > 0. Like before, α plays the role of a threshold, and β gives the width of the threshold region.When β → 0, the firing rule converges to a fixed threshold crossing.To estimate α and β in (33), we simulated 1000 spike times from the ML.The cumulative hazard A(t) = t 0 α(t) was then estimated from the simulated spike times by the standard empirical Nelson-Aalen estimator.The theoretical cumulative hazard using ( 27) and (33) can be calculated as where we have used the density f λs (r, 0) given in (29).Here, g(s) = √ 1 − e −2λs /β, and Φ(•) is the standard normal cumulative distribution function.Then, α and β were estimated by the least square distance between (34) and the estimated cumulative hazard from the simulated spike times.For σ * = 0.05 the estimates were α = 6.31 and β = 0.76.
The final model is where µ(R u− , du) is a Poisson measure with intensity α(R u− ), and R u− denotes the left limit of R u .Here, α(•) is either given by (32) or (33).The jump size is −R u− , thus giving the reset to 0 at spike times.A reasonable alternative to the soft threshold firing mechanism used here would be to use the firing rule defined by a threshold as in the classical LIF models, equation (26).A natural choice of threshold would be where the LIF process reaches a level corresponding to the unstable limit cycle.In fact, according to our estimates in Fig. 5 and Table 2, the firing probability of the ML at this threshold is around 1/2.However, the ISI distribution estimated from simulations using a hard threshold at the unstable limit cycle is shifted towards larger times, relative to the ML ISI distribution.This happens because the process might cycle many times inside the unstable limit cycle, so even if the probability of spiking in a single cycle is small, the total probability is not negligible.This is lost when only a hard threshold is considered.Instead we chose the threshold value such that the mean of the ML ISI distribution and the mean of the LIF ISI distribution were the same.In [12], the mean of T from (26) with X t = R t started at R 0 = 0 is given using a hypergeometric function, The average of the 1000 ML firing times for σ * = 0.05 was 447.Equating with (36) gives a value S = 2.97 for the hard threshold.Note that this is much smaller than the estimated α from (33).

Comparison of firing statistics
One of the major issues in computational neuroscience is to determine the ISI distribution.We therefore simulated the ML model given by ( 1), ( 3)-( 5) and ( 8) until spiking, and thereafter reset to the fixed point.This was done 1000 times, and the time of the firing was recorded.The ISI distribution from our approximate model is given by the density (31), or equivalently, from the survival function (30).Due to the law of large numbers and since we know the exact distribution of R u , for fixed t we can numerically determine (31) up to any desired precision by choosing n and M large enough through the expression λt ) are M realizations of R iλt/n , i = 0, 1, . . ., n, and the integral has been approximated by the trapezoidal rule.The hazard rate is either given by (32) or (33).
The results are illustrated in Figure 6 for σ * = 0.05, using M = 1000.The estimated ISI distributions from our approximate model with both firing mechanisms compare well with the estimated ISI distribution of ML reset to 0 after firings.On the contrary, the hard threshold does not reproduce the ISI distribution well, e.g. the right tail is too heavy.This is because the probability of firing during low subthreshold activity is set to 0, whereas we have seen it is not.

Discussion
A stochastic LIF model constructed with a radial OU process and firing mechanism of either logistic or exponential type has been shown to mimic the ISI statistics of a ML neuron model.It captures subthreshold dynamics, not of the membrane potential alone, but of a combination of the membrane potential and ion channels.This construction will allow us to answer several questions about ML models, which have been accessible only for LIF models, even though the latter have less biological motivation.
An example of such a question would be: Using ISI experimental data, the noise standard deviation σ can be estimated [18].In principle, this should also be possible from our new LIF model, even though we use a soft threshold.This will give an estimate of N , the number of ion channels involved, through the relation (σ * ) 2 ≈ 9/N .
A question we have not explored is: what is the best way to restart our new LIF model?In our simulations we restarted both our LIF and the ML at the fixed point of the ML.However, an uninterrupted stochastic ML produces continuous paths as in Black curve is estimated using (32), gray curve is estimated using (33), dashed curve is estimated using a fixed threshold, (26).means traversing the large stable limit cycle, possibly several times, they reenter a neighborhood of the fixed point from its edge.A further refinement of our LIF model will be obtained by introducing a reentry mechanism, which mimics this aspect of the ML.

2. 1 .
The stochastic Morris-Lecar model with channel noise It has long been known that the opening and closing of ion channels is an important part of neuron function.Channel activity is summarized, even in the comparatively detailed HH

FIGURE 2 :
FIGURE 2: Time series plots (black curves) of the stochastic Morris-Lecar model for different noise levels started inside the unstable limit cycle, but not at the fixed point.Upper left: σ * = 0.02, upper right: σ * = 0.03, lower left: σ * = 0.05, lower right: σ * = 0.1.Note different scales, in the upper left panel there is no firing.The magenta curves are the deterministic model, σ * = 0.

FIGURE 5 :
FIGURE 5: Conditional probability of spiking when crossing the line L = {(v, w) : v = Veq, w < Weq} for different values of σ * .(A) Original space.The circles, plus's and stars are individual nonparametric estimates obtained using σ * = 0.02, 0.05 and 0.08, respectively, with the fitted curves on top given by (21).The dashed line indicates where the unstable limit cycle crosses L, the full drawn line where the stable limit cycle crosses L. (B) The fitted curves in the transformed space for σ * = 0.02, 0.03, 0.04, 0.05, 0.06, 0.07 and 0.08 (right to left), as a function of the distance from the fixed point in the transformed coordinates.The crosses and boxed crosses indicate the crossing of the unstable and stable limit cycles of L, respectively, which depend on σ = 0.034σ * .

uu
are independent normal with mean 0 and equal variances.Thus, for fixed u, S is standard Cauchy distributed and θ u = arctan(S

Figure 6 :
Figure 6: Distribution of firing times for σ * = 0.05.The histogram is based on 1000 simulated firing times from the ML model, the vertical dotted line is the average.Curves are estimates of the probability density, equation (37).Black curve is estimated using (32), gray curve is estimated using (33), dashed curve is estimated using a fixed threshold, (26).

TABLE 1 :
Variables and parameter values used in the Morris-Lecar model Ca = 4.4 µS/cm 2 Maximal conductance associated with Ca 2+ current g K = 8 µS/cm 2 Maximal conductance associated with K + current g L = 2 µS/cm 2 Conductance associated with leak current V Ca = 120 mV Reversal potential for Ca 2+ current V K = -84 mV Reversal potential for K + current V L = -60 mV Reversal potential for leak current C = 20 µF/cm 2 Membrane capacitance φ = 0.04 1/ms Rate scaling parameter I = 90 µA/cm 2 Input current that a K + ion channel is open at time t.The non-linear model equations are

TABLE 2 :
Estimates of regression parameters for p(•) in the original space (first two rows), and in the transformed coordinates (last two rows).