{'fraction_non_alphanumeric': 0.038742432058208674, 'fraction_numerical': 0.010208237178617002, 'mean_word_length': 4.626947754353804, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stochastic Representations of Ion Channel Kinetics and Exact Stochastic Simulation of Neuronal Dynamics

In this paper we provide two representations for stochastic ion channel kinetics, and compare the performance of exact simulation with a commonly used numerical approximation strategy. The first representation we present is a random time change representation, popularized by Thomas Kurtz, with the second being analogous to a"Gillespie"representation. Exact stochastic algorithms are provided for the different representations, which are preferable to either (a) fixed time step or (b) piecewise constant propensity algorithms, which still appear in the literature. As examples, we provide versions of the exact algorithms for the Morris-Lecar conductance based model, and detail the error induced, both in a weak and a strong sense, by the use of approximate algorithms on this model. We include ready-to-use implementations of the random time change algorithm in both XPP and Matlab. Finally, through the consideration of parametric sensitivity analysis, we show how the representations presented here are useful in the development of further computational methods. The general representations and simulation strategies provided here are known in other parts of the sciences, but less so in the present setting.


Introduction
Fluctuations in membrane potential arise in part due to stochastic switching in voltagegated ion channel populations [11,28,47]. We consider a stochastic modeling, i.e. master equation, framework [5,10,12,29,48] for neural dynamics, with noise arising through the molecular fluctuations of ion channel states. We consider model nerve cells that may be represented by a single isopotential volume surrounded by a membrane with capacitance C > 0. Mathematically, these are hybrid stochastic models which include components, for example the voltage, that are continuous and piecewise differentiable and components, for example the number of open potassium channels, that make discrete transitions or jumps.
These models are typically described in the neuroscience literature by providing an ODE governing the absolutely continuous portion of the system, which is valid between jumps of the discrete components, and a chemical master equation providing the dynamics of the probability distribution of the jump portion, which itself depends upon the solution to the ODE. In this paper, we take a different perspective. We will introduce here two pathwise stochastic representations for these models that are similar to Itô SDEs or Langevin models. The difference between the models presented here and Langevin models is that here the noise arrises via stochastic counting processes as opposed to Brownian motions. These representations give a different type of insight into the models than master equation representations do, and, in particular, they imply different exact simulation strategies. These strategies are well known in some parts of the sciences, but less well known in the current context [5]. 1 From a computational standpoint the change in perspective from the master equation to pathwise representations is useful for a number of reasons. First, the different representations naturally imply different exact simulation strategies. Second, and perhaps more importantly, the different representations themselves can be utilized to develop new, highly efficient, computational methods such as finite difference methods for the approximation of parametric sensitivities, and multi-level Monte Carlo methods for the approximation of expectations [2,4]. Third, the representations can be utilized for the rigorous analysis of different computational strategies and for the algorithmic reduction of models with multiple scales [3,7].
We note that the types of representations and simulation strategies highlighted here are known in other branches of the sciences, especially operations research and queueing theory [19,22], and stochastically modeled biochemical processes [1,5]. See also [36] for a mathematical treatment of such stochastic hybrid systems and some corresponding approximation algorithms. However, the representations are not as well known in the context of computational neuroscience and as a result approximate methods for the simulation of sample paths including (a) fixed time step methods, or (b) piecewise constant propensity algorithms, are still utilized in the literature in situations where there is no need to make such approximations. Thus, the main contributions of this paper are: (i) the formulation of the two pathwise representations for the specific models under consideration, (ii) a presentation of the corresponding exact simulation strategies for the different representations, and (iii) a comparison of the error induced by utilizing an approximate simulation strategy on the Morris-Lecar model. Moreover, we show how to utilize the different representations in the development of methods for 1 For an early statement of an exact algorithm for the hybrid case in a neuroscience context see ( [9], Equations 2-3). Strassberg and DeFelice further investigated circumstances under which it is possible for random microscopic events (single ion channel state transitions) to generate random macroscopic events (action potentials) [45] using an exact simulation algorithm. Keener and Newby used an exact algorithm in a recent study of channel noise dependent action potential generation in the Morris-Lecar model [24,34]. parametric sensitivity analysis.
The outline of the paper is as follows. In Section, 2 we develop the different representations and provide the corresponding numerical methods for the two different representations. In Section 3, we present a stochastic representation for a particular conductance based model, the planar Morris-Lecar model, with a single ion channel treated stochastically. Here, we illustrate the numerical strategy on this example and provide in the appendix both XPP and Matlab code for this model. In Section 4, we present an example of a conductance based model with more than one stochastic gating variable, namely the Morris-Lecar model with both ion channel types (calcium channels and potassium channels) treated stochastically. To illustrate both the strong and weak divergence between the exact and approximate algorithms, in Section 5 we compare trajectories and histograms generated by the exact algorithms and the piecewise constant approximate algorithm. In Section 6, we show how to utilize the different representations presented here in the development of methods for parametric sensitivity analysis, which is a powerful tool for determining parameters to which a system output is most responsive. In Section 7, we provide conclusions and discuss avenues for future research.

Two stochastic representations
In Section 2.1, we heuristically present the random time change representation for the relevant models, whereas in Section 2.2 we present a "Gillespie" representation. In each section, we provide the corresponding numerical simulation strategies. See [5] for a development of the representations in the biochemical setting, and see [26,27] for a more general rigorous mathematical development. For a more recent addition to the literature, see also [36] in which the authors, among other things, explicitly prove that the models satisfying the representations below are Markov processes satisfying the corresponding master equation.

Random time change representations
We begin with an example. Consider a model of a system that can be in one of two states, A or B, which can represent "closed" and "open," respectively. We model the dynamics of the system by assuming that the dwell times in states A and B are exponential random variables with parameters α > 0 and β > 0, respectively. A graphical representation for this model is The usual formulation of the stochastic version of this model is to note that the probability that a closed channel opens in the next increment of time ∆s is α∆s + o(∆s), whereas the probability that an open channel closes is β∆s + o(∆s). This type of stochastic model is often described mathematically by providing the "chemical master equation," which for (1) is simply where p x 0 (x, t) is the probability of being in state x ∈ {A, B} at time t given an initial condition of state x 0 . Note that the chemical master equation is a linear ODE governing the dynamical behavior of the probability distribution of the model, and does not provide a stochastic representation for a particular realization of the process. In order to construct such a path-wise representation, let R 1 (t) be the number of times the transition A → B has taken place by time t and, similarly, let R 2 (t) be the number of times the transition B → A has taken place by time t. We let X 1 (t) ∈ {0, 1} be one if the channel is closed at time t, and zero otherwise, and let X 2 (t) = 1 − X 1 (t) take the value one if and only if the channel is open at time t. Then, denoting X(t) = (X 1 (t), X 2 (t)) T , we have We now consider how to represent the counting processes R 1 , R 2 in a useful fashion, and we do so with unit-rate Poisson processes as our mathematical building blocks. We recall that a unit-rate Poisson process is constructed in the following manner. Let be independent exponential random variables with a parameter of one. Then, let τ 1 = e 1 , τ 2 = τ 1 + e 2 , · · · , τ n = τ n−1 + e n , etc. The associated unit-rate Poisson process, Y (s), is simply the counting process determined by the number of points , that come before s ≥ 0. For example, if we let "x" denote the points τ n in the image below then Y (s) = 6. If instead of moving at constant rate, s, along the horizontal axis, we move instead at rate λ(s), then the number of points observed by time s is Y s 0 λ(r)dr . Further, from basic properties of exponential random variables, whenever λ(s) > 0 the probability of seeing a jump within the next small increment of time ∆s is Thus, the propensity for seeing another jump is precisely λ(s).
Combining all of the above implies that we can represent R 1 and R 2 via and so a pathwise representation for the stochastic model (1) is where Y 1 and Y 2 are independent, unit-rate Poisson processes. Slightly more generally, suppose now that X 1 (0) + X 2 (0) = N ≥ 1. For example, perhaps we are now modeling the number of open and closed ion channels out of a total of N , as opposed to simply considering a single such channel. Suppose further that the propensity, or rate, at which ion channels are opening can be modeled as and the rate at which they are closing has propensity where α(t), β(t) are non-negative functions of time, perhaps being voltage dependent. That is, suppose that for each i ∈ {1, 2}, the conditional probability of seeing the count- The expression analogous to (2) is now Having motivated the time dependent representation for the simple model above, we turn to the more general context. We now assume a jump model consisting of d chemical constituents (or ion channel states) undergoing transitions determined via M > 0 different reactions. For example, in the toy model above, the chemical constituents were {A, B}, and so d = 2, and the reactions were A → B and B → A, giving M = 2. We suppose that X i (t) determines the value of the ith constituent at time t, so that X(t) ∈ Z d , and that the propensity function of the kth reaction is λ k (t, X(t)). We further suppose that if the kth reaction channel takes place at time t, then the system is updated according to addition of the reaction vector ζ k ∈ Z d , The associated pathwise stochastic representation for this model is where the Y k are independent unit-rate Poisson processes. The chemical master equation for this general model is where P X 0 (x, t) is the probability of being in state x ∈ Z d ≥0 at time t ≥ 0 given an initial condition of X 0 .
When the variable X ∈ Z d represents the randomly fluctuating state of an ion channel in a single compartment conductance based neuronal model, we include the membrane potential V ∈ R as an additional dynamical variable. In contrast with neuronal models incorporating Gaussian noise processes, here we consider the voltage to evolve deterministically, conditional on the states of one or more ion channels. For illustration, suppose we have a single ion channel type with state variable X. Then, we supplement the pathwise representation (4) with the solution of a differential equation obtained from Kirchoff's current conservation law: Here, g o i is the conductance of an individual channel when it is the i th state, for 1 ≤ i ≤ d. The sum gives the total conductance associated with the channel represented by the vector X; the reversal potential for this channel is the constant V X . The term I V (V ) captures any deterministic voltage-dependent currents due to other channels besides channel type X, and I app represents a time-varying, deterministic applied current. In this case the propensity function will explicitly be a function of the voltage and we may replace λ k (s, X(s)) in (4) with λ k (V (s), X(s)). If multiple ion channel types are included in the model then, provided there are a finite number of types each with a finite number of individual channels, the vector X ∈ Z d represents the aggregated channel state. For specific examples of handling a single or multiple ion channel types, see Sections 3 and 4, respectively.

Simulation of the representation (4)-(5)
The construction of the representation (4) provided above implies a simulation strategy in which each point of the Poisson processes Y k , denoted τ n above, is generated sequentially and as needed. The time until the next reaction that occurs past time T is simply where τ k T is the first point associated with Y k coming after T 0 λ k (s, X(s))ds. The reaction that took place is indicated by the index at which the minimum is achieved. See [1] for more discussion on this topic, including the algorithm provided below in which T k will denote the value of the integrated intensity function t 0 λ k (s, X(s))ds and τ k will denote the first point associated with Y k located after T k .
All random numbers generated in the algorithm below are assumed to be independent.
Algorithm 1 (For the simulation of the representation (4)-(5)). .  (5) forward in time until one of the following equalities hold:

Numerically integrate
4. Set ∆ = min k {∆ k } and let µ be the index where the minimum was achieved.

Return to step 3 or quit.
Note that the above algorithm relies on us being able to calculate a hitting time for each of the T k (t) = t 0 λ k (s, X(s))ds exactly. Of course, in general this is not possible. However, making use of any reliable integration software will almost always be sufficient.

Gillespie representation
There are multiple alternative representations for the general stochastic process constructed in the previous section, with a "Gillespie" representation probably being the most useful in the current context. Following [5], we let Y be a unit rate Poisson process and let {ξ i , i = 0, 1, 2 . . . } be independent, uniform (0, 1) random variables that are independent of Y . Set Then the stochastic process (X, V ) defined via (6)-(8) is a Markov process that is equivalent to (4)-(5), see [5]. An intuitive way to understand the above is by noting that R 0 (t) simply determines the holding time in each state, whereas (7) simulates the embedded discrete time Markov chain (sometimes referred to as the skeletal chain) in the usual manner. Thus, this is the representation for the "Gillespie algorithm" [17] with time dependent propensity functions [1].

Simulation of the representation (6)-(8)
Simulation of (6)-(8) is the analog of using Gillespie's algorithm in the time-homogeneous case. All random numbers generated in the algorithm below are assumed to be independent.

1.
Initialize: set the initial number of molecules of each species, X. Set the initial voltage value V . Set t = 0.
5. Return to step 2 or quit.

Morris-Lecar
As a concrete illustration of the exact stochastic simulation algorithms, we will consider the well known Morris-Lecar system [37], developed as a model for oscillations observed in barnacle muscle fibers [32]. The deterministic equations, which correspond to an appropriate N → ∞ limit of the system, constitute a planar model for the evolution of the membrane potential v(t) and the fraction of potassium gates, n ∈ [0, 1], that are in the open or conducting state. In addition to a hyperpolarizing current carried by the potassium gates, there is a depolarizing calcium current gated by a rapidly equilibrating variable m ∈ [0, 1]. While a fully stochastic treatment of the Morris-Lecar system would include fluctuations in this calcium conductance, for simplicity in this section we will treat m as a fast, deterministic variable in the same manner as in the standard fast/slow decomposition, which we will refer to here as the planar Morris-Lecar model [13]. See Section 4 for a treatment of the Morris-Lecar system with both the potassium and calcium gates represented as discrete stochastic processes. The deterministic or mean field equations for the planar Morris-Lecar model are: The kinetics of the potassium channel may be specified either by the instantaneous time constant and asymptotic target, τ and n ∞ , or equivalently by the per capita transition rates α and β. The terms m ∞ , α, β, n ∞ and τ satisfy where for convenience we define for which the deterministic system has a stable limit cycle. For smaller values of the applied current (e.g. I app = 75) the system has a stable fixed point, that loses stability through a subcritical Hopf bifurcation as I app increases ( [13], §3).
In order to incorporate the effects of random ion channel gating, we will introduce a finite number of potassium channels, N tot , and treat the number of channels in the open state as a discrete random process, 0 ≤ N (t) ≤ N tot . In this simple model, each potassium channel switches between two states -closed or open -independently of the others, with voltage-dependent per capita transition rates α and β, respectively. The entire population conductance ranges from 0 to g o K N tot , where g o K = g K /N tot is the single channel conductance, and g K is the maximal whole cell conductance. For purposes of illustration and simulation we will typically use N tot = 40 individual channels. 2 Our random variables will therefore be the voltage, V ∈ (−∞, +∞), and the number of open potassium channels, N ∈ {0, 1, 2, · · · , N tot }. The number N tot ∈ N is taken to be a fixed parameter. We follow the usual capital/lowercase convention in the stochastic processes literature: N (t) and V (t) are the random processes and n and v are values they might take. In the random time change representation of Section 2.1, the opening and closing of the potassium channels are driven by two independent, unit rate Poisson processes, Y open (t) and Y close (t).
The evolution of V and N is closely linked. Conditioned on N having a specific value, say N = n, the evolution of V obeys a deterministic differential equation, Although its conditional evolution is deterministic, V is nevertheless a random variable. Meanwhile, N evolves as a jump process, i.e. N (t) is piecewise-constant, with transitions N → N ± 1 occurring with intensities that depend on the voltage V . Conditioned on V = v, the transition rates for N are Graphically, we may visualize the state space for N in the following manner: with the nodes of the above graph being the possible states for the process N , and the transition intensities located above and below the transition arrows. Adopting the random time change representation of Section 2.1 we write our stochastic Morris-Lecar system as follows (cf. equations (9-10)): 4 Models with more than one channel type The algorithms presented above readily generalize to systems with more than one type of stochastic channel, such as the Hodgkin-Huxley system [23,43]. As an illustration, we include code in Appendices A.2.1 and A.2.2, for XPP and Matlab, respectively, implementing a stochastic version of the original Morris-Lecar model, in a three dimensional phase space with both the calcium and potassium channels taken to be discrete, using the random time change algorithm ((4)- (5)).

Random Time Change Representation for the Morris-Lecar Model with Two Stochastic Channel Types
In Morris and Lecar's original treatment of voltage oscillations in barnacle muscle fiber [33] the calcium gating variable m is included as a dynamical variable. The full (deterministic) equations have the form: Here, rather than setting m to its asymptotic value m ∞ = α m /(α m + β m ), we allow the number of calcium gates to evolve according to (27). The planar form (Equations 9-10) is obtained by observing that m approaches equilibrium significantly more quickly than n and v. Using standard arguments from singular perturbation theory [37,38], one may approximate certain aspects of the full system (25)-(27) by setting m to m ∞ (v), and replacing F (v, n, m) and G(v, n, m) with f (v, n) = F (v, n, m ∞ (v)) and g(v, n) = G(v, n, m ∞ (v)), respectively. This reduction to the slow dynamics leads to the planar model (9)-(10).
To specify the full 3D equations, we introduce where the voltage falls along the activation curve for channel type x, relative to its half-activation point (v a for calcium and v c for potassium) and its slope (reciprocals of v b for calcium and v d for potassium). The per capita opening rates α m , α n and closing rates β m , β n for each channel type are given by with parameters v a = −1.
The asymptotic open probabilities for calcium and potassium are given, respectively, by the terms m ∞ , n ∞ , and the time constants by τ m and τ n . These terms satisfy the relations τ n (v) = 1/ (φ cosh(ξ n /2)) .
Assuming a population of M tot calcium gates and N tot potassium gates, the voltage evolves as according to the sum of the applied, leak, calcium, and potassium currents: while the number of open M and N remain constant except for unit increments and decrements. Figure 2 shows the results of the Matlab implementation for the 3D Morris-Lecar system, with both the potassium and calcium channel treated discretely, using Algorithm 1 (the random time change algorithm, (4)- (5)). Here M tot = N tot = 40 channels, and the applied current I app = 100 puts the deterministic system at a stable limit cycle close to a Hopf bifurcation.

Comparison of the Exact Algorithm with a Piecewise Constant Propensity Approximation.
Exact versions of the stochastic simulation algorithm for hybrid ion channel models have been known since at least the 1980s [9]. Nevertheless, the implementation one finds most often used in the literature is an approximate method in which the per capita reaction propensities are held fixed between channel events. That is, in step 3 of Algorithm 1 the integral leaving the remainder of the algorithm unchanged. Put another way, one generates the sequence of channel state jumps using the propensity immediately following the most recent jump, rather than taking into account the time dependence of the reaction propensities due to the continuously changing voltage. This piecewise constant propensity approximation is analogous, in a sense, to the forward Euler method for the numerical solution of ordinary differential equations. Figure 3 shows a direct comparison of pathwise numerical solutions obtained by the exact method provided in Algorithm 1 and the approximate forward method detailed above. In general, the solution of a stochastic differential equation with a given initial condition is a map from the sample space Ω to a space of trajectories. In the present context, the underlying sample space consists of one independent unit rate Poisson process per reaction. For the planar Morris-Lecar model a point in Ω amounts to fixing two Poisson processes, Y open and Y closed , to drive the transitions of the potassium channel. For the full 3D Morris-Lecar model we have four processes, The "exact" algorithm provides a numerical solution of the map from {Y k } 4 k=1 ∈ Ω and initial conditions (M 0 , N 0 , V 0 ) to the trajectory (M (t), N (t), V (t)). The approximate forward algorithm gives a map from the same domain to a different trajectory, (M (t),Ñ (t),Ṽ (t)). To make a pathwise comparison for the full Morris-Lecar model, therefore, we fix both the initial conditions and the four sample processes, and compare the resulting trajectories. Several features are evident in Figure 3. Both algorithms produce a sequence of noise-dependent voltage spikes, with similar firing rates. The trajectories (M, N, V ) and (M ,Ñ ,Ṽ ) initially remain close together, and the timing of the first spike (taken, e.g., as an upcrossing of V from negative to positive voltage) is similar for both algorithms. Over time, however, discrepancies between the trajectories accumulate. The timing of the second and third spikes is noticeably different, and before ten spikes have accumulated the spike trains have become effectively uncorrelated.
Even though trajectories generated by the exact and approximate algorithms diverge when driven by identical Poisson processes, the two processes could still generate sample paths with similar time-dependent or stationary distributions. That is, even though the two algorithms show strong divergence, they could still be close in a weak sense.
Given M tot calcium and N tot potassium channels, the density for the hybrid Markov process may be written ,Ṽ (t)) obtained using the forward approximation. Both algorithms were begun with identical initial conditions (M 0 , N 0 , V 0 ), and driven by the same four Poisson process streams Y 1 , · · · , Y 4 . Note the gradual divergence of the trajectories as differences between the exact and the forward approximate algorithms accumulate, demonstrating "strong" (or pathwise) divergence of the two methods. The exact and approximate trajectories diverge as time increases, even though they are driven by identical noise sources.
and obeys a master equation with initial condition ρ m,n (v, 0) ≥ 0 given by any (integrable) density such that v∈R m,n ρ m,n (v, 0) dv ≡ 1, and boundary conditions ρ → 0 as |v| → ∞ and ρ ≡ 0 for either m, n < 0 or m > M tot or n > N tot .
In contrast, the approximate algorithm with piecewise constant propensities does not generate a Markov process, since the transition probabilities depend on the past rather than the present values of the voltage component. Consequently they do not satisfy a master equation. Nevertheless it is plausible that they may have a unique stationary distribution. Figure 4 shows pseudocolor plots of the histograms viewed in the (v, n) plane, i.e. with entries summed over m, for Algorithm 1 ("Exact") and the approximate forward ("Euler") algorithm, with M tot = N tot = k channels for k = 1, 2, 5, 10, 20 and 40. The two algorithms were run with independent random number streams in the limit cycle regime (with I app = 100) for t max ≈ 200, 000 time units, sampled every 10 time units, which generated ≥ 17, 000 data points per histogram. For k < 5 the difference in the histograms is obvious at a glance. For k ≥ 10 the histograms appear increasingly similar. Figure 5 shows bar plots of the histograms projected on the voltage axis, i.e. with entries summed over m and n, for the same data as in Figure 4, with M tot = N tot = k channels ranging from k = 1 to k = 40. Again, for k ≤ 5, the two algorithms generate histograms that are clearly distinct. For k ≥ 20 they appear similar, while k = 10 appears to be a borderline case.
To quantify the similarity of the histograms we calculated the empirical L 1 difference between the histograms in two ways: first for the full (v, n, m) histograms, and then for the histograms collapsed to the voltage axis. Let ρ m,n (v) andρ m,n (v) denote the stationary distributions for the exact and the approximate algorithms, respectively (assuming these exist). To compare the two distributions we approximate the L 1 distance between them, i.e.
where v min and v max are chosen so that for all m, n we have F (v min , n, m) > 0 and F (v max , n, m) < 0. It is easy to see that such values of v min and v max must exist, since F (v, n, m) is a linear and monotonically decreasing function of v for any fixed pair (n, m), cf. equation (25). Therefore, for any exact simulation algorithm, once the voltage component of a trajectory falls in the interval v min ≤ v ≤ v max , it remains in this interval for all time. We approximate the voltage integral by summing over 100 evenly spaced bins along the voltage axis. Figure 6, Figure 6: L 1 differences between empirical histograms generated by exact and approximate algorithms. H j (v, n, m) is the number of samples for which V ∈ [v, v + ∆v), N = n and M = m, where ∆v is a discretization of the voltage axis into 100 equal width bins; j = 1 represents an empirical histogram obtained by an exact algorithm (Algorithm 1), and j = 2 represents an empirical histogram obtained by a piecewise-constant propensity approximation or "Euler" algorithm. Note the normalized L 1 difference of two histograms can range from 0 to 2. Top: Difference for the full histograms, Here v refers to binned voltage values. Bottom: Difference for the corresponding voltage histograms,

