A boundary preserving numerical scheme for the Wright-Fisher model
Ioannis S. Stamatiou
TL;DR
This paper introduces a new explicit numerical scheme for non-linear stochastic differential equations that preserves boundary conditions, with proven strong convergence and extensions to multidimensional cases, applicable in biological and physiological modeling.
Contribution
It presents a boundary-preserving explicit numerical scheme for non-linear SDEs, extending previous semi-discrete methods and proving strong convergence for models in population and cellular dynamics.
Findings
The scheme is strongly convergent for the class of SDEs considered.
The method preserves boundary conditions in numerical simulations.
Extension to multidimensional SDEs demonstrated.
Abstract
We are interested in the numerical approximation of non-linear stochastic differential equations (SDEs) with solution in a certain domain. Our goal is to construct explicit numerical schemes that preserve that structure. We generalize the semi-discrete method \emph{Halidias N. and Stamatiou I.S. (2016), On the numerical solution of some non-linear stochastic differential equations using the Semi-Discrete method, Computational Methods in Applied Mathematics,16(1)} and propose a numerical scheme, for which we prove a strong convergence result, to a class of SDEs that appears in population dynamics and ion channel dynamics within cardiac and neuronal cells. We furthermore extend our scheme to a multidimensional case.
| Step | SD | BISS | HYB |
|---|---|---|---|
| 0.008997 | |||
| 0.004216 | |||
| 0.002072 | |||
| 0.001015 | |||
| 0.000509 | |||
| 0.000257 | |||
| 0.000130 | |||
| 0.000067 | |||
| 0.000032 | |||
| 0.000012 |
| Step | SD | BISS | HYB |
|---|---|---|---|
| 0.009030 | |||
| 0.004210 | |||
| 0.002054 | |||
| 0.000999 | |||
| 0.000498 | |||
| 0.000243 | |||
| 0.000120 | |||
| 0.000057 | |||
| 0.000026 | |||
| 0.000011 |
| Step | SD | BISS |
|---|---|---|
| 0.004860 | ||
| 0.002181 | ||
| 0.001169 | ||
| 0.000653 | ||
| 0.000376 | ||
| 0.000231 | ||
| 0.000138 | ||
| 0.000083 | ||
| 0.000051 | ||
| 0.000030 |
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
A boundary preserving numerical scheme for the Wright-Fisher model
I. S. Stamatiou
Abstract.
We are interested in the numerical approximation of non-linear stochastic differential equations (SDEs) with solution in a certain domain. Our goal is to construct explicit numerical schemes that preserve that structure. We generalize the semi-discrete method Halidias N. and Stamatiou I.S. (2016), On the numerical solution of some non-linear stochastic differential equations using the Semi-Discrete method, Computational Methods in Applied Mathematics,16(1) and propose a numerical scheme, for which we prove a strong convergence result, to a class of SDEs that appears in population dynamics and ion channel dynamics within cardiac and neuronal cells. We furthermore extend our scheme to a multidimensional case.
Key words and phrases:
Explicit Numerical Scheme; Semi-Discrete Method; non-linear SDEs; Stochastic Differential Equations; Strong Approximation Error; Boundary Preserving Numerical Algorithm; Wright-Fisher Model.
AMS subject classification 2010: 60H10, 60H35, 65C20, 65C30, 65J15, 65L20, 92D99
1. Introduction
Let and be a complete probability space and let be a one-dimensional Wiener process adapted to the filtration We are interested in the numerical approximation of the following scalar stochastic differential equation (SDE),
[TABLE]
where A boundary classification result, see Appendix A, implies that a.s. when and We therefore aim for a numerical scheme which apart from strongly converging to the true solution of (1.1), produces values in the same domain, i.e. in In other words, we are interested in numerical schemes that have an eternal life time.
**Definition 1.1 **[Eternal Life time of numerical solution] Let and consider a process well defined on the domain with initial condition and such that
[TABLE]
for all A numerical solution has an eternal life time if
[TABLE]
In [1] the main interest is in the domain Moreover, it is clear that the Euler-Maruyama scheme has always a finite life time.
The proposed semi-discrete (SD) iterative scheme for the numerical approximation of (1.1) reads
[TABLE]
where
[TABLE]
and are the increments of the Wiener process and the discretization step is such that (1.2) is well-defined. By construction, the SD scheme (1.2) possesses an eternal life time. To get (1.2) we use an additive semi-discretization of the drift coefficient.. Briefly saying, a part of the SDE is discretized in a certain way such that the resulting SDE to be solved has an analytical solution (see details in Section 2). This is also a special feature of the method, since in the derivation of it, instead of an algebraic equation a new SDE has to be solved. The SD method can also reproduce the Euler scheme. The semi-discrete method was originally proposed in [2].
An attempt in that direction, i.e. in constructing explicit numerical schemes with an eternal life time, has been made in [3] where a class of one-dimensional SDEs with non-negative solutions is treated, which covers cases like that of the Heston -model, a popular model in the field of financial mathematics which is super-linear. The case of sub-linearities is also treated in [4] where the domain is still
The purpose of this paper is to generalize further the method to preserve the structure of the original SDE. In the previous works, the suggested schemes preserve positivity; all the quantities appearing belong to the field of finance and are meant to be non-negative. The application that motivated us now, is used in population dynamics to describe fluctuations in gene frequency of reproducing individuals among finite populations [5] and in a different setting for the description of the random behavior of ion channels within cardiac and neuronal cells (cf. [6], [7], [8] and references therein). We are able in that case to preserve the domain of the original process. In fact, in applications in biology we have to solve systems of SDEs. The extension of the Wright-Fisher model to the multidimensional case has been proposed in [9] and [10]. Here, we will treat a three-state system as in [6, Sec. 6] given by the following system of SDEs
[TABLE]
where is the proportion of alleles or channels in state and is the proportion in state
In Section 2 we provide the setting and the main goal which concerns the mean-square convergence of the proposed structure-preserving SD scheme (1.2) for the approximation of a modification of (1.1) with dynamics described by (see (2.6)). We also discuss the multidimensional case (1.3)-(1.4).
In Section 3 we treat a more general class of SDEs. We further extend the analysis of the semi-discrete method introduced in [3]. We cover the sub-linear diffusion case and show as in [3, Th. 2.1] the strong convergence of the proposed numerical scheme to the true solution.
Section 4 is devoted to numerical experiments. The proofs of all the results are given in the sections to follow, that is in Sections 5 and 6.
2. The setting and the main goal.
Consider the partition with uniform discretization step and the following process
[TABLE]
for with a.s.
[TABLE]
and
[TABLE]
where Process (2.1) has jumps at nodes of order and the solution in each step is given by, see Appendix B,
[TABLE]
which has the pleasant feature that when Process (2.4) is well defined when i.e. when
[TABLE]
Therefore, we assume the following condition for the well-posedness of the SD scheme (2.4).
*Assumption 2.2 *** Let the discretization step be such that (2.5) holds.
*Remark 2.3 *** Note that in general the discretization step satisfying (2.5) is a r.v. depending on The -dependence is inherited through the increments which in turn affect the sequence . Nevertheless under the assumptions on the parameters considered later on the step is not a r.v. but a fixed sufficiently small number.
Now, we consider the process
[TABLE]
which is a martingale with quadratic variation and thus a standard Brownian motion w.r.t. its own filtration, justified by Lévy’s theorem [11, Th. 3.16, p.157] and consequently (2.1) becomes
[TABLE]
Moreover, consider
[TABLE]
The process of (1.1) and the process of (2.8) have the same distribution. Our main goal is to deduce an estimate of the form
[TABLE]
In Theorem 2 below, we deduce that
[TABLE]
By a simple application of the triangle inequality we deduce an analogous result for the unique solution of (1.1), i.e. We present in Appendix C the details. To simplify notation we write as respectively.
**Theorem 2.4 **[Strong convergence] Let Assumption 2 hold. Then, the semi-discrete scheme (2.1) converges strongly in the mean-square sense to the true solution of (1.1), that is
[TABLE]
As already noted in Remark 2, in order to apply Theorem 2 we have to find a sufficiently small step-size such that (2.5) holds, i.e.
[TABLE]
where
[TABLE]
To simplify the conditions on , when necessary, we may adopt the following procedure. We consider a perturbation of order in the initial condition of (2.7), that is
[TABLE]
for with a.s. and
[TABLE]
The following result is a consequence of Theorem 2.
**Proposition 2.5 **[Strong convergence] Let and where is given by (2.10). Then, the semi-discrete scheme (2.11) converges strongly in the mean-square sense to the true solution of (1.1), that is
[TABLE]
Now, we turn to the approximation of the solution of system (1.3)-(1.4). This time we also discretize the diffusion coefficient in a multiplicative way such that the resulting SDEs for each component can be solved analytically. In particular we consider the following processes
[TABLE]
[TABLE]
for with a.s. and
[TABLE]
[TABLE]
where are the deterministic parts of (2.13) and (2.14) respectively and We set for with whenever for a tolerance Finally
Processes (2.13) and (2.14) have jumps at nodes of order and their solution in each step is given respectively by, see Appendix B,
[TABLE]
and
[TABLE]
which has the pleasant feature that when Processes (2.15) and (2.16) are well defined when
Working as before, considering this time the processes
[TABLE]
we conclude to the following result.
**Theorem 2.6 **[Strong convergence] Let the discretization step be such that Then, the semi-discrete scheme (2.13)-(2.14) converges strongly in the mean-square sense to the true solution of (1.3)-(1.4), that is
[TABLE]
where for a -dimensional vector .
We apply the above results in Section 4 and prove them in Section 6.
3. An extension of the semi-discrete method.
Throughout, let and be a complete probability space, meaning that the filtration satisfies the usual conditions, i.e. is right continuous and includes all -null sets. Let be a one-dimensional Wiener process adapted to the filtration Consider the following stochastic differential equation (SDE),
[TABLE]
where the coefficients are measurable functions such that (3.1) has a unique strong solution and is independent of all SDE (3.1) has non-autonomous coefficients, i.e. depend explicitly on
To be more precise, we assume the existence of a predictable stochastic process such that ([12, Def. 2.1]),
[TABLE]
and
[TABLE]
SDEs of the form (3.1) have rarely explicit solutions, thus numerical approximations are necessary for simulations of the paths or for approximation of functionals of the form where We are interested in strong approximations (mean-square) of (3.1), in the case of super- or sub-linear drift and diffusion coefficients and cover cases not included in the previous work [3].
The purpose of this section is to further generalize the semi-discrete (SD) method covering cases of sub-linear diffusion coefficients such as the Cox-Ingersoll-Ross model (CIR) or the Constant Elasticity of Variance model (CEV) (cf. [3, (1.2) and (1.4)]), where also an additive discretization is considered.
*Assumption 3.7 *** Let and be such that where satisfy the following conditions
[TABLE]
for some appropriate (we take in Theorem 3)
[TABLE]
and
[TABLE]
for any such that where the positive parameter the quantity depends on and denotes the maximum of (By the fact that we want the problem (3.1) to be well-posed and by the conditions on and we get that are bounded on bounded intervals.)
Let the equidistant partition and We propose the following semi-discrete numerical scheme
[TABLE]
where we assume that for every (3.2) has a unique strong solution and a.s
[TABLE]
In order to compare with the exact solution which is a continuous time process, we consider the following interpolation process of the semi-discrete approximation, in a compact form,
[TABLE]
where when Process (3.3) has jumps at nodes The first and third variable in denote the discretized part of the original SDE. We observe from (3.3) that in order to solve for , we have to solve an SDE and not an algebraic equation, thus in this context, we cannot reproduce implicit schemes, but we can reproduce the Euler scheme if we choose and
The numerical scheme (3.3) converges to the true solution of SDE (3.1) and this is stated in the following, which is our main result.
**Theorem 3.8 **[Strong convergence] Suppose Assumption 3 holds and (3.2) has a unique strong solution for every where Let also
[TABLE]
for some and Then the semi-discrete numerical scheme (3.3) converges to the true solution of (3.1) in the -sense, that is
[TABLE]
4. Numerics
Here, we make numerical tests to study the strong convergence of the proposed semi-discrete methods (2.1) and (2.11) for the Wright-Fisher model described by the Itô SDE (1.1) with where the parameters and are positive and We take as initial condition the steady state of the deterministic part, that is . This setting has been used for the approximation of ion channels within cardiac and neuronal cells, see [13, Sec. 2.1]; the ion channel occupies one of two positions (open and closed states) with transition rates and respectively and is the total number of ion channels within a cell, see [6, (2.3)]. We consider two set of parameters
- •
SET I:
- •
SET II:
as in [6, Sec. 6 and 7.1], where the Balance Implicit Split Step (BISS) method is suggested [6, (4.8)]
[TABLE]
where is the step-size of the equidistant discretization of the interval , the control function is given by
[TABLE]
and
[TABLE]
The hybrid (HYB) scheme as proposed in [13, (11)] is the result of a splitting method and reads
[TABLE]
It works only for the parameter SET I since we have to assume that
[TABLE]
Finally, the proposed semi-discrete (SD) scheme reads
[TABLE]
for SET I and
[TABLE]
for parameter SET II. The parameters of SET I are chosen in a way that the probability of the Euler-Maruyama (EM) scheme
[TABLE]
leaving the interval is very small, whereas in the case of SET II this probability is high. The paths of the solutions of EM exiting the boundaries [math] and are only a few in the first case and one may reject them. Nevertheless, since such an approach induces bias to the solution obtained by the EM method (much more evident in the second case) we choose only to compare our method with BISS and HYB.
We estimate the endpoint -norm of the difference between the numerical scheme evaluated at step size and the exact solution of (1.1). To do so, we compute batches of simulation paths, where each batch is estimated by and the Monte Carlo estimator of the error is
[TABLE]
and requires Monte Carlo sample paths. The reference solution is calculated using the method at a very fine time grid, We have shown in Theorem 2 and Proposition 2 that the SD numerical schemes converge strongly to the exact solution, so we use the SD method as a reference solution, and the HYB method when applicable. The BISS method converges in the -norm to the true solution [6, Th. 5.1], so we choose not to consider it as a reference solution even though we conjecture that a similar technique may be used to show an -convergence result.
We simulate paths, where the choice of the number of Monte Carlo paths is adequately large, so as not to significantly hinder the mean-square errors. We compute the approximation error (4.5) with -confidence intervals. We present the results for the parameter SET I and II in a - scale in Figures 1 - 3 and Tables 1 - 3 respectively.
Furthermore, we study the strong convergence of the proposed semi-discrete methods (2.13)-(2.14) for the Wright-Fisher model described by the system of Itô SDEs (1.3)-(1.4) with This setting has been used for the approximation of the proportion of alleles or channels in a 3 state system and we consider the set of parameters
- •
SET III:
as in [6, Sec. 6], where the Balance Implicit Split Step (BISS) method is suggested; they use the splitting method, solving first the stochastic system with a balanced implicit scheme
[TABLE]
[TABLE]
and then the deterministic part with the one-step EM scheme; here the control functions are given by
[TABLE]
[TABLE]
[TABLE]
In case then set and if then set The proposed semi-discrete scheme reads, see (2.15)-(2.16)
[TABLE]
[TABLE]
where is the step-size of the equidistant discretization of the interval Note that for the parameter SET III for any The parameters of SET III are chosen in a way that the probability of the Euler-Maruyama (EM) scheme
[TABLE]
leaving the region is very small. The paths of the solutions of EM exiting the boundaries [math] and are extremely few when we take as initial condition the steady state of the deterministic part of the system (which for the parameter SET III is away from the boundaries) and consider a very fine time discretization (). Thus, we take as the exact solution the EM method and reject the paths, if any, outside the region We compare our method with BISS and present the results in Figure 4.
We make the following comments.
- •
The performance of the HYB and SD schemes for parameter Set I is quite similar with SD producing smaller errors. They both perform better than BISS. Nevertheless, HYB is not boundary preserving for parameter SET II. The performance of EM and SD schemes is almost identical for parameter Set III. They both perform better than BISS. Nevertheless, EM is not boundary preserving (we just used it in this experiment for comparative reasons as in [6, Sec. 6].)
- •
The numerical results suggest that the SD schemes converge in the mean-square sense with order for parameter SET I and SET III and at least for parameter SET II.
- •
The proposed SD schemes perform better w.r.t. the computational time required to achieve a desired level of accuracy, since there is no need to calculate a control function.
- •
For the implementation of the SD method, we have to assume that for the parameters of SET I, for the parameters of SET II, and for the parameters of SET III, so the step-size is sufficient. There is not a step-size restriction in the BISS method; nevertheless we propose the SD method since it is fast and more accurate.
5. Proof of Theorem 3.
We split the proof is three steps. First, we prove a general estimate of the error of the SD method for any Then, we show the -convergence of the semi-discrete method and finally the -convergence (3.4). We denote the indicator function of a set by The quantity may vary from line to line and it may depend apart from on other quantities, like time for example, which are all constant, in the sense that we don’t let them grow to infinity.
**Lemma 5.9 **[Error bound for the semi-discrete scheme] Let the assumption of Theorem 3 hold. Let and set the stopping time Then the following estimate holds
[TABLE]
for any where does not depend on implying as
Proof of Lemma 5.
We fix a Let integer such that It holds that
[TABLE]
where we have used Cauchy-Schwarz inequality and Assumption 3 for the function . Taking expectations in the above inequality gives
[TABLE]
where in the third step we have used the Burkholder-Davis-Gundy (BDG) inequality [12, Th. 1.7.3], [11, Th. 3.3.28] on the diffusion term, in the last step Assumption 3 for the function Thus,
[TABLE]
which justifies the notation. Now for we have that
[TABLE]
where we have used Jensen inequality for the concave function ∎
In the next result we estimate the -error using the Yamada-Watanabe approach. We denote the difference
**Proposition 5.10 **[Convergence of the semi-discrete scheme in ] Let the assumptions of Theorem 3 hold. Let and set the stopping time Then we have
[TABLE]
for any where
[TABLE]
and does not depend on It holds that
Proof of Proposition 5.
Let the non-increasing sequence with and We introduce the following sequence of smooth approximations of (method of Yamada and Watanabe, [14])
[TABLE]
where the existence of the continuous function with and support in is justified by The following relations hold for with
[TABLE]
[TABLE]
We have that
[TABLE]
Applying Itô’s formula to the sequence we get
[TABLE]
where in the second step we have used Assumption 3 for the functions the subadditivity property of and the properties of and
[TABLE]
Using the estimate
[TABLE]
we get
[TABLE]
Taking expectations in the above inequality yields
[TABLE]
where we have used Lemma 5 and the fact that .111The function belongs to the space of real-valued measurable -adapted processes such that thus ([12, Th. 1.5.8]) implies . Thus (5.2) becomes
[TABLE]
where in the last step we have used the Gronwall inequality ([15, (7)]) and Taking the supremum over all gives (5.1). ∎
Convergence of the semi-discrete scheme in .
Let the events be defined by and the stopping time for some big enough. We have that
[TABLE]
where is such that the moments of and are bounded by the constant We want to estimate each term of the right hand side of (5.3). It holds that
[TABLE]
for any where the first step comes from the subadditivity of the measure and the second step from Markov inequality. Thus for we get
[TABLE]
We estimate the difference Itô’s formula implies that
[TABLE]
where It holds that
[TABLE]
thus we get that
[TABLE]
Assumption 3 implies that
[TABLE]
Moreover, it holds that
[TABLE]
Taking the supremum over all and then expectation we have
[TABLE]
where we have used Lemma 5 for Using Assumption 3 again we get that
[TABLE]
where we have used the subadditivity property of thus taking expectation we have
[TABLE]
where we have applied again Lemma 5 for and Jensen inequality with We get the following estimate
[TABLE]
where we have used Proposition 5 and
[TABLE]
Plugging the estimates (5.5), (5.6) into (5.4) gives
[TABLE]
where we have applied the Gronwall inequality. Note that, given the bound can be made arbitrarily small by choosing big enough and small enough Relation (5.3) becomes,
[TABLE]
Given any we may first choose such that then choose big enough and small enough such that a concluding as required to verify (3.4).
6. Proof of Theorem 2, Proposition 2 and Theorem 2.
In this Section we prove our main strong convergence result. First, we provide uniform moment bounds for the original SDE and the SD scheme. We remind here that for notational reasons the processes stand for
**Lemma 6.11 **[Moment bounds for original problem and SD approximation] Let Assumption 2 hold. Then
[TABLE]
for any
Proof of Lemma 6.
The result is trivial since we already know that satisfying (1.1) has the property when by Appendix A and regarding the bounds for the SD approximation it is clear, by its form (2.4), that they are valid. ∎
Now, let us rewrite the approximation process
[TABLE]
In the general setting of (3.2) we have
[TABLE]
By the above representation, the form of the discretization becomes apparent. We only discretized the drift coefficient of (1.1) in an additive way. Therefore, by an application of Theorem 3 we have the strong convergence result of Theorem 2
[TABLE]
Now, we briefly sketch the proof of Proposition 2. The process (2.11) is well-defined when or equivalently when and using (2.12) and (2.10). The strong convergence result of Proposition 2 is a consequence of the triangle inequality and the following regularity-type result
[TABLE]
for any where is finite positive.
Theorem 2 is an application of a slight generalization of Theorem 3 including multidimensional noise (see also [16]). Therefore, we omit the proof since one essentially repeats the steps of the proof of Theorem 3. The auxiliary functions in the sense of (3.2) are,
[TABLE]
for the evolution of the first component (1.3), where denotes the discretized part of the SDE, and accordingly for the second component (1.4)
[TABLE]
, see (2.13) and (2.14). By the above representation, the form of the discretization of (1.3) and (1.4) becomes apparent. We discretized the drift coefficient in an additive and multiplicative way and the diffusion coefficient in a multiplicative way.
Acknowledgements
The author would like thank the anonymous referees for their helpful comments.
Appendix A Boundary classification of one-dimensional time-homogeneous SDEs.
Let us now recall some results [11, Sec. 5.5] concerning the boundary behavior of SDEs of the form,
[TABLE]
Let be an interval with and define the exit time from to be
[TABLE]
Let also the coefficients of (A.1) satisfy the following conditions
[TABLE]
[TABLE]
Then for we can define the scale function
[TABLE]
whose behavior at the endpoints of determines the boundary behavior of [11, Prop. 5.5.22]. In particular, we get that the dynamics (1.1) have a boundary behavior which is determined by the scale function
[TABLE]
for any where Let and take We compute
[TABLE]
when and
[TABLE]
when thus by [11, Prop. 5.5.22a] we have that that is
Appendix B Solution process of stochastic integral equations (2.1), (2.13), (2.14).
We will show that the process (2.4) for is the solution of the stochastic integral equation (2.1) for that is
[TABLE]
satisfies
[TABLE]
for with
[TABLE]
Relations (2.6) and (2.3) yield
[TABLE]
where
[TABLE]
The cases for follow with the appropriate modifications.
We can write the increment of the Wiener process as
[TABLE]
and view (B.1) as a function of i.e. with
[TABLE]
and
[TABLE]
Application of Itô’s formula implies
[TABLE]
For the derivation of (2.13) and (2.14) we now write the multidimensional Wiener process as
[TABLE]
where is the zero matrix and the identity matrix and apply appropriately the multidimensional Itô formula.
Appendix C Uniform moment estimate for .
In Theorem 2 we actually proved that In order to finish the proof we have to find a uniform moment bound for In particular applying the triangle inequality
[TABLE]
thus it suffices to show
[TABLE]
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] H. Schurz. Numerical regularization for SD Es: Construction of nonnegative solutions. Dynamic Systems and Applications , 5(3):323–351, 1996.
- 2[2] N. Halidias. Semi-discrete approximations for stochastic differential equations and applications. International Journal of Computer Mathematics , 89(6):780–794, 2012.
- 3[3] N. Halidias and I.S. Stamatiou. On the Numerical Solution of Some Non-Linear Stochastic Differential Equations Using the Semi-Discrete Method. Computational Methods in Applied Mathematics , 16(1):105–132, 2016.
- 4[4] N. Halidias and I.S. Stamatiou. Approximating Explicitly the Mean-Reverting CEV Process. Journal of Probability and Statistics , Article ID 513137, 20 pages, 2015.
- 5[5] W. J. Ewens. Mathematical Population Genetics 1: Theoretical Introduction , volume 27. Springer Science & Business Media, 2012.
- 6[6] C.E. Dangerfield, D. Kay, S. Mac Namara, and K Burrage. A boundary preserving numerical algorithm for the wright-fisher model with mutation. BIT Numerical Mathematics , 52(2):283–304, 2012.
- 7[7] J. H. Goldwyn, N. S. Imennov, M. Famulare, and E. Shea-Brown. Stochastic differential equation models for ion channel noise in hodgkin-huxley neurons. Physical Review E , 83(4):041908, 2011.
- 8[8] C. E. Dangerfield, D. Kay, and K. Burrage. Modeling ion channel dynamics through reflected stochastic differential equations. Physical Review E , 85(5):051907, 2012.