Efficient stochastic simulation of piecewise-deterministic Markov processes and its application to the Morris-Lecar model of neural dynamics

TL;DR

The paper tackles efficient simulation of piecewise-deterministic Markov processes with state-dependent jump rates by reformulating the next-event problem in terms of a cumulative rate $\Phi$ satisfying $\frac{d\Phi}{dt}=\Lambda$. This transforms inter-event sampling into fixed-interval ODE integration (e.g., solving up to $\Phi=\Delta$ with $\Delta=-\ln r_1$), enabling robust and fast construction of jump times using standard numerical methods. The authors demonstrate the approach on a stochastic Morris–Lecar neuron model, showing controllable accuracy and substantial speed advantages over traditional Gillespie-style event locating and thinning methods, particularly as system size grows. Overall, the method provides a general, efficient framework for simulating PDMPs with time- or state-dependent rates, with proven practicality for neural dynamics and potential extension to related Gillespie variants.

Abstract

Piecewise-deterministic Markov processes combine continuous in time dynamics with jump events, the rates of which generally depend on the continuous variables and thus are not constants. This leads to a problem in a Monte-Carlo simulation of such a system, where, at each step, one must find the time instant of the next event. The latter is determined by an integral equation and usually is rather slow in numerical implementation. We suggest a reformulation of the next event problem as an ordinary differential equation where the independent variable is not the time but the cumulative rate. This reformulation is similar to the Hénon approach to efficiently constructing the Poincaré map in deterministic dynamics. The problem is then reduced to a standard numerical task of solving a system of ordinary differential equations with given initial conditions on a prescribed interval. We illustrate the method with a stochastic Morris-Lecar model of neuron spiking with stochasticity in the opening and closing of voltage-gated ion channels.

Figures (5)

• Figure 1: Time series $V(t),n(t)$ of stochastic simulations of the Morris-Lecar model with different numbers of channels: (a) $N_K=20$; (b) $N_K=40$; (c) $N_K=100$; (d) $N_K=1000$. Right column shows the trajectories on the plane $(V,n)$.
• Figure 2: Maximal errors for different values of $h_0$ in dependence on the number of channels $N_K$. Squares: errors for variable $V$; circles: errors for the times $t$. Red color: $h_0=0.01$; green color: $h_0=0.001$; blue color: $h_0=0.0001$.
• Figure 3: Comparison of trajectories created with the same sequences of random numbers $r_1,r_2$ with the method described in Section \ref{['sec:ssm']} (red curves) and with the approximate method where variations of the rates are neglected on the interval between the events (blue curves). Numbers of channels are given on the panels. Notice the different time ranges of the panels. Integration was performed with the Dormand-Prince method with a fixed step $10^{-3}\cdot N_K$.
• Figure 4: Accuracies in voltages (filled markers) and time differences (open markers) for $N_K=20$ (circles) and $N_K=100$ (squares). The slope of the dashed line is $5.6$, the slope of the dotted line is $4.9$.
• Figure 5: Average CPU times vs accuracies $Err(V)$ (left panel) and $Err(t)$ (right panel) for our method and the methods based on the solution of ODEs Eqs. \ref{['eq:bode']},\ref{['eq:lambdaode']}, for different numbers of iteration steps $m$. The errors are defined as the averaged decimal logarithms, eo, e.g., the value $Err=-8$ corresponds to the absolute error $\approx 10^{-8}$.