Stochastic representations of ion channel kinetics and exact stochastic simulation of neuronal dynamics

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

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.

Competition of two time scales determines the performance of a voltage-gated potassium channel

The dynamics of a voltage-gated potassium channel in real environments represents a crucial bridge between its molecular structure and functions. However, it is still missing due to the mathematical difficulty that arises from the high dimensionality and nonlinear interregulation. Here we present a method for solving the stationary distribution of a hybrid process that contains two negatively interregulating kinetics: channel gating and voltage decay. The results can be summarized as follows: first, the voltage distribution is determined by the competition of their time scales; second, the fluctuation structures in parameter space illustrate that, to perform the voltage-controlling task, the channel gating is elastic while the membrane produces the stabilizing function; third, the power dissipated by the capacitive currents and the internal battery current are calculated and explained. Based on these findings, we examine the manner in which macroscopic functions of potassium channels are manifested. Our methodology provides an accurate characterization of hybrid processes that are pervasive in the life sciences.

Revisiting kinetic Monte Carlo algorithms for time-dependent processes: From open-loop control to feedback control

Simulating stochastic systems with feedback control is challenging due to the complex interplay between the system's dynamics and the feedback-dependent control protocols. We present a single-step-trajectory probability analysis to time-dependent stochastic systems. Based on this analysis, we revisit several time-dependent kinetic Monte Carlo (KMC) algorithms designed for systems under open-loop-control protocols. Our analysis provides a unified alternative proof to these algorithms, summarized into a pedagogical tutorial. Moreover, with the trajectory probability analysis, we present a novel feedback-controlled KMC algorithm that accurately captures the dynamics systems controlled by an external signal based on the measurements of the system's state. Our method correctly captures the system dynamics and avoids the artificial Zeno effect that arises from incorrectly applying the direct Gillespie algorithm to feedback-controlled systems. This work provides a unified perspective on existing open-loop-control KMC algorithms and also offers a powerful and accurate tool for simulating stochastic systems with feedback control.

Heteroclinic cycling and extinction in May-Leonard models with demographic stochasticity

May and Leonard (SIAM J Appl Math 29:243-253, 1975) introduced a three-species Lotka-Volterra type population model that exhibits heteroclinic cycling. Rather than producing a periodic limit cycle, the trajectory takes longer and longer to complete each "cycle", passing closer and closer to unstable fixed points in which one population dominates and the others approach zero. Aperiodic heteroclinic dynamics have subsequently been studied in ecological systems (side-blotched lizards; colicinogenic Escherichia coli), in the immune system, in neural information processing models ("winnerless competition"), and in models of neural central pattern generators. Yet as May and Leonard observed "Biologically, the behavior (produced by the model) is nonsense. Once it is conceded that the variables represent animals, and therefore cannot fall below unity, it is clear that the system will, after a few cycles, converge on some single population, extinguishing the other two." Here, we explore different ways of introducing discrete stochastic dynamics based on May and Leonard's ODE model, with application to ecological population dynamics, and to a neuromotor central pattern generator system. We study examples of several quantitatively distinct asymptotic behaviors, including total extinction of all species, extinction to a single species, and persistent cyclic dominance with finite mean cycle length.

DeepGANnel: Synthesis of fully annotated single molecule patch-clamp data using generative adversarial networks

Development of automated analysis tools for "single ion channel" recording is hampered by the lack of available training data. For machine learning based tools, very large training sets are necessary with sample-by-sample point labelled data (e.g., 1 sample point every 100microsecond). In an experimental context, such data are labelled with human supervision, and whilst this is feasible for simple experimental analysis, it is infeasible to generate the enormous datasets that would be necessary for a big data approach using hand crafting. In this work we aimed to develop methods to generate simulated ion channel data that is free from assumptions and prior knowledge of noise and underlying hidden Markov models. We successfully leverage generative adversarial networks (GANs) to build an end-to-end pipeline for generating an unlimited amount of labelled training data from a small, annotated ion channel "seed" record, and this needs no prior knowledge of theoretical dynamical ion channel properties. Our method utilises 2D CNNs to maintain the synchronised temporal relationship between the raw and idealised record. We demonstrate the applicability of the method with 5 different data sources and show authenticity with t-SNE and UMAP projection comparisons between real and synthetic data. The model would be easily extendable to other time series data requiring parallel labelling, such as labelled ECG signals or raw nanopore sequencing data.

