Stochastic Hybrid Systems in Cellular Neuroscience

Abstract

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.

Fig. 1

Schematic illustration of a piecewise deterministic Markov process

Fig. 2

Sketch of a double-well potential of a bistable deterministic system in $\mathbb{R}$

Fig. 3

Opening and closing of ion channels underlying initiation and propagation of an action potential

Fig. 4

Deterministic phase plane dynamics (adapted from [114]). Thick curves show the nullclines: $\dot{v}=0$ as grey and $\dot{w}=0$ as black. Black stream lines represent deterministic trajectories. Green/blue curves represent an action potential trajectory in the limit of slow w. Parameter values are $C_m=20\text{mF}$, $V_{\mathrm{Ca}}=120\text{mV}$, $V_{\mathrm{K}}=-84\text{mV}$, $V_{\mathrm{L}}=-60\text{mV}$, $g_{\mathrm{Ca}}=4.4\text{mS}$, $g_{\mathrm{K}}=8\text{mS}$, $g_{\mathrm{L}}=2.2\text{mS}$, ${\beta }_{\mathrm{Ca}}=0.8{\text{s}}^{-1}$, $v_{\mathrm{Ca},1}=-1.2\text{mV}$, $v_{\mathrm{Ca},2}=18\text{mV}$, $v_{\mathrm{K},1}=2\text{mV}$, $v_{\mathrm{K},2}=30\text{mV}$, and $\varphi =0.04{\text{ms}}^{-1}$

Fig. 5

Sketch of deterministic potential $\Psi (v)$ as a function of voltage v for different values of the applied stimulus current $I_{\mathrm{app}}$. At a critical current $I_{\ast }$, the deterministic system switches from a bistable to a monostable regime, that is, $I^{\ast }$ is the threshold current for action potential generation

Fig. 6

Phase portrait of Hamilton’s equations of motion for the ion channel model with Hamiltonian given by the Perron eigenvalue (3.22). (x and q are taken to be dimensionless.) The zero energy solution representing the maximum likelihood path of escape from $x_{-}$ is shown as the gray curve. (The corresponding path from $x_{+}$ is not shown.) Same parameter values as Fig. 4

Fig. 7

Schematic diagram comparing MFPT calculated using the diffusion approximation with the MFPT of the full system. (Redrawn from [80].) The scales of the axes are based on numerical results for $N=10$. Other parameter values as in Fig. 4

Fig. 8

Electrical coupling via gap junctions. (a) Schematic diagram of gap junction coupling between two cells. (b) Schematic illustration of a Cx43 gap junction channel containing fast (arrow with square) and slow (arrow with hexagon) gates. Voltage gating is mediated by both fast and slow gating mechanisms. Chemical gating is mediated by the slow gating mechanism in both hemichannels

Fig. 9

One-dimensional diffusion in a domain with a randomly switching gate on the right-hand side

Fig. 10

One-dimensional line of cells coupled by gap junctions. At steady-state there is a uniform flux $J_0$ through each cell but a jump discontinuity $\Delta U=-J_0/\mu$ in the concentration across each gap junction, where μ is the permeability of each junction. See text for details

Fig. 11

Pair of cells coupled by a stochastically gated gap junction

Fig. 12

Schematic diagram illustrating volume neurotransmission. Stimulation of an axon terminal contacting a specific synapse on the dendrite of one neuron leads to the release of neurotransmitter within the corresponding synaptic cleft. (A) If neurotransmitter uptake is weak, then it is possible for neurotransmitters to diffuse in the extracellular space and subsequently bind to receptors at other synaptic locations of the same neuron (B) or of another neuron (C)

Fig. 13

Schematic diagram illustrating mRNA granule mobility in dendrites. Under basal conditions, most granules are either stationary (or exhibit localized oscillations), whereas a minority exhibit bidirectional transport. Depolarization by bathing in extracellular KCl solution activates transcription of mRNA at the cell body and converts existing stationary granules into anterograde granules [125]

Fig. 14

Transition diagram of “stop-and-go” model for the slow axonal transport of neurofilaments. See text for definition of different states

Fig. 15

Schematic diagram of an asymmetric tug-of-war model. Two kinesin and two dynein motors transport a cargo in opposite directions along a single polarized microtubule track. Transitions between two possible motor states are shown

Fig. 16

Diagram showing (a) the effective potential well created by a region of tau coating an MT, and (b) a representative trajectory showing random oscillations within the well. (Adapted from [113])

Fig. 17

Schematic illustration of how MAP2 regulation of kinesin motor activities leads to cargo sorting and trafficking in axons. (Redrawn from [68])