Coupling, variance reduction, and parametric sensitivities
It is relatively straightforward to derive the exact simulation strategies from the different representations provided here. Less obvious, and perhaps even more useful, is the fact that the random time change formalism also lends itself to the development of new computational methods that potentially allow for significantly greater computational speed, with no loss in accuracy. They do so by providing ways to couple two processes in order to reduce the variance, and thereby increase the speed, of different natural estimators.
For an example of a common computational problem where such a benefit can be had, consider the computation of parametric sensitivities, which is a powerful tool in determining parameters to which a system output is most responsive. Specifically, suppose that the intensity/propensity functions are dependent on some vector of parameters θ. For instance, θ may represent a subset of the system's mass action kinetics constants, the cell capacitance, or the underlying number of channels of each type. We may wish to know how sensitively a quantity such as the average firing rate or interspike interval variance depends on the parameters θ. We then consider a family of models V θ , X θ , parameterized by θ, with stochastic equations where f is some function, and all other notation is as before. Some quantities arising in neural models depend on the entire sample path, such as the mean firing rate and the interspike interval variance. Let g θ, V θ , X θ by a path functional capturing the quantity of interest. In order to efficiently and accurately evaluate the relative shift in the expectation of g due to a perturbation of the parameter vector, we would estimatẽ represents the ith path generated with a parameter choice of θ, and N is the number of sample paths computed for the estimation. That is, we would use a finite difference and Monte Carlo sampling to approximate the change in the expectation.
If the paths V θ are generated independently, then the variance of the estimator (38) fors is O(N −1 −2 ), and the standard deviation of the estimator is O(N −1/2 −1 ). This implies that in order to reduce the confidence interval of the estimator to a target level of ρ > 0, we require which can be prohibitive. Reducing the variance of the estimator can be achieved by coupling the processes V θ , so that they are correlated, by constructing them on the same probability space. This process is analogous to the construction of simulation paths according to the random time change algorithm and a piecewise constant propensity approximation algorithm used to illustrate strong divergence ( Figure 3). The representations we present here lead rather directly to schemes for reducing the variance of the estimator. We discuss different possibilities.
1. The common random number (CRN) method. The CRN method simulates both processes according to the Gillespie representation (6)-(8) with the same Poisson process Y and the same stream of uniform random variables {ξ i }. In terms of implementation, one simply reuses the same two streams of uniform random variables in Algorithm 2 for both processes.
2. The common reaction path method (CRP). The CRP method simulates both processes according to the random time change representation with the same Poisson processes Y k . In terms of implementation, one simply makes one stream of uniform random variables for each reaction channel, and then uses these streams for the simulation of both processes. See [35], where this method was introduced in the context of biochemical models.
3. The coupled finite difference method (CFD). The CFD method utilizes a "split coupling" introduced in [2]. Specifically, it splits the counting process for each of the reaction channels into three pieces: one counting process that is shared by X θ and X θ (and has propensity equal to the minimum of their respective intensities for that reaction channel), one that only accounts for the jumps of X θ , and one that only accounts for the jumps of X θ . Letting a ∧ b = min{a, b}, the precise coupling is where , X θ (s)), and where {Y k,1 , Y k,2 , Y k,3 } are independent unit-rate Poisson processes. Implementation is then carried out by Algorithm 1 in the obvious manner.
While the different representations provided in this paper imply different exact simulation strategies (i.e. Algorithms 1 and 2), those strategies still produce statistically equivalent paths. This is not the case for the methods for parametric sensitivities provided in this section. To be precise, each method constructs a coupled pair of processes (V θ , X θ ), (V θ , X θ ) , and the marginal processes V θ , X θ and V θ , X θ are all statistically equivalent no matter the method used. However, the covariance Cov (V θ (t), X θ (t)), (V θ (t), X θ (t)) can be drastically different. This is important, since it is variance reduction we are after, and for any component X j , Thus, minimizing variance is equivalent to maximizing covariance. Typically, the CRN method does the worst job of maximizing the covariance, even though it is the most widely used method [44]. The CFD method typically does the best job of maximizing the covariance, though examples exist in which the CRP method is the most efficient [2,44].