Quantitative comparison of the mean-return-time phase and the stochastic asymptotic phase for noisy oscillators

Seminal work by A. Winfree and J. Guckenheimer showed that a deterministic phase variable can be defined either in terms of Poincaré sections or in terms of the asymptotic (long-time) behaviour of trajectories approaching a stable limit cycle. However, this equivalence between the deterministic notions of phase is broken in the presence of noise. Different notions of phase reduction for a stochastic oscillator can be defined either in terms of mean-return-time sections or as the argument of the slowest decaying complex eigenfunction of the Kolmogorov backwards operator. Although both notions of phase enjoy a solid theoretical foundation, their relationship remains unexplored. Here, we quantitatively compare both notions of stochastic phase. We derive an expression relating both notions of phase and use it to discuss differences (and similarities) between both definitions of stochastic phase for (i) a spiral sink motivated by stochastic models for electroencephalograms, (ii) noisy limit-cycle systems-neuroscience models, and (iii) a stochastic heteroclinic oscillator inspired by a simple motor-control system.

Mechanisms and Physiological Implications of Cooperative Gating of Ion Channels Clusters

Ion channels play a central role in the regulation of nearly every cellular process. Dating back to the classic 1952 Hodgkin-Huxley model of the generation of the action potential, ion channels have always been thought of as independent agents. A myriad of recent experimental findings exploiting advances in electrophysiology, structural biology, and imaging techniques, however, have posed a serious challenge to this long-held axiom as several classes of ion channels appear to open and close in a coordinated, cooperative manner. Ion channel cooperativity ranges from variable-sized oligomeric cooperative gating in voltage-gated, dihydropyridine-sensitive Ca_v1.2 and Ca_v1.3 channels to obligatory dimeric assembly and gating of voltage-gated Na_v1.5 channels. Potassium channels, transient receptor potential channels, hyperpolarization cyclic nucleotide-activated channels, ryanodine receptors (RyRs), and inositol trisphosphate receptors (IP₃Rs) have also been shown to gate cooperatively. The implications of cooperative gating of these ion channels range from fine tuning excitation-contraction coupling in muscle cells to regulating cardiac function and vascular tone, to modulation of action potential and conduction velocity in neurons and cardiac cells, and to control of pace-making activity in the heart. In this review, we discuss the mechanisms leading to cooperative gating of ion channels, their physiological consequences and how alterations in cooperative gating of ion channels may induce a range of clinically significant pathologies.

Resolving molecular contributions of ion channel noise to interspike interval variability through stochastic shielding

Molecular fluctuations can lead to macroscopically observable effects. The random gating of ion channels in the membrane of a nerve cell provides an important example. The contributions of independent noise sources to the variability of action potential timing have not previously been studied at the level of molecular transitions within a conductance-based model ion-state graph. Here we study a stochastic Langevin model for the Hodgkin-Huxley (HH) system based on a detailed representation of the underlying channel state Markov process, the "[Formula: see text]D model" introduced in (Pu and Thomas in Neural Computation 32(10):1775-1835, 2020). We show how to resolve the individual contributions that each transition in the ion channel graph makes to the variance of the interspike interval (ISI). We extend the mean return time (MRT) phase reduction developed in (Cao et al. in SIAM J Appl Math 80(1):422-447, 2020) to the second moment of the return time from an MRT isochron to itself. Because fixed-voltage spike detection triggers do not correspond to MRT isochrons, the inter-phase interval (IPI) variance only approximates the ISI variance. We find the IPI variance and ISI variance agree to within a few percent when both can be computed. Moreover, we prove rigorously, and show numerically, that our expression for the IPI variance is accurate in the small noise (large system size) regime; our theory is exact in the limit of small noise. By selectively including the noise associated with only those few transitions responsible for most of the ISI variance, our analysis extends the stochastic shielding (SS) paradigm (Schmandt and Galán in Phys Rev Lett 109(11):118101, 2012) from the stationary voltage clamp case to the current clamp case. We show numerically that the SS approximation has a high degree of accuracy even for larger, physiologically relevant noise levels. Finally, we demonstrate that the ISI variance is not an unambiguously defined quantity, but depends on the choice of voltage level set as the spike detection threshold. We find a small but significant increase in ISI variance, the higher the spike detection voltage, both for simulated stochastic HH data and for voltage traces recorded in in vitro experiments. In contrast, the IPI variance is invariant with respect to the choice of isochron used as a trigger for counting "spikes."