Fig. 18

Three-state model of the bidirectional transport of a single motor-cargo complex. The particle switches between an anterograde state ($n=+$) of speed $v_{+}$, a stationary or slowly diffusing state ($n=0$), and a retrograde state ($n=-$) of speed $v_{-}$. The motor-complex can only deliver a vesicle to a presynaptic target in the state $n=0$

Fig. 19

Numerical solutions for steady-state vesicle concentration as a function of axonal distance for different values of $\varphi =k_{-}/{\gamma }_c$ and $J_0=1.5$. (Adapted from [20].) For comparison, the corresponding concentration profile when $J_0=0$ (which is insensitive to ϕ) is shown by the thick line (red line in color online). We have also set $\gamma ={10}^{-2}{\text{s}}^{-1}$, $\hat{J}=1.5$, $k_{+}=0.5{\text{s}}^{-1}$, $k_{-}=1.0\mu {\text{ms}}^{-1}$, $v=\hat{v}=1\mu {\text{ms}}^{-1}$ and $D=\hat{D}=0.1\mu {\text{ms}}^{-2}$

Fig. 20

Different choices of amplitude-phase decomposition. Two possibilities are orthogonal projection with phase ${\theta }^{\prime }(t)$ and isochronal projection with phase $\theta (t)$. In the latter case, the response to perturbations depends on the phase response curve $R(\theta )$, which is normal to the isochron at the point of intersection with the limit cycle

Fig. 21

ML model for subthreshold oscillations. (Adapted from [27].) (a) Bifurcation diagram of the deterministic ML model. As $I_{\mathrm{app}}$ is increased, the system undergoes a supercritical Hopf bifurcation (H) at $I_{\mathrm{app}}^{\ast }=183$ pA, which leads to the generation of stable oscillations. The maximum and minimum values of oscillations are plotted as black (solid) curves. Oscillations disappear via another supercritical Hopf bifurcation. (b), (c) Phase plane diagrams of the deterministic model for (b) $I_{\mathrm{app}}=170\text{pA}$ (below the Hopf bifurcation point) and (c) $I_{\mathrm{app}}=190\text{pA}$ (above the Hopf bifurcation point). The red (dashed) curve is the w-nullcline and the solid (gray) curve represents the v-nullcline. (d), (e) Corresponding voltage time courses. In contrast to Sect. 3.1, we now take ${\alpha }_{\mathrm{K}}={\beta }_{\mathrm{K}}{\mathrm{e}}^{2[v-v_{\mathrm{K},1}]/v_{\mathrm{K},2}}$. Sodium parameters: $g_{\mathrm{Na}}=4.4\text{mS}$, $V_{\mathrm{Na}}=55\text{mV}$, ${\beta }_{\mathrm{Na}}=100{\text{ms}}^{-1}$, $v_{\mathrm{Na},1}=-1.2\text{mV}$, $v_{\mathrm{Na},2}=18\text{mV}$. Leak parameters: $g_{\mathrm{L}}=2.2\text{mS}$, $V_{\mathrm{L}}=-60\text{mV}$. Potassium parameters: $g_{\mathrm{K}}=8\text{mS}$, $V_{\mathrm{K}}=-84\text{mV}$, ${\beta }_{\mathrm{K}}=0.35{\text{ms}}^{-1}$, $v_{\mathrm{K},1}=2\text{mV}$, $v_{\mathrm{K},2}=30\text{mV}$. Also $C_m=1\text{mF}$

Fig. 22

Simulation of the stochastic Morris–Lecar model for subthreshold Na⁺ oscillations with $N=10$ and $\epsilon =0.01$. (Adapted from Ref. [21].) Other parameter values as in Fig. 21. (a) Plot of the approximate phase ${\theta }_t-t{\omega }_0$ in green (with ${\theta }_t$ satisfying equation (6.26) and the exact variational phase (satisfying (6.19)) ${\beta }_t-t{\omega }_0$ in black. On the scale $[-\pi ,\pi ]$ the two phases are in strong agreement. However, zooming in, we can see that the phases slowly drift apart as noise accumulates. The diffusive nature of the drift in both phases can be clearly seen with the typical deviation of the phase from ${\omega }_{0}t$ increasing in time. (b) Stochastic trajectory around limit cycle (dashed curve) in the $v,w$-plane. The stable attractor of the deterministic limit cycle is quite large, which is why the system can tolerate quite substantial stochastic perturbations

Fig. 23

Pair of noninteracting limit cycle oscillators with phases ${\theta }_j(t)$, $j=1,2$, driven by a common switching external input $I(t)$