Discussion
We have provided two general representations for stochastic ion channel kinetics, one based on the random time change formalism, and one extending Gillespie's algorithm to the case of ion channels driven by time-varying membrane potential. We believe that the random time change representation (Algorithm 1, (4)-(5)) will be particularly useful to the computational neuroscience community as it allows for generalizations of computational methods developed in the context of biochemistry, in which the propensities depend upon the state of the jump process only. For example, variance reduction strategies for the efficient computation of first and second order sensitivities [2,6,35], as discussed in Section 6, and for the efficient computation of expectations using multilevel Monte Carlo [4,16] now become feasible. The random time change approach avoids several approximations that are commonly found in the literature. In simulation algorithms based on a fixed time step chemical Langevin approach, it is necessary to assume that the increments in channel state are approximately Gaussian distributed over an appropriate time interval [15,20,21,31,46]. However, in exact simulations with small membrane patches corresponding to M tot = 40 calcium channels, the exact algorithm routinely visits the states M (t) = 0 and M (t) = 40, for which the Gaussian increment approximation is invalid regardless of time step size. The main alternative algorithm typically found in the literature is the piecewise constant propensity or approximate forward algorithm [14,25,41]. However, this algorithm ignores changes to membrane potential during the intervals between channel state changes. As the sensitivity of ion channel opening and closing to voltage is fundamental to the neurophysiology of cellular excitability, these algorithms are not appropriate unless the time between openings and closing is especially small. The exact algorithm [9,24,34] is straightforward to implement and avoids these approximations and pitfalls.