Dynamical consequences of sensory feedback in a half-center oscillator coupled to a simple motor system

We investigate a simple model for motor pattern generation that combines central pattern generator (CPG) dynamics with a sensory feedback (FB) mechanism. Our CPG comprises a half-center oscillator with conductance-based Morris-Lecar model neurons. Output from the CPG drives a push-pull motor system with biomechanics based on experimental data. A sensory feedback conductance from the muscles allows modulation of the CPG activity. We consider parameters under which the isolated CPG system has either "escape" or "release" dynamics, and we study both inhibitory and excitatory feedback conductances. We find that increasing the FB conductance relative to the CPG conductance makes the system more robust against external perturbations, but more susceptible to internal noise. Conversely, increasing the CPG conductance relative to the FB conductance has the opposite effects. We find that the "closed-loop" system, with sensory feedback in place, exhibits a richer repertoire of behaviors than the "open-loop" system, with motion determined entirely by the CPG dynamics. Moreover, we find that purely feedback-driven motor patterns, analogous to a chain reflex, occur only in the inhibition-mediated system. Finally, for pattern generation systems with inhibition-mediated sensory feedback, we find that the distinction between escape- and release-mediated CPG mechanisms is diminished in the presence of internal noise. Our observations support an anti-reductionist view of neuromotor physiology: Understanding mechanisms of robust motor control requires studying not only the central pattern generator circuit in isolation, but the intact closed-loop system as a whole.

Toward Closed-Loop Electrical Stimulation of Neuronal Systems: A Review

Biological neuronal cells communicate using neurochemistry and electrical signals. The same phenomena also allow us to probe and manipulate neuronal systems and communicate with them. Neuronal system malfunctions cause a multitude of symptoms and functional deficiencies that can be assessed and sometimes alleviated by electrical stimulation. Our working hypothesis is that real-time closed-loop full-duplex measurement and stimulation paradigms can provide more in-depth insight into neuronal networks and enhance our capability to control diseases of the nervous system. In this study, we review extracellular electrical stimulation methods used in in vivo, in vitro, and in silico neuroscience research and in the clinic (excluding methods mainly aimed at neuronal growth and other similar effects) and highlight the potential of closed-loop measurement and stimulation systems. A multitude of electrical stimulation and measurement-based methods are widely used in research and the clinic. Closed-loop methods have been proposed, and some are used in the clinic. However, closed-loop systems utilizing more complex measurement analysis and adaptive stimulation systems, such as artificial intelligence systems connected to biological neuronal systems, do not yet exist. Our review promotes the research and development of intelligent paradigms aimed at meaningful communications between neuronal and information and communications technology systems, "dialogical paradigms," which have the potential to take neuroscience and clinical methods to a new level.

Fast and Accurate Langevin Simulations of Stochastic Hodgkin-Huxley Dynamics

Fox and Lu introduced a Langevin framework for discrete-time stochastic models of randomly gated ion channels such as the Hodgkin-Huxley (HH) system. They derived a Fokker-Planck equation with state-dependent diffusion tensor D and suggested a Langevin formulation with noise coefficient matrix S such that SS⊤=D. Subsequently, several authors introduced a variety of Langevin equations for the HH system. In this article, we present a natural 14-dimensional dynamics for the HH system in which each directed edge in the ion channel state transition graph acts as an independent noise source, leading to a 14 × 28 noise coefficient matrix S. We show that (1) the corresponding 14D system of ordinary differential equations is consistent with the classical 4D representation of the HH system; (2) the 14D representation leads to a noise coefficient matrix S that can be obtained cheaply on each time step, without requiring a matrix decomposition; (3) sample trajectories of the 14D representation are pathwise equivalent to trajectories of Fox and Lu's system, as well as trajectories of several existing Langevin models; (4) our 14D representation (and those equivalent to it) gives the most accurate interspike interval distribution, not only with respect to moments but under both the L1 and L∞ metric-space norms; and (5) the 14D representation gives an approximation to exact Markov chain simulations that are as fast and as efficient as all equivalent models. Our approach goes beyond existing models, in that it supports a stochastic shielding decomposition that dramatically simplifies S with minimal loss of accuracy under both voltage- and current-clamp conditions.

Ion channel noise shapes the electrical activity of endocrine cells

Endocrine cells in the pituitary gland typically display either spiking or bursting electrical activity, which is related to the level of hormone secretion. Recent work, which combines mathematical modelling with dynamic clamp experiments, suggests the difference is due to the presence or absence of a few large-conductance potassium channels. Since endocrine cells only contain a handful of these channels, it is likely that stochastic effects play an important role in the pattern of electrical activity. Here, for the first time, we explicitly determine the effect of such noise by studying a mathematical model that includes the realistic noisy opening and closing of ion channels. This allows us to investigate how noise affects the electrical activity, examine the origin of spiking and bursting, and determine which channel types are responsible for the greatest noise. Further, for the first time, we address the role of cell size in endocrine cell electrical activity, finding that larger cells typically display more bursting, while the smallest cells almost always only exhibit spiking behaviour.

Clusters of cooperative ion channels enable a membrane-potential-based mechanism for short-term memory

Across biological systems, cooperativity between proteins enables fast actions, supra-linear responses, and long-lasting molecular switches. In the nervous system, however, the function of cooperative interactions between voltage-dependent ionic channels remains largely unknown. Based on mathematical modeling, we here demonstrate that clusters of strongly cooperative ion channels can plausibly form bistable conductances. Consequently, clusters are permanently switched on by neuronal spiking, switched off by strong hyperpolarization, and remain in their state for seconds after stimulation. The resulting short-term memory of the membrane potential allows to generate persistent firing when clusters of cooperative channels are present together with non-cooperative spike-generating conductances. Dynamic clamp experiments in rodent cortical neurons confirm that channel cooperativity can robustly induce graded persistent activity - a single-cell based, multistable mnemonic firing mode experimentally observed in several brain regions. We therefore propose that ion channel cooperativity constitutes an efficient cell-intrinsic implementation for short-term memories at the voltage level.

Fluctuation Scaling of Neuronal Firing and Bursting in Spontaneously Active Brain Circuits

We employed high-density microelectrode arrays to investigate spontaneous firing patterns of neurons in brain circuits of the primary somatosensory cortex (S1) in mice. We recorded from over 150 neurons for 10min in each of eight different experiments, identified their location in S1, sorted their action potentials (spikes), and computed their power spectra and inter-spike interval (ISI) statistics. Of all persistently active neurons, 92% fired with a single dominant frequency - regularly firing neurons (RNs) - from 1 to 8Hz while 8% fired in burst with two dominant frequencies - bursting neurons (BNs) - corresponding to the inter-burst (2-6Hz) and intra-burst intervals (20-160Hz). RNs were predominantly located in layers 2/3 and 5/6 while BNs localized to layers 4 and 5. Across neurons, the standard deviation of ISI was a power law of its mean, a property known as fluctuation scaling, with a power law exponent of 1 for RNs and 1.25 for BNs. The power law implies that firing and bursting patterns are scale invariant: the firing pattern of a given RN or BN resembles that of another RN or BN, respectively, after a time contraction or dilation. An explanation for this scale invariance is discussed in the context of previous computational studies as well as its potential role in information processing.

Special Issue from the 2017 International Conference on Mathematical Neuroscience

The ongoing acquisition of large and multifaceted data sets in neuroscience requires new mathematical tools for quantitatively grounding these experimental findings. Since 2015, the International Conference on Mathematical Neuroscience (ICMNS) has provided a forum for researchers to discuss current mathematical innovations emerging in neuroscience. This special issue assembles current research and tutorials that were presented at the 2017 ICMNS held in Boulder, Colorado from May 30 to June 2. Topics discussed at the meeting include correlation analysis of network activity, information theory for plastic synapses, combinatorics for attractor neural networks, and novel data assimilation methods for neuroscience-all of which are represented in this special issue.

Stochastic Mechanochemical Description of a Bioinspired Polymerization Process

We present a theoretical investigation of a polymerization process catalyzed by an enzyme. A structural model of enzyme, sliding along the polymer chain as a Brownian particle, is proposed, and a stochastic approach is employed to describe the kinetics of the whole process. The key point of this work is the coupling mechanics/chemistry obtained by assuming that (1) some rates of chemical reaction depend on the position of the enzyme with respect to the polymer chain and (2) the potential energy and the friction coefficient in the Langevin equation depend on the chemical state of the polymerizing complex. We describe an algorithm for computing our stochastic model and a methodology to solve the Langevin equation numerically. We predict in particular: (1) the sudden arrest of the polymerization, (2) the decrease in the relative polydispersity with the increase in the length of the polymer chain, (3) the occurrence of four regimes, (4) the manifestation of the coupling mechanics/chemistry for one regime and (5) the possibility to evaluate the mechanical variables through classical chemical analysis. Although essentially devoted to the elongation phase, this work also briefly addresses the problem of phase termination and we propose a new device aimed at reducing the polydispersity of technical origin in actual polymerization processes.

A discrete in continuous mathematical model of cardiac progenitor cells formation and growth as spheroid clusters (Cardiospheres)

We propose a discrete in continuous mathematical model describing the in vitro growth process of biophsy-derived mammalian cardiac progenitor cells growing as clusters in the form of spheres (Cardiospheres). The approach is hybrid: discrete at cellular scale and continuous at molecular level. In the present model, cells are subject to the self-organizing collective dynamics mechanism and, additionally, they can proliferate and differentiate, also depending on stochastic processes. The two latter processes are triggered and regulated by chemical signals present in the environment. Numerical simulations show the structure and the development of the clustered progenitors and are in a good agreement with the results obtained from in vitro experiments.

Stochastic Hybrid Systems in Cellular Neuroscience

We review recent work on the theory and applications of stochastic hybrid systems in cellular neuroscience. A stochastic hybrid system or piecewise deterministic Markov process involves the coupling between a piecewise deterministic differential equation and a time-homogeneous Markov chain on some discrete space. The latter typically represents some random switching process. We begin by summarizing the basic theory of stochastic hybrid systems, including various approximation schemes in the fast switching (weak noise) limit. In subsequent sections, we consider various applications of stochastic hybrid systems, including stochastic ion channels and membrane voltage fluctuations, stochastic gap junctions and diffusion in randomly switching environments, and intracellular transport in axons and dendrites. Finally, we describe recent work on phase reduction methods for stochastic hybrid limit cycle oscillators.

Stochastic shielding and edge importance for Markov chains with timescale separation

Nerve cells produce electrical impulses (“spikes”) through the coordinated opening and closing of ion channels. Markov processes with voltage-dependent transition rates capture the stochasticity of spike generation at the cost of complex, time-consuming simulations. Schmandt and Galán introduced a novel method, based on the stochastic shielding approximation, as a fast, accurate method for generating approximate sample paths with excellent first and second moment agreement to exact stochastic simulations. We previously analyzed the mathematical basis for the method’s remarkable accuracy, and showed that for models with a Gaussian noise approximation, the stationary variance of the occupancy at each vertex in the ion channel state graph could be written as a sum of distinct contributions from each edge in the graph. We extend this analysis to arbitrary discrete population models with first-order kinetics. The resulting decomposition allows us to rank the “importance” of each edge’s contribution to the variance of the current under stationary conditions. In most cases, transitions between open (conducting) and closed (non-conducting) states make the greatest contributions to the variance, but there are exceptions. In a 5-state model of the nicotinic acetylcholine receptor, at low agonist concentration, a pair of “hidden” transitions (between two closed states) makes a greater contribution to the variance than any of the open-closed transitions. We exhaustively investigate this “edge importance reversal” phenomenon in simplified 3-state models, and obtain an exact formula for the contribution of each edge to the variance of the open state. Two conditions contribute to reversals: the opening rate should be faster than all other rates in the system, and the closed state leading to the opening rate should be sparsely occupied. When edge importance reversal occurs, current fluctuations are dominated by a slow noise component arising from the hidden transitions.

Discrete state, continuous time Markov processes occur throughout cell biology, neuroscience, and ecology, representing the random dynamics of processes transitioning among multiple locations or states. Complexity reduction for such models aims to capture the essential dynamics and stochastic properties via a simpler representation, with minimal loss of accuracy. Classical approaches, such as aggregation of nodes and elimination of fast variables, lead to reduced models that are no longer Markovian. Stochastic shielding provides an alternative approach by simplifying the description of the noise driving the process, while preserving the Markov property, by removing from the model those fluctuations that are not directly observable. We previously applied the stochastic shielding approximation to several Markov processes arising in neuroscience and processes on random graphs. Here we explore the range of validity of stochastic shielding for processes with nonuniform stationary probabilities and multiple timescales, including ion channels with “bursty” dynamics. We show that stochastic shielding is robust to the introduction of timescale separation, for a class of simple networks, but it can break down for more complex systems with three distinct timescales. We also show that our related edge importance measure remains valid for arbitrary networks regardless of multiple timescales.

A variational method for analyzing limit cycle oscillations in stochastic hybrid systems

Many systems in biology can be modeled through ordinary differential equations, which are piece-wise continuous, and switch between different states according to a Markov jump process known as a stochastic hybrid system or piecewise deterministic Markov process (PDMP). In the fast switching limit, the dynamics converges to a deterministic ODE. In this paper, we develop a phase reduction method for stochastic hybrid systems that support a stable limit cycle in the deterministic limit. A classic example is the Morris-Lecar model of a neuron, where the switching Markov process is the number of open ion channels and the continuous process is the membrane voltage. We outline a variational principle for the phase reduction, yielding an exact analytic expression for the resulting phase dynamics. We demonstrate that this decomposition is accurate over timescales that are exponential in the switching rate ϵ^-1. That is, we show that for a constant C, the probability that the expected time to leave an O(a) neighborhood of the limit cycle is less than T scales as T exp (-Ca/ϵ).

On two diffusion neuronal models with multiplicative noise: The mean first-passage time properties

Two diffusion processes with multiplicative noise, able to model the changes in the neuronal membrane depolarization between two consecutive spikes of a single neuron, are considered and compared. The processes have the same deterministic part but different stochastic components. The differences in the state-dependent variabilities, their asymptotic distributions, and the properties of the first-passage time across a constant threshold are investigated. Closed form expressions for the mean of the first-passage time of both processes are derived and applied to determine the role played by the parameters involved in the model. It is shown that for some values of the input parameters, the higher variability, given by the second moment, does not imply shorter mean first-passage time. The reason for that can be found in the complete shape of the stationary distribution of the two processes. Applications outside neuroscience are also mentioned.

Feynman-Kac formula for stochastic hybrid systems

We derive a Feynman-Kac formula for functionals of a stochastic hybrid system evolving according to a piecewise deterministic Markov process. We first derive a stochastic Liouville equation for the moment generator of the stochastic functional, given a particular realization of the underlying discrete Markov process; the latter generates transitions between different dynamical equations for the continuous process. We then analyze the stochastic Liouville equation using methods recently developed for diffusion processes in randomly switching environments. In particular, we obtain dynamical equations for the moment generating function, averaged with respect to realizations of the discrete Markov process. The resulting Feynman-Kac formula takes the form of a differential Chapman-Kolmogorov equation. We illustrate the theory by calculating the occupation time for a one-dimensional velocity jump process on the infinite or semi-infinite real line. Finally, we present an alternative derivation of the Feynman-Kac formula based on a recent path-integral formulation of stochastic hybrid systems.

Numerical simulations of piecewise deterministic Markov processes with an application to the stochastic Hodgkin-Huxley model

The stochastic Hodgkin-Huxley model is one of the best-known examples of piecewise deterministic Markov processes (PDMPs), in which the electrical potential across a cell membrane, V(t), is coupled with a mesoscopic Markov jump process representing the stochastic opening and closing of ion channels embedded in the membrane. The rates of the channel kinetics, in turn, are voltage-dependent. Due to this interdependence, an accurate and efficient sampling of the time evolution of the hybrid stochastic systems has been challenging. The current exact simulation methods require solving a voltage-dependent hitting time problem for multiple path-dependent intensity functions with random thresholds. This paper proposes a simulation algorithm that approximates an alternative representation of the exact solution by fitting the log-survival function of the inter-jump dwell time, H(t), with a piecewise linear one. The latter uses interpolation points that are chosen according to the time evolution of the H(t), as the numerical solution to the coupled ordinary differential equations of V(t) and H(t). This computational method can be applied to all PDMPs. Pathwise convergence of the approximated sample trajectories to the exact solution is proven, and error estimates are provided. Comparison with a previous algorithm that is based on piecewise constant approximation is also presented.

Channel-noise-induced critical slowing in the subthreshold Hodgkin-Huxley neuron

The dynamics of a spiking neuron approaching threshold is investigated in the framework of Markov-chain models describing the random state-transitions of the underlying ion-channel proteins. We characterize subthreshold channel-noise-induced transmembrane potential fluctuations in both type-I (integrator) and type-II (resonator) parametrizations of the classic conductance-based Hodgkin-Huxley equations. As each neuron approaches spiking threshold from below, numerical simulations of stochastic trajectories demonstrate pronounced growth in amplitude simultaneous with decay in frequency of membrane voltage fluctuations induced by ion-channel state transitions. To explore this progression of fluctuation statistics, we approximate the exact Markov treatment with a 12-variable channel-based stochastic differential equation (SDE) and its Ornstein-Uhlenbeck (OU) linearization and show excellent agreement between Markov and SDE numerical simulations. Predictions of the OU theory with respect to membrane potential fluctuation variance, autocorrelation, correlation time, and spectral density are also in agreement and illustrate the close connection between the eigenvalue structure of the associated deterministic bifurcations and the observed behavior of the noisy Markov traces on close approach to threshold for both integrator and resonator point-neuron varieties.

Hybrid pathwise sensitivity methods for discrete stochastic models of chemical reaction systems

Stochastic models are often used to help understand the behavior of intracellular biochemical processes. The most common such models are continuous time Markov chains (CTMCs). Parametric sensitivities, which are derivatives of expectations of model output quantities with respect to model parameters, are useful in this setting for a variety of applications. In this paper, we introduce a class of hybrid pathwise differentiation methods for the numerical estimation of parametric sensitivities. The new hybrid methods combine elements from the three main classes of procedures for sensitivity estimation and have a number of desirable qualities. First, the new methods are unbiased for a broad class of problems. Second, the methods are applicable to nearly any physically relevant biochemical CTMC model. Third, and as we demonstrate on several numerical examples, the new methods are quite efficient, particularly if one wishes to estimate the full gradient of parametric sensitivities. The methods are rather intuitive and utilize the multilevel Monte Carlo philosophy of splitting an expectation into separate parts and handling each in an efficient manner.

Asymptotic phase for stochastic oscillators

Oscillations and noise are ubiquitous in physical and biological systems. When oscillations arise from a deterministic limit cycle, entrainment and synchronization may be analyzed in terms of the asymptotic phase function. In the presence of noise, the asymptotic phase is no longer well defined. We introduce a new definition of asymptotic phase in terms of the slowest decaying modes of the Kolmogorov backward operator. Our stochastic asymptotic phase is well defined for noisy oscillators, even when the oscillations are noise dependent. It reduces to the classical asymptotic phase in the limit of vanishing noise. The phase can be obtained either by solving an eigenvalue problem, or by empirical observation of an oscillating density's approach to its steady state.