# Large Deviations of Piecewise-Deterministic-Markov-Processes with
Application to Stochastic Calcium Waves

Gaetan Barbet, James MacLaurin and Moshe Silverstein

###### Abstract

We prove a Large Deviation Principle for Piecewise Deterministic Markov
Processes (PDMPs). This is an asymptotic estimate for the probability of a
trajectory in the large size limit. Explicit Euler-Lagrange equations are
determined for computing optimal first-hitting-time trajectories. The results
are applied to a model of stochastic calcium dynamics. It is widely
conjectured that the mechanism of calcium puff generation is a multiscale
process: with microscopic stochastic fluctuations in the opening and closing
of individual channels generating cell-wide waves via the diffusion of calcium
and other signaling molecules. We model this system as a PDMP, with $N\gg 1$
stochastic calcium channels that are coupled via the ambient calcium
concentration. We employ the Large Deviations theory to estimate the
probability of cell-wide calcium waves being produced through microscopic
stochasticity.

This paper determines a Large Deviation Principle for Piecewise-Deterministic
Markov Processes (PDMPs), in the large population size limit. Piecewise-
Deterministic Markov Processes (also known as stochastic hybrid systems) are
used to model systems with multiple timescales [25, 4]. They consist of a
Markov Chain flipping between discrete states, coupled to an ordinary
differential equation that evolves deterministically in-between the flips (and
whose evolution typically depends on the state of the Markov Chain). The are
also referred to as slow-fast systems. They enjoy numerous applications,
particularly in biology [47, 12, 8, 10]. They have been used (for instance) to
model excitable membranes in neuroscience [46, 39], population dynamics in
ecology [27], run-and-tumble dynamics of bacteria [5] and stochastic models of
calcium signals [33, 44, 45].

Often, one wishes to estimate the probability of rare-events in these systems
[35]. For example, one may be interested in the probability of extinction of a
species in population dynamics or a chemical reaction networks [11], or the
spontaneous and stochastic production of waves / pulses [31]. Another well-
known application - which is thoroughly explored in this paper - is stochastic
calcium dynamics (this will be discussed in more detail below). One of the
best tools for estimating the probability of rare events is the theory of
Large Deviations [13]. Large Deviations analysis has been performed on many
stochastic differential equations and jump-Markov Processes, including [26,
7]. However it is comparatively less-developed for piecewise-deterministic
Markov Processes.

To the author’s knowledge, the first scholar to study the Large Deviations of
PDMPs was Kifer [34]. This was followed by the work of Faggionato, Crivellari
and Gabrielli [22], who determined a Large Deviation Principle for piecewise-
deterministic Markov Processes. They require that (conditionally on the value
of the ‘slow variable’), the Markov Process is ergodic. Kumar and Popovic also
determined a Large Deviation Principle for stochastic hybrid systems [37]. To
do this, they used an ergodicity assumption to prove that the logarithmic
moment generating function is sufficiently regular. Bressloff and Faugeras [9]
demonstrated that the Large Deviations rate function can be written to have a
Hamiltonian structure. [36] proved a general Large Deviation Principle for
slow-fast systems, by proving convergence of the Hamilton-Jacobi-Bellman
equation.

On a related note, there has been much recent interest in the Large Deviations
of Chemical Reaction Networks [43]. These are high-dimensional jump Markovian
Processes, similar to the subject of this paper, except without the ‘slow’
dynamics. Recent works include those of Dupuis, Ramanan and Wu [20], Pardoux
and Samegni-Kepgnou [42], Agazzi, Eckmann and Dembo [2], and Patterson and
Renger [43]. Patterson and Renger [43] and Agazzi, Patterson, Renger and [1]
prove the Large Deviations Principle in a general setting by also studying the
convergence of the reaction fluxes (like in this paper). Budhiraja and
Friedlander [14] determine the Large Deviations of jump-Markov processes.

One of the novelties of the proof in this paper is its simplicity. We map the
PDMP system to a set of homogeneous Poisson Processes via a time-rescaling.
The Large Deviation Principle is then an almost immediate result of the
Inverse Contraction Principle [18]. In addition, we explicitly determine the
structure of the Euler-Lagrange equations (needed for computing optimal first-
hitting times). This becomes an optimization problem over the time interval
$[0,T]$, subject to a constraint (given by the slow processes). We then solve
these equations to compute the probability of noise-induced calcium waves.

###  0.1 Overview of Stochastic Models of Calcium Signaling

Classical models of calcium signalling are almost entirely deterministic [30].
These models typically assumed that the calcium concentration (and signalling
molecules such as IP3) within the cell can be well-approximated as
homogeneous, and therefore the dynamics can be accurately described by
ordinary differential equations for the evolution of their concentrations in
time [19]. However recent experimental evidence / mathematical analysis has
called this into question [32, 30, 50, 29, 41, 48, 44, 45]. It has been
demonstrated that (i) the interspike interval can show significant
variability, and (ii) intracellular calcium concentrations show steep
gradients, and therefore it is widely postulated that much of the emergent
phenomena are stochastic in nature [44, 45].

One of our fundamental aims is to develop an accurate microscopic model of
calcium signalling (which is thoroughly stochastic) and then use statistical
mechanical techniques to determine effective macroscopic equations. Calcium
signaling is organized in a hierarchical manner [15, 19, 23]. At the lowest
(most microscopic) level, stochastic release of $Ca^{2+}$ upon binding of IP3
receptor results in a small localized increase in cytoplasmic $Ca^{2+}$
concentration, which is often termed a ‘blip’. The IP3 receptors are typically
tightly clustered, which means that a $Ca^{2+}$ blip can stimulate the release
of additional $Ca^{2+}$ through neighboring receptors, so that the entire
cluster emits a localized ‘puff’. At the highest level of organization, if
enough puffs are generated, they can form a propagating wave of increased
$Ca^{2+}$ across an entire cell. There is considerable evidence that calcium
puffs and waves are nonlinear stochastic phenomena: it has been observed for
instance that IP3 channels can be highly active even when the average open
probability is less than half its maximum [49]. The very presence of abortive
calcium waves is indicative of the stochastic origins of the phenomena. Indeed
clusters have a diameter of $60-100$ nanometers whereas the distance between
clusters is much larger: $3-7$ micrometers [19].

There is a well-established literature that examines how spatially-distributed
calcium waves and patterns can arise out of microscopic stochastic models.
Some of the earliest work is in the Fire-Diffuse-Fire model of Keizer and
Smith [33, 32]. In this model, Ryanodine Receptors open (and release calcium)
once the ambient calcium exceeds a threshold; before closing and entering a
refractory state before being ready to open again. Keizer and Smith were able
to demonstrate that this model exhibits spatially-distributed waves. Coombes,
Hinch and Timofeeva [17, 16] developed a similar model. They estimated the
release probability for a cluster of channels and determined that it is
approximately a sigmoidal function of the local calcium concentration. Keener
[30] further extended these works by taking the opening and closing of
individual channels to be stochastic. An advantage of all of the previous
works is that the models are all spatially extended, unlike this paper.
However the model in this paper fully couples the calcium concentration to the
opening / closing of the channel. It thus allows a more detailed estimation of
how the feedback between the opening / closing of channels and the local
calcium concentration leads to calcium waves.

There have already been some efforts towards a detailed derivation of
effective macroscopic equations from a microscopic model. Hinch and Chapman
[28] used exponential asymptotic methods were employed to determine the
relative frequency of calcium sparks. Falcke et al employed approximations to
determine estimates for the approximate probability of calcium puffs (given
the frequency of blips), and also estimates for the probability of a wave
throughout the cell [41, 19].

This paper is structured as follows. In Section I, we introduce the model and
outline our main results. In Section II, we determine the Euler-Lagrange
equations that follow from the Large Deviation Principle. These can be used to
estimate optimal trajectories for first-hitting-time problems. In Section III
we apply the theory to estimate the probability of the spark-to-wave
transition in calcium dynamics. Finally in Section IV, we outline the proofs.

Notation:
Let $\left\|\cdot\right\|_{t}$ denote the supremum norm over $[0,t]$, i.e.

$\left\|z\right\|_{t}=\sup_{s\in[0,t]}|z_{s}|.$

Let $D([0,t],\mathbb{R})$ be the cadlag space of real functions from $[0,t]$
to $\mathbb{R}$ that are right-continuous and have left-hand limits. Define
$\Lambda_{t}$ to be the class of strictly increasing, continuous mappings of
$[0,t]$ onto itself (so that any $\lambda\in\Lambda_{t}$ is such that
$\lambda(0)=0$ and $\lambda(t)=t$. Following Billingsley [6], for any
$\lambda\in\Lambda$, define

$\left\|\lambda\right\|^{\circ}_{t}=\sup_{s<t}\bigg{|}\log\frac{\lambda(t)-\lambda(s)}{t-s}\bigg{|}$
(1)

Define the Skorohod Metric on $D_{t}$,

$d^{\circ}_{t}(x,y)=\inf_{\lambda\in\Lambda_{t}}\big{\\{}\left\|\lambda\right\|^{\circ}_{t}\wedge\sup_{s\in[0,t]}|x(s)-y(\lambda_{j}(s))|\big{\\}}.$
(2)

It is proved in Billingsley [6] that $d^{\circ}_{t}$ is a metric, and that
$D_{t}$ is complete and separable with respect to the topology induced by
$d^{\circ}_{t}$. Next, for $a\in\mathbb{Z}^{+}$ define the pseudo-metric
$\tilde{d}^{\circ}_{a}:\mathcal{D}([0,\infty),\mathbb{R})\times\mathcal{D}([0,\infty),\mathbb{R})\to\mathbb{R}^{+}$
to be

$\displaystyle\tilde{d}^{\circ}_{a}(x,y)=$ $\displaystyle
d^{\circ}_{a}(g_{a}x,g_{a}y)\text{ where }$ (3) $\displaystyle g_{a}(t)=$
$\displaystyle\left\\{\begin{array}[]{c}1\text{ if }t\leq a-1\\\ 0\text{ if
}t\geq a\\\ a-t\text{ if }a-1\leq t\leq a\end{array}\right.$ (7)

Then define the metric
$d^{\circ}_{\infty}:\mathcal{D}([0,\infty),\mathbb{R})\times\mathcal{D}([0,\infty),\mathbb{R})\to\mathbb{R}^{+}$
to be

$d^{\circ}_{\infty}(x,y)=\sum_{a=1}^{\infty}\inf\big{\\{}2^{-a},\tilde{d}^{\circ}_{a}(x,y)\big{\\}}.$
(8)

It is proved in Billingsley [6] that $d^{\circ}_{\infty}$ metrizes convergence
over $\mathcal{D}([0,\infty),\mathbb{R})$.

##  1 Model Outline and Large Deviations Principle

A chemical reaction network [3] consists of the following triple:

  * •

species $\mathscr{S}$, which are the chemical components whose counts we wish
to model dynamically. There are $d$ of these.

  * •

complexes $\mathscr{C}$ which are nonnegative linear combinations of the
species that describe how the species can interact and

  * •

reactions $\mathscr{R}$ which describe how to convert one complex to another.

Let $X(t)\in\mathbb{Z}_{+}^{d}$ be a vector that counts the number of
particles of each species. For a reaction $\alpha\in\mathscr{R}$, let
$Z_{\alpha}(t)\in\mathbb{Z}_{+}$ count the number of times the reaction has
occured. We thus have that

$X(t)=X(0)+\sum_{\alpha\in\mathscr{R}}Z_{\alpha}(t)\xi_{\alpha}$ (9)

where $\xi_{\alpha}\in\mathbb{Z}^{d}$ is the reaction vector corresponding to
the $\alpha^{th}$ reaction. We define the concentration vector to be the
scaled number of particles, i.e.

$\displaystyle x(t)$
$\displaystyle=N^{-1}X(t)=N^{-1}X(0)+\sum_{\alpha\in\mathscr{R}}z_{\alpha}(t)\xi_{\alpha},$
(10)

where $N\gg 1$ is a scaling parameter that indicates the system size (for
instance, the total number of particles at time $0$).

The intensity of the $\alpha^{th}$ reaction is defined to be

$N\lambda_{\alpha}\big{(}x(t),u(t)\big{)},$ (11)

where $\lambda_{\alpha}:\mathbb{R}^{d}\times\mathbb{R}^{m}\to\mathbb{R}^{+}$
is continuous and bounded. It is assumed that reactions cannot produce a
negative concentration: that is, it is assumed that if
$x\in\mathbb{R}_{+}^{d}$ and $u\in\mathbb{R}^{m}$, and $\alpha\in\mathscr{R}$
are such that $x+N^{-1}\xi_{\alpha}\notin\mathbb{R}_{+}^{d}$, then necessarily
$\lambda_{\alpha}(x,u)=0$.

The ordinary differential equation takes values in a state space
$\mathbb{R}^{m}$, and is such that

$\displaystyle\frac{du}{dt}=A\big{(}u(t),x(t)\big{)},$ (12)

where $A:\mathbb{R}^{m}\times\mathbb{R}^{d}\to\mathbb{R}^{d}$ is continuous.
The initial conditions are taken to be non-random constants i.e.

$\displaystyle u(0)$ $\displaystyle:=\hat{u}^{N}_{0}\text{ and
}x(0):=\hat{x}^{N}_{0}.$ (13)

and it is assumed that the following limits exist

$\displaystyle\lim_{N\to\infty}\hat{u}^{N}_{0}$ $\displaystyle=\hat{u}_{0}$
(14) $\displaystyle\lim_{N\to\infty}\hat{x}^{N}_{0}$
$\displaystyle=\hat{x}_{0}.$ (15)

We assume the following uniform growth bounds: for some constant $K>0$, for
all $t\geq 0$, and $u\in\mathbb{R}^{m}$, and all $\alpha\in\mathscr{R}$,

$\displaystyle\left\|A(u,X)\right\|$ $\displaystyle\leq
K\big{(}1+\left\|u\right\|+\left\|X\right\|\big{)}$ (16)
$\displaystyle\big{|}\lambda_{\alpha}\big{(}X,u\big{)}\big{|}$
$\displaystyle\leq K.$ (17)

Second, we assume that $A$ and $\lambda_{n}$ are both uniformly Lipschitz in
their arguments: that is,

$\displaystyle\left\|A\big{(}u,X\big{)}-A\big{(}\tilde{u},\tilde{X}\big{)}\right\|$
$\displaystyle\leq
C\big{\\{}\left\|u-\tilde{u}\right\|+\left\|X-\tilde{X}\right\|\big{\\}}$ (18)
$\displaystyle\sup_{\alpha\in\mathscr{R}}\|\lambda_{\alpha}\big{(}u,X\big{)}-\lambda_{\alpha}\big{(}\tilde{u},\tilde{X}\big{)}\|$
$\displaystyle\leq
C\big{\\{}\left\|u-\tilde{u}\right\|+\left\|X-\tilde{X}\right\|\big{\\}}.$
(19)

Standard ODE theory dictates that for any piecewise-continuous $x$, there is a
unique solution to the ODE in (12).

The linearity of Poisson Processes ensures that $Z_{\alpha}(t)$ can be
represented as, writing $\\{Y_{\alpha}(t)\\}_{\alpha\in\mathscr{R}}$ to be
independent unit-intensity Poisson Processes [3],

$\displaystyle Z_{\alpha}(t)=$ $\displaystyle
Y_{\alpha}\bigg{(}N\int_{0}^{t}\lambda_{\alpha}\big{(}x(s),u(s)\big{)}ds\bigg{)}$
(20) $\displaystyle z_{\alpha}(t)=$ $\displaystyle N^{-1}Z_{\alpha}(t)\text{
with initial condition}$ (21) $\displaystyle z_{\alpha}(0)=$ $\displaystyle
0.$ (22)

Writing $y_{\alpha}(t)=N^{-1}Y_{\alpha}(Nt)$, it may be observed that

$z_{\alpha}(t)=y_{\alpha}\bigg{(}\int_{0}^{t}\lambda_{\alpha}\big{(}x(s),u(s)\big{)}ds\bigg{)}.$
(23)

The main result of this paper is the following theorem. Further below in (27)
we specify $\Upsilon$ to be the state space that the triple $(z,x,u)$ lives in
(and we define a topology too). The sequence of probability laws of
$\big{(}z,x,u\big{)}$ satisfies the following asymptotic estimate on the space
$\Upsilon$ as $N\to\infty$.

######  Theorem 1.1.

Supponse that $\mathcal{O},\mathcal{A}\subset\Upsilon$, with $\mathcal{O}$
open and $\mathcal{A}$ closed, are such that for any $T>0$,

$\displaystyle\inf_{(z,x,u)\in\mathcal{A}\cup\mathcal{O}}\inf_{\alpha\in\mathscr{R}}\inf_{s\in[0,T]}\lambda_{\alpha}(x(s),u(s))>0.$
(24)

Then there exists a function $\mathcal{J}:\Upsilon\to\mathbb{R}$ (specified in
Section 1.1) that is (i) lower-semi-continuous and (ii) has compact level
sets, such that

$\displaystyle\underset{N\to\infty}{\overline{\lim}}N^{-1}\log\mathbb{P}\big{(}(z,x,u)\in\mathcal{A}\big{)}$
$\displaystyle\leq-\inf_{(z,x,u)\in\mathcal{A}}\mathcal{J}(z,x,u)$ (25)
$\displaystyle\underset{N\to\infty}{\underline{\lim}}N^{-1}\log\mathbb{P}\big{(}(z,x,u)\in\mathcal{O}\big{)}$
$\displaystyle\geq-\inf_{(z,x,u)\in\mathcal{O}}\mathcal{J}(z,x,u).$ (26)

###  1.1 Topological Definitions

Before we state our main result we must briefly note the topological space
that our variables belong to. Write the state space for $\big{(}z,x,u\big{)}$
as $\Upsilon$, i.e.

$\Upsilon:=\big{\\{}(z,x,u)\in\tilde{\mathcal{D}}([0,\infty),\mathbb{R}^{+})^{M}\times\mathcal{D}([0,\infty),\mathbb{R}^{+})^{d}\times\mathcal{C}^{1}([0,\infty),\mathbb{R})^{m}\\\
\text{ where
(i)}\sum_{\alpha\in\mathscr{R}}z_{\alpha}(t)\zeta_{\alpha}=x(t),\;\text{(ii)
for all }t\geq 0,\;\;\frac{du}{dt}=A\big{(}u(t),x(t)\big{)}\text{ and
}u(0)=\hat{u}\big{\\}}$ (27)

and endow $\Upsilon$ with the product topology. Since the values of $x$ and
$u$ are effectively determined by $z$, we only need to metrize the convergence
of $z$. This is captured by the following lemma.

######  Lemma 1.2.

$\Upsilon$ is a closed subset of
$\tilde{\mathcal{D}}([0,\infty),\mathbb{R}^{+})^{M}\times\mathcal{D}([0,\infty),\mathbb{R}^{+})^{d}\times\mathcal{C}^{1}([0,\infty),\mathbb{R})^{m}$.
Also, $\Upsilon$ is separable and metrizable. Furthermore for any $t>0$, there
exists $c_{t}>0$ such that for all $\delta>0$,

$\big{\\{}(z,x,u),(\tilde{z},\tilde{x},\tilde{u})\in\Upsilon\times\Upsilon:\sup_{\alpha\in\mathscr{R}}d^{\circ}_{t}(z_{\alpha},\tilde{z}_{\alpha})\leq\delta\big{\\}}\subseteq\big{\\{}(z,x,u),(\tilde{z},\tilde{x},\tilde{u})\in\Upsilon\times\Upsilon:d^{\circ}_{t}(z,\tilde{z})\leq\delta\text{
and }\\\ \sup_{0\leq s\leq t}\left\|x(s)-\tilde{x}(s)\right\|\leq
c_{t}\delta\text{ and }\sup_{0\leq s\leq
t}\left\|u(s)-\tilde{u}(s)\right\|\leq c_{t}\delta\big{\\}}$ (28)

We can now define the rate function

$\displaystyle\mathcal{J}(z,x,u):\Upsilon\mapsto\mathbb{R}^{+},$ (29)

starting by stipulating that it is infinite in the case that
$z_{\alpha}:\mathbb{R}^{+}\to\mathbb{R}^{+}$ is not absolutely continuous for
any $\alpha\in\mathscr{R}$. Otherwise, if $z_{\alpha}$ is absolutely
continuous for every $\alpha\in\mathscr{R}$ (which means that it must have a
derivative $\dot{z}_{\alpha}$ for Lebesgue almost every time), for any
realization of $z$, define the set

$\mathfrak{V}(z,x,u)=\big{\\{}t\in\mathbb{R}^{+}:\text{For any
}\alpha\in\mathscr{R},\;\dot{z}_{\alpha}(t)\neq 0\text{ and
}\lambda_{\alpha}(x(t),u(t))=0\big{\\}}$ (30)

and define the rate function

$\displaystyle\mathcal{J}\big{(}z,x,u\big{)}$ $\displaystyle=\infty\text{ if
}\mathfrak{V}(z,x,u)\text{ is of nonzero Lebesgue Measure, else otherwise }$
(31)
$\displaystyle=\sum_{\alpha\in\mathscr{R}}\int_{\mathbb{R}^{+}/\mathfrak{V}(z,x,u)}\ell\bigg{(}\dot{z}_{\alpha}(r)/\lambda_{\alpha}(x(r),u(r))\bigg{)}\lambda_{\alpha}(x(r),u(r))dr\text{
where }$ (32) $\displaystyle\ell(a)$ $\displaystyle=a\log a-a+1.$ (33)

##  2 Euler-Lagrange Equations

In this section we determine the structure of the Euler-Lagrange equations
that determine the most likely trajectory followed by the system in attaining
a specified state. Computing estimates for first-hitting-times is one of the
most important applications of Large Deviations Theory. For the application of
Large Deviations Theory to first-hitting-time estimates, see for example [21,
24, 40]. The ‘fast processes’ (i.e. the ODE in $u(t)$) are implemented as a
constraint (with Lagrange multipliers). We start by framing the problem in
terms of the reaction fluxes; that is, we we wish to determine the most likely
trajectory followed by the system in attaining a target reaction flux
$z^{*}\in(\mathbb{R}^{+})^{|\mathscr{R}|}$ at time $T$. We must determine the
trajectory

$(z,x,u)\in\mathcal{C}^{1}(\mathbb{R})^{M}\times\mathcal{C}^{1}(\mathbb{R})^{d}\times\mathcal{C}^{1}(\mathbb{R})^{m}$

such that $\mathcal{J}_{T}$ is minimized, where

$\displaystyle\mathcal{J}_{T}=\sum_{\alpha\in\mathscr{R}}\int_{0}^{T}L\big{(}\dot{z}(s),x(s),u(s)\big{)}ds,$
(34)

where $L:\mathbb{R}^{M}\times\mathbb{R}^{d}\times\mathbb{R}^{m}\
\to\mathbb{R}^{+}$ is the integrand of the rate function in (32), i.e.

$\displaystyle L\big{(}q,x,u\big{)}$
$\displaystyle=\sum_{\alpha\in\mathscr{R}}\ell\big{(}q_{\alpha}/\lambda_{\alpha}(x,u)\big{)}\lambda_{\alpha}(x,u)$
(35) $\displaystyle\ell(a)$ $\displaystyle=a\log a-a+1.$ (36)

subject to the following constraints and boundary conditions

$\displaystyle\frac{du}{dt}=$ $\displaystyle A(u,x)$ (37) $\displaystyle
x(t)=$ $\displaystyle
x(0)+\sum_{\alpha\in\mathscr{R}}z_{\alpha}(t)\xi_{\alpha}.$ (38)
$\displaystyle u(0)=$ $\displaystyle\hat{u}$ (39) $\displaystyle
z_{\alpha}(0)=$ $\displaystyle 0$ (40) $\displaystyle z_{\alpha}(T)=$
$\displaystyle z_{\alpha}^{*}.$ (41)

We introduce Lagrange Multipliers $\\{\eta^{i}(t)\\}_{1\leq i\leq m}$
corresponding to the $m$ constraints

$\displaystyle\frac{du^{i}_{t}}{dt}=A^{i}(u_{t},x_{t}).$ (42)

Our task is now to find the critical points of the functional

$\displaystyle\tilde{\mathcal{J}}_{T}(z,u,\eta)=$
$\displaystyle\int_{0}^{T}\tilde{L}\big{(}\dot{z}(s),x(s),u(s)\big{)}ds\text{
where }$ (43) $\displaystyle\tilde{L}\big{(}\dot{z}(s),x(s),u(s)\big{)}=$
$\displaystyle
L\big{(}\dot{z}(s),x(s),u(s)\big{)}-\sum_{i=1}^{m}\eta^{i}(s)\big{\\{}\dot{u}^{i}(s)-A^{i}(u(s),x(s))\big{\\}}\text{
such that }$ (44) $\displaystyle x(t)=$
$\displaystyle\hat{x}(0)+\sum_{\alpha\in\mathscr{R}}z_{\alpha}(t)\xi_{\alpha}.$
(45)

We now take the Frechet Derivative of (43), i.e. for
$w\in\mathcal{C}^{1}([0,T],\mathbb{R})^{M}$ and
$v\in\mathcal{C}^{1}([0,T],\mathbb{R})^{m}$, and
$\kappa\in\mathcal{C}^{1}([0,T],\mathbb{R})^{m}$, with boundary conditions
$w(0)=0,w(T)=0,v(0)=0$, and the boundary conditions for $\kappa$ are not yet
clear. Thanks to the constraint in (45), we have that

$\displaystyle\frac{d}{dz_{\alpha}}=\sum_{i=1}^{d}\frac{\partial}{\partial
x^{i}}\xi^{i}_{\alpha}.$ (46)

We therefore find that

$D\tilde{\mathcal{J}}_{T}(z,u,\eta)\cdot(w,v,\kappa):=\lim_{\epsilon\to
0^{+}}\epsilon^{-1}\big{\\{}\tilde{\mathcal{J}}_{T}(z+\epsilon w,u+\epsilon
v,\eta+\epsilon\kappa)-\tilde{\mathcal{J}}_{T}(z,u,\eta)\big{\\}}\\\
=\sum_{\alpha\in\mathscr{R}}\int_{0}^{T}w_{\alpha}(t)\bigg{\\{}\sum_{i=1}^{d}\bigg{(}\frac{\partial
L}{\partial x^{i}}\xi^{i}_{\alpha}+\sum_{k=1}^{m}\eta^{k}\frac{\partial
A^{k}}{\partial x^{i}}\xi^{i}_{\alpha}\bigg{)}-\frac{d}{dt}\frac{\partial
L}{\partial\dot{z}_{\alpha}}\bigg{\\}}dt\\\
+\sum_{i=1}^{m}\int_{0}^{T}v^{i}(t)\bigg{\\{}\frac{\partial L}{\partial
u^{i}}+\frac{d\eta^{i}}{dt}+\sum_{k=1}^{m}\eta^{k}\frac{\partial
A^{k}}{\partial
u^{i}}\bigg{\\}}-\kappa^{i}\bigg{\\{}\frac{du^{i}}{dt}-A^{i}(u(t),x(t))\bigg{\\}}dt\\\
+\sum_{\alpha\in\mathscr{R}}\bigg{\\{}w_{\alpha}(T)\frac{\partial}{\partial\dot{z}_{\alpha}}L(\dot{z}(T),x(T),u(T))-w_{\alpha}(0)\frac{\partial}{\partial\dot{z}_{\alpha}}L(\dot{z}(0),x(0),u(0))\bigg{\\}}\\\
-\sum_{i=1}^{m}\big{\\{}\eta^{i}(T)v^{i}(T)-\eta^{i}(0)v^{i}(0)\big{\\}},$
(47)

and to obtain the above expression, we have used integration by parts to write

$\sum_{\alpha\in\mathscr{R}}\int_{0}^{T}\dot{w}_{\alpha}(t)\frac{\partial
L}{\partial\dot{z}_{\alpha}}dt=\sum_{\alpha\in\mathscr{R}}\bigg{\\{}w_{\alpha}(T)\frac{\partial}{\partial\dot{z}_{\alpha}}L(\dot{z}(T),x(T),u(T))-w_{\alpha}(0)\frac{\partial}{\partial\dot{z}_{\alpha}}L(\dot{z}(0),x(0),u(0))\\\
-\int_{0}^{T}w_{\alpha}(t)\frac{d}{dt}\frac{\partial
L}{\partial\dot{z}_{\alpha}}dt\bigg{\\}}$ (48)

and

$\sum_{i=1}^{m}\int_{0}^{t}\eta^{i}(t)\dot{v}^{i}(t)dt=\sum_{i=1}^{m}\big{\\{}\eta^{i}(T)v^{i}(T)-\eta^{i}(0)v^{i}(0)-\int_{0}^{t}\dot{\eta}^{i}(t)v^{i}(t)dt\big{\\}}.$
(49)

Setting the coefficients of $w_{\alpha}$, $v$ and $\kappa^{i}$ to zero, we
obtain the equations

$\displaystyle\sum_{i=1}^{d}\bigg{(}\frac{\partial L}{\partial
x^{i}}\xi^{i}_{\alpha}+\sum_{k=1}^{m}\eta^{k}\frac{\partial A^{k}}{\partial
x^{i}}\xi^{i}_{\alpha}\bigg{)}-\frac{d}{dt}\frac{\partial
L}{\partial\dot{z}_{\alpha}}$ $\displaystyle=0$ (50)
$\displaystyle\frac{\partial L}{\partial
u^{i}}+\frac{d\eta^{i}}{dt}+\sum_{k=1}^{m}\eta^{k}\frac{\partial
A^{k}}{\partial u^{i}}$ $\displaystyle=0$ (51)
$\displaystyle\frac{du^{i}}{dt}-A^{i}(u(s),x(s))$ $\displaystyle=0$ (52)
$\displaystyle\eta^{i}(T)$ $\displaystyle=0.$ (53)

Noting that $\dot{\ell}(a)=\log a$, we compute that $\frac{\partial
L}{\partial\dot{z}_{\alpha}}=\log\frac{\dot{z}_{\alpha}}{\lambda_{\alpha}(x(t),u(t))}$,
and therefore

$\displaystyle\frac{d}{dt}\frac{\partial L}{\partial\dot{z}_{\alpha}}=$
$\displaystyle\frac{d}{dt}\left(\log\frac{\dot{z}_{\alpha}}{\lambda_{\alpha}(x(t),u(t))}\right)$
(54) $\displaystyle=$
$\displaystyle\frac{\ddot{z}_{\alpha}}{\dot{z}_{\alpha}}-\frac{1}{\lambda_{\alpha}(x(t),u(t))}\frac{d}{dt}\lambda_{\alpha}(x(t),u(t)).$
(55)

Thus in the case that $\lambda_{\alpha}(x,u)>0$, the Euler Lagrange equations
are such that for $1\leq j\leq m$, $\alpha\in\mathscr{R}$ and $1\leq i\leq d$,

$\displaystyle\frac{\ddot{z}_{\alpha}}{\dot{z}_{\alpha}}=$
$\displaystyle\frac{1}{\lambda_{\alpha}}\frac{d\lambda_{\alpha}}{dt}+\sum_{i=1}^{d}\xi^{i}_{\alpha}\frac{\partial
L}{\partial x_{i}}-\sum_{k=1}^{d}\sum_{j=1}^{m}\eta^{j}(t)\frac{\partial
A^{j}}{\partial x^{k}_{t}}\xi^{k}_{\alpha}$ (56)
$\displaystyle\frac{d\eta^{j}}{dt}=$ $\displaystyle-\frac{\partial L}{\partial
u^{j}}-\sum_{k=1}^{m}\eta^{k}\frac{\partial A^{k}}{\partial u^{j}}$ (57)
$\displaystyle x^{i}(t)=$ $\displaystyle
x^{i}(0)+\sum_{\alpha\in\mathscr{R}}z_{\alpha}(t)\xi^{i}_{\alpha}\text{ for
all }t\geq 0.$ (58) $\displaystyle\eta^{j}(T)=$ $\displaystyle 0$ (59)
$\displaystyle z_{\alpha}(0)=$ $\displaystyle 0$ (60) $\displaystyle z(T)=$
$\displaystyle z_{*}.$ (61)

For the sake of completeness, we note that

$\displaystyle\frac{d\lambda_{\alpha}}{dt}=$
$\displaystyle\sum_{\beta\in\mathscr{R}}\sum_{i=1}^{d}\frac{\partial\lambda_{\alpha}}{\partial
x^{i}}\xi^{i}_{\beta}\dot{z}_{\beta}+\sum_{i=1}^{m}\frac{\partial\lambda_{\alpha}}{\partial
u^{i}}A^{i}(u(t),x(t))$ (62) $\displaystyle\frac{\partial L}{\partial u^{i}}=$
$\displaystyle\sum_{\alpha\in\mathscr{R}}\frac{\partial\lambda_{\alpha}}{\partial
u^{i}}\bigg{(}1-\frac{\dot{z}_{\alpha}}{\lambda_{\alpha}}\bigg{)}$ (63)

It can be seen that this system amounts to an ODE boundary value problem. One
way that one may attempt to solve this problem is via the shooting method
[21]. At time $0$, the unknowns are
$\\{\dot{z}_{\alpha}(0)\\}_{\alpha\in\mathscr{R}}$ and
$\\{\eta^{j}(0)\\}_{1\leq j\leq m}$. The shooting method requires one to guess
values for the unknowns at time $0$, (ii) integrate forward the equations to
time $T$, and then (iii) implement the constraints
$\\{\dot{z}_{\alpha}(T)=\dot{z}_{*,\alpha}\\}_{\alpha\in\mathscr{R}}$ and
$\\{\eta^{j}(T)=0\\}_{1\leq j\leq m}$.

###  2.1 Simplified Euler-Lagrange Equations

Usually, one desires to know the first-hitting-time for the concentrations
$\\{x^{i}(T)\\}_{1\leq i\leq d}$, rather than the first hitting time for the
reaction fluxes. That is, for a fixed $x_{*}\in\mathbb{R}^{d}$, one wishes to
estimate the probability that $x_{T}\simeq x_{*}$, and also determine the
optimal (most likely) path followed by the system in attaining this point. For
this reason, it can be simpler to eliminate the reaction fluxes and reduce the
problem to one purely in terms of the concentrations and the ODE variables
$\\{u^{i}(T)\\}_{1\leq i\leq m}$. To this end, for any
$\dot{x}\in\mathbb{R}^{d}$, $x\in(\mathbb{R}^{+})^{d}$, and
$u\in\mathbb{R}^{m}$, define the set

$\displaystyle\Xi(\dot{x},x,u)=\big{\\{}\dot{z}\;:\;\dot{x}^{p}=\sum_{\alpha\in\mathscr{R}}\xi^{p}_{\alpha}\dot{z}_{\alpha}\big{\\}}\subset(\mathbb{R}^{+})^{M}$
(64)

Then define the function

$\hat{L}\big{(}\dot{x}(s),x(s),u(s)\big{)}=\inf\big{\\{}L(\dot{z},x(s),u(s))\;:\;\dot{z}\in\Xi(\dot{x},x,u)\big{\\}}$
(65)

In the case that $\Xi(\dot{x},x,u)=\emptyset$, define
$\hat{L}\big{(}\dot{x}(s),x(s),u(s)\big{)}=\infty$. One can thus define the
contracted rate function

$\displaystyle\tilde{\mathcal{J}}_{T}:$
$\displaystyle\mathcal{C}_{ac}([0,T],(\mathbb{R}^{+})^{d}\big{)}\times\mathcal{C}_{ac}([0,T],\mathbb{R}^{m}\big{)}\to\mathbb{R}$
(66) $\displaystyle\tilde{\mathcal{J}}_{T}(x,u)$
$\displaystyle=\int_{0}^{T}\tilde{L}\big{(}\dot{x}(s),x(s),u(s)\big{)}ds,$
(67)

and it is immediate from the Contraction Principle [18] that
$\tilde{\mathcal{J}}_{T}$ governs the Large Deviations of $(x,u)$. The
following convexity property makes it easier to compute optimal trajectories.

######  Lemma 2.1.

$\tilde{L}$ is convex in its first argument.

###### Proof.

This is because $\Xi$ is linear in $\dot{x}$, and $L$ is convex in its first
argument. ∎

Next, we outline the Euler-Lagrange equations for $\tilde{\mathcal{J}}_{T}$.
The derivation parallels that of the previous section. Our task is to find the
critical points of the functional

$\displaystyle\grave{\mathcal{J}}_{T}(x,u,\eta)=$
$\displaystyle\int_{0}^{T}\grave{L}\big{(}\dot{x}(s),x(s),u(s)\big{)}ds\text{
where }$ (68)
$\displaystyle\grave{L}\big{(}\dot{x}(s),x(s),u(s),\eta(s)\big{)}=$
$\displaystyle\tilde{L}\big{(}\dot{x}(s),x(s),u(s)\big{)}-\sum_{i=1}^{m}\eta^{i}(s)\big{\\{}\dot{u}^{i}(s)-A^{i}(u(s),x(s))\big{\\}},$
(69)

and $\eta^{i}\in\mathcal{C}^{1}([0,T],\mathbb{R})$ is the Lagrange Multiplier.

We now take the Frechet Derivative of (68), i.e. for
$w\in\mathcal{C}^{1}([0,T],\mathbb{R})^{M}$ and
$v\in\mathcal{C}^{1}([0,T],\mathbb{R})^{m}$, and
$\kappa\in\mathcal{C}^{1}([0,T],\mathbb{R})^{m}$, with boundary conditions
$w(0)=0,w(T)=0,v(0)=0$, and the boundary conditions for $\kappa$ are not yet
clear,

$D\grave{\mathcal{J}}_{T}(x,u,\eta)\cdot(w,v,\kappa):=\lim_{\epsilon\to
0^{+}}\epsilon^{-1}\big{\\{}\grave{\mathcal{J}}_{T}(x+\epsilon w,u+\epsilon
v,\eta+\epsilon\kappa)-\tilde{\mathcal{J}}_{T}(x,u,\eta)\big{\\}}\\\
=\int_{0}^{T}\sum_{i=1}^{d}w_{i}(t)\bigg{\\{}\bigg{(}\frac{\partial\hat{L}}{\partial
x^{i}}+\sum_{k=1}^{m}\eta^{k}\frac{\partial A^{k}}{\partial
x^{i}}\bigg{)}-\frac{d}{dt}\frac{\partial
L}{\partial\dot{x}^{i}}\bigg{\\}}dt\\\
+\sum_{i=1}^{m}\int_{0}^{T}v^{i}(t)\bigg{\\{}\frac{\partial\hat{L}}{\partial
u^{i}}+\frac{d\eta^{i}}{dt}+\sum_{k=1}^{m}\eta^{k}\frac{\partial
A^{k}}{\partial
u^{i}}\bigg{\\}}-\kappa^{i}\bigg{\\{}\frac{du^{i}}{dt}-A^{i}(u(t),x(t))\bigg{\\}}dt\\\
-\sum_{i=1}^{m}\big{\\{}\eta^{i}(T)v^{i}(T)-\eta^{i}(0)v^{i}(0)\big{\\}}$ (70)

The Euler-Lagrange equations are thus, for $1\leq i\leq d$ and $1\leq j\leq
m$,

$\displaystyle\frac{\partial\hat{L}}{\partial
x^{i}}+\sum_{k=1}^{m}\eta^{k}\frac{\partial A^{k}}{\partial
x^{i}}-\frac{d}{dt}\frac{\partial\hat{L}}{\partial\dot{x}^{i}}$
$\displaystyle=0$ (71) $\displaystyle\frac{\partial\hat{L}}{\partial
u^{j}}+\frac{d\eta^{j}}{dt}+\sum_{k=1}^{m}\eta^{k}\frac{\partial
A^{k}}{\partial u^{j}}$ $\displaystyle=0$ (72)
$\displaystyle\frac{du^{i}}{dt}-A^{i}(u(s),x(s))$ $\displaystyle=0$ (73)
$\displaystyle\eta^{i}(T)$ $\displaystyle=0.$ (74)

For the sake of completeness, we note that the derivative assumes the form

$\displaystyle\frac{d}{dt}\frac{\partial\hat{L}}{\partial\dot{x}^{i}}=\sum_{j=1}^{d}\bigg{\\{}\frac{\partial^{2}\tilde{L}}{\partial\dot{x}^{i}\partial\dot{x}^{j}}\ddot{x}^{j}+\frac{\partial^{2}\hat{L}}{\partial\dot{x}^{i}\partial
x^{j}}\dot{x}^{j}\bigg{\\}}+\sum_{j=1}^{m}\frac{\partial^{2}\hat{L}}{\partial\dot{x}^{i}\partial
u^{j}}\frac{du^{j}}{dt}$ (75)

##  3 Application to Calcium Signalling

We apply this theory to determine estimates for the typical time it takes for
cells with stochastically opening and closing calcium channels to change from
most of the channels being closed to most being open [19]. There are thus two
components to the model: the opening and closing channels, and the calcium
concentration (that evolves deterministically between the opening and closing
of the channels). The model is simplified to be such that the calcium
concentration throughout the cell is spatially homogeneous.

To begin with, we must outline a simple microscopic model of stochastic
calcium effects. Following for instance [17, 30, 44, 45], we employ a hybrid
model: with the calcium diffusion modeled deterministically, the IP3
concentration constant throughout the cell, and the channel opening and
closing modeled stochastically.

Our basic aim is to determine the probability of calcium puffs and waves by
studying a similar stochastic model to [17, 30, 44, 45]. These papers
performed detailed numerical simulations, and analytically estimated the
probabilities. To calculate the probability of a calcium puff in a cluster,
the calcium and IP3 concentration throughout the cluster is assumed to be
homogeneous. However there is a feedback effect of the opening / closing of
the channels: when these bind or unbind the calcium / IP3, they alter the
local concentration. Thus the opening and closing of channels within a cluster
is not independent because the channels communicate with each other via the
calcium / IP3 concentrations.

###  3.1 Piecewise-Deterministic-Markov-Process Model of Calcium Dynamics

We suppose that there are $N$ channels distributed on the cell membrane. It is
assumed that spatial effects are negligible, so that the calcium
concentrations can be modelled as approximately spatially homogeneous using
ODEs. Let $u_{1}$ be the calcium concentration in the cellular cytosol, and
$u_{2}$ be the calcium concentration in the endoplasmic reticulum.

We take the number of stochastic variables corresponding to each channel to be
$1$. That is, $Z^{i}_{t}=1$ if the channel is open, and $Z^{i}_{t}=0$ if the
channel is closed. Employing a standard model [19] for the calcium dynamics,
we obtain that

$\displaystyle\frac{du_{1}}{dt}=$ $\displaystyle
k_{f}x(t)(u_{2}-u_{1})-J_{serca}(u)+J_{leak}:=A^{1}(u,x)$ (76)
$\displaystyle\frac{du_{2}}{dt}=$ $\displaystyle\gamma J_{serca}(u_{1})-\gamma
k_{f}x(t)(u_{2}-u_{1})-\gamma J_{leak}:=A^{2}(u,x),$ (77) $\displaystyle
J_{serca}(u_{1})=$
$\displaystyle\frac{V_{s}(u_{1})^{2}}{K_{s}^{2}+(u_{1})^{2}}$ (78)
$\displaystyle J_{leak}(u)=$ $\displaystyle k_{leak}(u_{2}-u_{1})$ (79)

One may observe that

$\displaystyle\frac{d}{dt}\big{(}\gamma u_{1}+u_{2}\big{)}=0,$ (80)

which means that we can eliminate $u_{2}$ from the dynamics, and write
$u:=u_{1}(t)$, and $A(t):=A^{1}(u,x)$. Now

$\displaystyle x(t)=\frac{1}{N}\sum_{j=1}^{N}Z^{j}_{t}.$ (81)

The Markovian switching of the channel is given by

  * •

Reaction 1: $1$ to $0$ with intensity $\lambda_{1}(u,x)=\alpha_{-1}x$ (a
constant). The stoichiometric constant is $\zeta_{1}=-1$.

  * •

Reaction 2: $0$ to $1$ with intensity $\lambda_{2}(u,x)=\alpha_{1}u_{1}(1-x)$.
The stoichiometric constant is $\zeta_{2}=1$.

Let $z_{1}(t)$ count the number of times reaction $1$ happens, and $z_{2}(t)$
count the number of times that reaction $2$ happens. Then it must be that

$\displaystyle x(t)=x(0)+z_{2}(t)-z_{1}(t).$ (82)

#### Large $N$ Limiting Dynamics.

As $N\to\infty$ it is classical that the concentration of open channels
converges to the following ODE [38].

$\displaystyle\frac{dz}{dt}=$
$\displaystyle\lambda_{2}(u,x)-\lambda_{1}(u,x).$ (83)
$\displaystyle\frac{du_{1}}{dt}=$ $\displaystyle
k_{f}x(t)(u_{2}-u_{1})-J_{serca}(u)+J_{leak}:=A^{1}(u,x)$ (84)
$\displaystyle\frac{du_{2}}{dt}=$ $\displaystyle\gamma J_{serca}(u_{1})-\gamma
k_{f}x(t)(u_{2}-u_{1})-\gamma J_{leak}:=A^{2}(u,x),$ (85)

It turns out that there exists a unique fixed point of (83) - (85). This point
is written as $(z^{*},u_{1}^{*},u_{2}^{*})$.

#### Euler-Lagrange Equations

It seems easiest to work with the contracted rate function in (65). We first
find a much simpler form for the contracted rate function.

######  Lemma 3.1.

$\displaystyle\hat{L}\big{(}\dot{x},x,u\big{)}=\ell\bigg{(}\frac{\dot{z}_{1}}{\lambda_{1}(x,u)}\bigg{)}\lambda_{1}(x,u)+\ell\bigg{(}\frac{\dot{x}+\dot{z}_{1}}{\lambda_{2}(x,u)}\bigg{)}\lambda_{2}(x,u)$
(86)

where $\dot{z}_{1}$ is such that

$\displaystyle\dot{z}_{1}\big{(}\dot{x}+\dot{z}_{1}\big{)}=\lambda_{1}(x,u)\lambda_{2}(x,u).$
(87)

One must make sure to choose the root of the above quadratic that is such that
$\dot{z}_{1}\geq 0$.

###### Proof.

It must be that

$\dot{x}=\dot{z}_{2}-\dot{z}_{1}.$

Thus, upon substitution, the Lagrangian assumes the form in (86). We must
choose $\dot{z}_{1}$ to minimize (86), subject to the constraint
$0\leq\dot{z}_{1}\leq\dot{x}$. Differentiating (86) with respect to
$\dot{z}_{1}$, and setting the derivative to zero, we obtain that

$\displaystyle\log\bigg{(}\frac{\dot{z}_{1}}{\lambda_{1}(x,u)}\bigg{)}+\log\bigg{(}\frac{\dot{z}_{1}+\dot{x}}{\lambda_{2}(x,u)}\bigg{)}=0.$
(88)

Exponentiating both sides, we obtain (87). ∎

Adapted to the calcium-signalling model of this section, the original Euler-
Lagrange equations in (71)-(74) are thus, for $i\in\\{1,2\\}$,

$\displaystyle\frac{\partial\hat{L}}{\partial
x}+\sum_{k=1}^{2}\eta^{k}\frac{\partial A^{k}}{\partial x}$
$\displaystyle=\frac{d}{dt}\frac{\partial\hat{L}}{\partial\dot{x}}$ (89)
$\displaystyle\frac{\partial\hat{L}}{\partial
u}+\frac{d\eta^{i}}{dt}+\sum_{k=1}^{2}\eta^{k}\frac{\partial A^{k}}{\partial
u}$ $\displaystyle=0$ (90) $\displaystyle\frac{du}{dt}-A(u(t),x(t))$
$\displaystyle=0$ (91) $\displaystyle\eta^{i}(T)$ $\displaystyle=0.$ (92)

Upon making substitutions, (89) leads to a second order ODE of the form
$\frac{d^{2}x}{dt^{2}}=\ldots$. Next, we compute that the first order
derivatives are

######  Lemma 3.2.

$\displaystyle\frac{\partial\hat{L}}{\partial\dot{x}}$
$\displaystyle=\log\bigg{(}\frac{\dot{z}_{1}}{\lambda_{1}(x,u)}\bigg{)}\frac{\partial\dot{z}_{1}}{\partial\dot{x}}+\log\bigg{(}\frac{\dot{x}+\dot{z}_{1}}{\lambda_{2}(x,u)}\bigg{)}\bigg{(}1+\frac{\partial\dot{z}_{1}}{\partial\dot{x}}\bigg{)}\text{
where }$ (93) $\displaystyle\frac{\partial\dot{z}_{1}}{\partial\dot{x}}$
$\displaystyle=-\frac{\dot{z}_{1}}{2\dot{z}_{1}+\dot{x}}.$ (94)

Notice also that the derivative assumes the form

$\displaystyle\frac{d}{dt}\frac{\partial\hat{L}}{\partial\dot{x}}=\frac{\partial^{2}\hat{L}}{\partial\dot{x}^{2}}\ddot{x}+\frac{\partial^{2}\hat{L}}{\partial\dot{x}\partial
x}\dot{x}+\sum_{j=1}^{2}\frac{\partial^{2}\hat{L}}{\partial\dot{x}\partial
u_{j}}\frac{du_{j}}{dt}$ (95)

In the next lemma we collect some more identities. The proof is neglected (one
just applies Calculus to compute the derivatives).

######  Lemma 3.3.

The second derivatives are

$\displaystyle\frac{\partial^{2}\dot{z}_{1}}{\partial\dot{x}^{2}}$
$\displaystyle=\frac{\dot{z}_{1}}{(2\dot{z}_{1}+\dot{x})^{2}}\bigg{(}1+\frac{\dot{x}}{2\dot{z}_{1}+\dot{x}}\bigg{)}$
(96) $\displaystyle\frac{\partial^{2}\hat{L}}{\partial\dot{x}^{2}}$
$\displaystyle=2\frac{\partial^{2}\dot{z}_{1}}{\partial\dot{x}^{2}}\log\bigg{(}\frac{\dot{z}_{1}(\dot{x}+\dot{z}_{1})}{\lambda_{1}(x,u)\lambda_{2}(x,u)}\bigg{)}+\frac{\lambda_{1}(x,u)}{\dot{x}+\dot{z}_{1}}\bigg{(}1+\frac{\partial\dot{z}_{1}}{\partial\dot{x}}\bigg{)}$
(97) $\displaystyle\frac{\partial^{2}\hat{L}}{\partial\dot{x}\partial x}$
$\displaystyle=-\frac{\partial\lambda_{1}(x,u)}{\partial
x}\frac{1}{\lambda_{1}(x,u)}\frac{\partial\dot{z}_{1}}{\partial\dot{x}}-\frac{\partial\lambda_{2}(x,u)}{\partial
x}\frac{1}{\lambda_{2}(x,u)}\bigg{(}1+\frac{\partial\dot{z}_{1}}{\partial\dot{x}}\bigg{)}$
(98) $\displaystyle\frac{\partial^{2}\hat{L}}{\partial\dot{x}\partial u_{i}}$
$\displaystyle=-\frac{\partial\lambda_{1}(x,u)}{\partial
u_{i}}\frac{1}{\lambda_{1}(x,u)}\frac{\partial\dot{z}_{1}}{\partial\dot{x}}-\frac{\partial\lambda_{2}(x,u)}{\partial
u_{i}}\frac{1}{\lambda_{2}(x,u)}\bigg{(}1+\frac{\partial\dot{z}_{1}}{\partial\dot{x}}\bigg{)}$
(99)

###  3.2 Optimal Trajectories

We wish to predict the probability of the proportion of open calcium channels
being significantly greater than its equilibrium value at some time $T$. We
write the ‘target’ proportion of channels to be $\hat{x}\in(x_{*},1]$. We thus
wish to compute the trajectory
$(x,u)\in\mathcal{A}\subset\mathcal{C}^{2}([0,T],\mathbb{R})\times\mathcal{C}^{1}([0,T],\mathbb{R}^{2})$
such that

$\displaystyle\tilde{\mathcal{J}}_{T}(x,u)=\inf\big{\\{}\tilde{\mathcal{J}}_{T}(y,v)\;:\;(y,v)\in\mathcal{A}\big{\\}}$
(100)

and
$\mathcal{A}\subset\mathcal{C}^{2}([0,T],\mathbb{R})\times\mathcal{C}^{1}([0,T],\mathbb{R}^{2})$
consists of all $(x,u)$ such that $x(0)=x^{*}(0)$,

$\displaystyle x(T)=\hat{x},$ (101)

$u(0)=u^{*}$ and for all $t\in[0,T]$,

$\displaystyle\frac{du}{dt}-A(u(t),x(t)).$ (102)

The trajectory that minimizes $\tilde{\mathcal{J}}_{T}$ indicates the most
likely means that a significant proportion of calcium channels open over the
time interval $[0,T]$ due to stochastic effects. The second order ODE for
$x(t)$ may be inferred from the previous section. The dynamics of the Lagrange
multiplier $\eta(t)$ is given by, and the boundary value at time $T$ is

$\displaystyle\eta^{1}(T)$ $\displaystyle=0$ (103) $\displaystyle\eta^{2}(T)$
$\displaystyle=0.$ (104)

Our system thus amounts to a boundary value ODE. It is solved using the
shooting method. The unknown derivatives at time $0$ (which must be solved for
using the shooting method) are $\dot{\eta}^{1}(0)$ , $\dot{\eta}^{2}(0)$,
$\dot{x}(0)$. There are three constraints at time $T$ that are used to solve
for the unknowns at time $0$: (101), (103) and (104).

###  3.3 Discussion

We have demonstrated that Large Deviation Theory can be used to compute the
most likely trajectory followed by $N$ stochastic channels / diffusing calcium
in making the spark to wave transition [19]. This model is relatively simple,
particularly insofar as the calcium concentration is spatially homogeneous. In
a subsequent work, we will employ a spatially-distributed model of the calcium
to compute the optimal trajectories.

##  4 Proofs

To prove Theorem 1.1, we apply the Contraction Principle [18] to the Large
Deviation Principle for the independent system. In other words, we are going
to perform a change of variable, which will transfer known results for the
Large Deviations of independent Poisson Processes [51, Chapter 6]. We first
recount the Large Deviations for independent Poisson Processes (recalling the
definition of $\\{y_{\alpha}(t)\\}_{\alpha\in\mathscr{R}}$ in (23)).

######  Theorem 4.1.

The sequence of probability laws of
$(y_{\alpha})_{\alpha\in\mathscr{R}}\in\tilde{\mathcal{D}}([0,\infty),\mathbb{R}^{+})^{M}$
satisfy a Large Deviation Principle with good rate function. That is, for
$\mathcal{O},\mathcal{A}\subset\tilde{\mathcal{D}}([0,\infty),\mathbb{R}^{+})^{M}$,
with $\mathcal{O}$ open and $\mathcal{A}$ closed,

$\displaystyle\underset{N\to\infty}{\overline{\lim}}N^{-1}\log\mathbb{P}\big{(}y\in\mathcal{A}\big{)}$
$\displaystyle\leq-\inf_{\alpha\in\mathcal{A}}\mathcal{I}(\alpha)$ (105)
$\displaystyle\underset{N\to\infty}{\underline{\lim}}N^{-1}\log\mathbb{P}\big{(}y\in\mathcal{O}\big{)}$
$\displaystyle\geq-\inf_{\alpha\in\mathcal{O}}\mathcal{I}(\alpha),$ (106)

and the rate function is

$\displaystyle\mathcal{I}(y)$
$\displaystyle=\left\\{\begin{array}[]{c}\infty\text{ in the case that
}y_{\alpha}\text{ is not absolutely continuous for some
}\alpha\in\mathscr{R}\\\
\sum_{\alpha\in\mathscr{R}}\int_{0}^{\infty}\ell\big{(}\dot{y}_{\alpha}(r)\big{)}dr\text{
otherwise, }\end{array}\right.$ (109) $\displaystyle\ell(r)$
$\displaystyle=r\log r-r+1.$ (110)

Furthermore the level sets of $\mathcal{I}$ are compact.

###### Proof.

The Large Deviations Principle for arbitrarily long finite time intervals is
proved in [51, Chapter 5], with corresponding rate function

$\mathcal{I}(y)=\sum_{\alpha\in\mathscr{R}}\int_{0}^{\infty}\ell\big{(}\dot{y}_{\alpha}(r)\big{)}dr.$
(111)

∎

The first result that we must prove is that the Large Deviations Principle
must hold for sets with reaction rates bounded away from zero.

######  Lemma 4.2.

Suppose that $\mathcal{O},\mathcal{A}\subseteq\Upsilon$, with $\mathcal{O}$
open and $\mathcal{A}$ closed, are such that for any $T>0$,

$\displaystyle\inf_{(z,x,u)\in\mathcal{A}\cup\mathcal{O}}\inf_{\alpha\in\mathscr{R}}\inf_{s\in[0,T]}\lambda_{\alpha}(x(s),u(s))>0.$
(112)

Then the sequence of probability laws of $\big{(}z,x,u\big{)}$ satisfy a Large
Deviation Principle on the space $\Upsilon$, i.e.

$\displaystyle\underset{N\to\infty}{\overline{\lim}}N^{-1}\log\mathbb{P}\big{(}(z,x,u)\in\mathcal{A}\big{)}$
$\displaystyle\leq-\inf_{\beta\in\mathcal{A}}\mathcal{J}(\beta)$ (113)
$\displaystyle\underset{N\to\infty}{\underline{\lim}}N^{-1}\log\mathbb{P}\big{(}(z,x,u)\in\mathcal{O}\big{)}$
$\displaystyle\geq-\inf_{\beta\in\mathcal{O}}\mathcal{J}(\beta).$ (114)

$\mathcal{J}$ is lower-semi-continuous and the level sets of $\mathcal{J}$ are
compact.

###### Proof.

We are going to prove the Large Deviation Principle using a contraction
principle. To this end, define a mapping
$\Psi:\Upsilon\to\tilde{\mathcal{D}}([0,\infty),\mathbb{R}^{+})^{M}$ as
follows. For any $(z,x,u)\in\Upsilon$, define

$\displaystyle\Lambda_{\alpha}(t)=\int_{0}^{t}\lambda_{\alpha}(x(s),u(s))ds.$
(115)

Let $\Lambda^{-1}_{\alpha}\in\tilde{\mathcal{D}}([0,\infty),\mathbb{R}^{+}$ be
such that

$\displaystyle\Lambda^{-1}_{\alpha}(t)=\inf\big{\\{}s\geq
0\;:\Lambda_{\alpha}(s)=t\big{\\}}.$ (116)

Note that, since $\inf_{s\in[0,t]}\lambda_{\alpha}(x(s),u(s))>0$, it must be
that $\Lambda^{-1}_{\alpha}$ is the function inverse of $\Lambda_{\alpha}$.
Write

$\displaystyle\tau_{\alpha}=\lim_{t\to\infty}\Lambda_{\alpha}(t),$ (117)

and note that $\tau_{\alpha}$ could be $\infty$. We define
$\Psi(z,x,u)=\big{(}w_{\alpha}\big{)}_{\alpha\in\mathscr{R}}$, where for any
$t<\tau_{\alpha}$,

$\displaystyle
w_{\alpha}(t)=z_{\alpha}\big{(}\Lambda^{-1}_{\alpha}(t)\big{)}.$ (118)

In the case that $\tau_{\alpha}<\infty$, for all $t\geq\tau_{\alpha}$, we
define $w_{\alpha}(t)=\lim_{s\to\tau_{\alpha}^{-}}w_{\alpha}(s)$. One observes
that over sets of the form (112), $\Psi$ is continuous. Furthermore,
$w_{\alpha}(t)=y_{\alpha}(t)$. We wish to apply the Inverse Contraction
Principle [18], (noting also Lemmas 4.1 and 4.4). In order that this theorem
applies, we must also prove that $\Psi$ is one-to-one. This is proved in Lemma
4.3.

∎

######  Lemma 4.3.

Over any set $\mathcal{A}_{\epsilon}$ of the form, for $\epsilon>0$,

$\displaystyle\mathcal{A}_{\epsilon}=\inf_{(z,x,u)\in\Upsilon}\inf_{\alpha\in\mathscr{R}}\inf_{s\in[0,T]}\lambda_{\alpha}(x(s),u(s))\geq\epsilon.$
(119)

$\Psi$ is one-to-one.

###### Proof.

Let $(z,x,u),(\hat{z},\hat{x},\hat{u})\in\Upsilon$. Write

$\displaystyle\Lambda_{\alpha}(t)$
$\displaystyle=\int_{0}^{t}\lambda_{\alpha}(x(s),u(s))ds$ (120)
$\displaystyle\hat{\Lambda}_{\alpha}(t)$
$\displaystyle=\int_{0}^{t}\lambda_{\alpha}(\hat{x}(s),\hat{u}(s))ds$ (121)
$\displaystyle w_{\alpha}(t)$
$\displaystyle=z_{\alpha}\big{(}\Lambda^{-1}_{\alpha}(t)\big{)}$ (122)
$\displaystyle\hat{w}_{\alpha}(t)$
$\displaystyle=\hat{z}_{\alpha}\big{(}\hat{\Lambda}^{-1}_{\alpha}(t)\big{)}.$
(123)

It is noted in Lemma that there is a constant $C_{t}>0$ such that

$\displaystyle\sup_{\alpha\in\mathscr{R}}\Lambda_{\alpha}(t)\leq C_{t}$ (124)

We are going to demonstrate that for any $t\geq 0$, there exists a constant
$c_{\epsilon,t}$ such that for all $s\leq t$,

$\displaystyle\tilde{d}^{\circ}_{s}(z_{\alpha},\hat{z}_{\alpha})\leq\tilde{d}^{\circ}_{C_{s}}(w_{\alpha},\hat{w}_{\alpha})+c_{\epsilon,t}\tilde{d}^{\circ}_{s}(z_{\alpha},\hat{z}_{\alpha}).$
(125)

An application of Gronwall’s Inequality to (125) implies that if
$w_{\alpha}=\hat{w}_{\alpha}$ for all $\alpha\in\mathscr{R}$, then necessarily
$z_{\alpha}=\hat{z}_{\alpha}$. This implies the Lemma.

With an aim of proving (125), it follows from Calculus that

$\displaystyle\Lambda^{-1}_{\alpha}(t)$
$\displaystyle=\int_{0}^{t}\lambda^{-1}_{\alpha}(x(s),u(s))ds$ (126)
$\displaystyle\hat{\Lambda}^{-1}_{\alpha}(t)$
$\displaystyle=\int_{0}^{t}\lambda^{-1}_{\alpha}(\hat{x}(s),\hat{u}(s))ds.$
(127)

Our assumptions on the set $\mathcal{A}_{\epsilon}$ imply that

$\displaystyle\Lambda^{-1}_{\alpha}(t)$ $\displaystyle\leq t\epsilon^{-1}$
(128) $\displaystyle\hat{\Lambda}^{-1}_{\alpha}(t)$ $\displaystyle\leq
t\epsilon^{-1}.$ (129)

Define

$\displaystyle\tilde{z}_{\alpha}(t)=w_{\alpha}\big{(}\hat{\Lambda}(t)\big{)}.$
(130)

The triangle inequality implies that

$\displaystyle\tilde{d}^{\circ}_{t}(z_{\alpha},\hat{z}_{\alpha})\leq\tilde{d}^{\circ}_{t}(z_{\alpha},\tilde{z}_{\alpha})+\tilde{d}^{\circ}_{t}(\tilde{z}_{\alpha},\hat{z}_{\alpha}).$
(131)

The definition of the Skorohod metric implies that

$\displaystyle\tilde{d}^{\circ}_{t}(z_{\alpha},\tilde{z}_{\alpha})\leq\sup_{s\leq
t}\bigg{|}\int_{0}^{s}\lambda_{\alpha}(x(r),u(r))/\lambda_{\alpha}(\hat{x}(r),\hat{u}(r))dr-s\bigg{|}.$
(132)

Our assumption that the reaction rates and functions are Lipschitz, together
with Lemma 1.2 implies that there is a constant $\bar{C}_{\epsilon,t}>0$ such
that

$\displaystyle\sup_{s\in[0,t]}\sup_{\alpha\in\mathscr{R}}\big{|}\lambda_{\alpha}(x(s),u(s))-\lambda_{\alpha}(\hat{x}(s),\hat{u}(s))\big{|}\leq$
$\displaystyle\bar{C}_{\epsilon,t}\sup_{s\in[0,t]}\sup_{\alpha\in\mathscr{R}}\big{|}z_{\alpha}(s)-\hat{z}_{\alpha}(s)\big{|}$
(133) $\displaystyle\leq$
$\displaystyle\bar{C}_{\epsilon,t}\sup_{\alpha\in\mathscr{R}}\tilde{d}^{\circ}_{t}(z_{\alpha},\hat{z}_{\alpha}).$
(134)

Since the reaction rates are (by assumption in the statement of the Lemma),
bounded from below, we thus obtain that there is a constant
$\grave{C}_{\epsilon,t}$ such that for all $s\leq\bar{t}$,

$\displaystyle\sup_{\alpha\in\mathscr{R}}\tilde{d}^{\circ}_{s}(z_{\alpha},\tilde{z}_{\alpha})\leq\grave{C}_{\epsilon,\bar{t}}\sup_{\alpha\in\mathscr{R}}\tilde{d}^{\circ}_{s}(z_{\alpha},\hat{z}_{\alpha}).$
(135)

The definition of the Skorohod Metric, together with the bounds in (128) also
implies that

$\displaystyle\tilde{d}^{\circ}_{t}(\tilde{z}_{\alpha},\hat{z}_{\alpha})\leq\tilde{d}^{\circ}_{C_{t}}(w_{\alpha},\hat{w}_{\alpha}).$
(136)

We have thus proved (125). ∎

The following exponential tightness result is a standard requirement for Large
Deviation Principles.

######  Lemma 4.4.

For any $L>0$, there exists a compact set $\mathcal{K}_{L}\subset\Upsilon$
such that

$\displaystyle\underset{N\to\infty}{\overline{\lim}}N^{-1}\log\mathbb{P}\big{(}(z,x,u)\notin\mathcal{K}_{L}\big{)}\leq-L.$
(137)

###### Proof.

We can equivalently formulate the system in terms of an empirical measure, and
the exponential tightness is an immediate consequence. In more detail, write

$\displaystyle\hat{\mu}^{N}=N^{-1}\sum_{j\in
I_{N}}\delta_{w^{j}}\in\mathcal{P}\big{(}\mathcal{D}\big{(}[0,T],\mathbb{Z}^{+}\big{)}^{|\mathscr{R}|}\big{)}.$
(138)

Here
$w^{j}:=\big{(}w^{j}_{\alpha}\big{)}_{\alpha\in\mathscr{R}}\subset\mathcal{D}\big{(}[0,T],\mathbb{Z}^{+}\big{)}^{|\mathscr{R}|}$
are inhomogeneous counting processes, i.e. they are such that

$\displaystyle\mathbb{P}\big{(}w^{j}_{\alpha}(t+\Delta)=w^{j}_{\alpha}(t)+1\;|\;\mathcal{F}_{t}\big{)}$
$\displaystyle\simeq\Delta\lambda_{\alpha}(x(t),u(t))+O(\Delta^{2})$ (139)
$\displaystyle\mathbb{P}\big{(}w^{j}_{\alpha}(t+\Delta)=w^{j}_{\alpha}(t)\;|\;\mathcal{F}_{t}\big{)}$
$\displaystyle\simeq 1-\Delta\lambda_{\alpha}(x(t),u(t))+O(\Delta^{2}).$ (140)

Here

$\displaystyle x(t)=$ $\displaystyle x(0)+N^{-1}\sum_{j\in
I_{N}}\sum_{\alpha\in\mathscr{R}}\xi_{\alpha}w^{j}_{\alpha}(t)$ (141)
$\displaystyle\frac{du}{dt}=$ $\displaystyle A(u(t),x(t))$ (142)

We then find that, substituting
$z_{\alpha}(t)=N^{-1}\sum_{j=1}^{N}w^{j}_{\alpha}(t)$, the above system has
the same probability law as the original system. Write
$\mathcal{D}_{*}\big{(}[0,T],\mathbb{Z}^{+}\big{)}^{|\mathscr{R}|}\big{)}\subseteq\mathcal{D}\big{(}[0,T],\mathbb{Z}^{+}\big{)}^{|\mathscr{R}|}$
to consist of all processes that are (i) equal to $0$ at time $0$, (ii) non-
decreasing.

More precisely, we see that there exists a continuous mapping
$\Psi:\mathcal{P}\big{(}\mathcal{D}_{*}\big{(}[0,T],\mathbb{Z}^{+}\big{)}^{|\mathscr{R}|}\big{)}$
such that $(z,x,u)=\Psi\big{(}\hat{\mu}^{N}\big{)}$, with unit probability.
For a positive number $L>0$, we are going to define a compact set
$\mathcal{K}_{L}\subseteq\mathcal{P}\big{(}\mathcal{D}_{*}\big{(}[0,T],\mathbb{Z}^{+}\big{)}^{|\mathscr{R}|}\big{)}$
such that

$\displaystyle\underset{N\to\infty}{\overline{\lim}}N^{-1}\log\mathbb{P}\big{(}\hat{\mu}^{N}\notin\mathcal{K}_{L}\big{)}\leq-L.$
(143)

This suffices for the lemma because the continuity of $\Psi$ then implies that
$\Psi(\mathcal{K}_{L})$ is compact.

For a positive integer $p\geq 1$, write
$\mathcal{U}_{p,T}\subseteq\mathcal{D}_{*}\big{(}[0,T],\mathbb{Z}^{+}\big{)}^{|\mathscr{R}|}$
to consist of all paths that are less than or equal to $p$ at time $T$. Its
easy to check that $\mathcal{U}_{p}$ is compact with respect to the Skorohod
Topology (there are at most $p$ ‘spike times’, and these spike times must be
in the compact time interval $[0,T]$). Write $Y^{p,j}_{\alpha}$ to be
independent counting processes of unit intensity. Thus the time-rescaled
representation of Poisson Processes [3] means that we can write

$\displaystyle
w^{j}_{\alpha}(t)=Y^{j}_{\alpha}\bigg{(}\int_{0}^{t}\lambda_{\alpha}(x(s),u(s))ds\bigg{)}$
(144)

For another number $b_{p}>0$, we find that

$\displaystyle\mathbb{P}\big{(}\hat{\mu}^{N}(\mathbf{w})(\mathcal{U}_{p,T})<1-b_{p}\big{)}\leq\mathbb{P}\big{(}\hat{\mu}^{N}(\mathbf{Y})\big{(}\mathcal{U}_{p,KT}\big{)}<1-b_{p}\big{)}.$
(145)

since $\lambda\leq K$ uniformly. Thanks to Chernoff’s Inequality, for a
constant $c>0$,

$\displaystyle\mathbb{P}\big{(}\hat{\mu}^{N}(\mathbf{Y})\big{(}\mathcal{U}_{p,KT}\big{)}<1-b_{p}\big{)}\leq$
$\displaystyle\mathbb{P}\bigg{(}\sup_{\alpha\in\mathscr{R}}N^{-1}\sum_{j=1}^{N}\chi\big{\\{}w^{j}_{\alpha}(KT)>p\big{\\}}\geq
b_{p}\bigg{)}$ (146) $\displaystyle\leq$
$\displaystyle\big{|}\mathscr{R}\big{|}\mathbb{E}\bigg{[}\exp\bigg{(}c\sum_{j=1}^{N}\chi\big{\\{}w^{j}_{\alpha}(KT)>p\big{\\}}-Ncb_{p}\bigg{)}\bigg{]}$
(147) $\displaystyle=$
$\displaystyle\big{|}\mathscr{R}\big{|}\bigg{\\{}1+\mathbb{P}\big{(}w^{j}_{\alpha}(KT)>p\big{)}\big{(}\exp(c)-1\big{)}\bigg{\\}}^{N}\exp\big{(}-Ncb_{p}\big{)}$
(148)

Through taking $p$ large enough, and $1\ll
c\ll-\log\mathbb{P}\big{(}w^{j}_{\alpha}(KT)>p\big{)}$, we find that

$\displaystyle\mathbb{P}\big{(}\hat{\mu}^{N}(\mathbf{w})(\mathcal{U}_{p,T})<1-b_{p}\big{)}\leq\exp\big{(}-pN\big{)}.$
(149)

Now define $\mathcal{K}_{L}$ to consist of all measures $\mu$ such that for
all $p\geq p_{L}$ (for an integer $p_{L}$ to be specified below),

$\displaystyle\mu(\mathcal{U}_{p,KT})\geq b_{p}.$ (150)

We thus find through a union of events bound that

$\displaystyle\mathbb{P}\bigg{(}\hat{\mu}^{N}\notin\mathcal{K}_{L}\bigg{)}\leq$
$\displaystyle\sum_{p=p_{L}}^{\infty}\mathbb{P}\bigg{(}\hat{\mu}^{N}(\mathbf{Y})\big{(}\mathcal{U}_{p,KT}\big{)}<1-b_{p}\bigg{)}$
(151) $\displaystyle\leq$
$\displaystyle\sum_{p=p_{L}}^{\infty}\exp\big{(}-pN\big{)}.$ (152)

Thus for large enough $p_{L}$, it must be that for all $N\geq 1$,

$\displaystyle
N^{-1}\log\mathbb{P}\bigg{(}\hat{\mu}^{N}\notin\mathcal{K}_{L}\bigg{)}\leq-L.$
(153)

∎

We finally note the proof of Lemma 1.2.

###### Proof.

It is clear from the definition of $x$ that there is a universal constant such
that

$\left\|x(s)-\tilde{x}(s)\right\|\leq
C\sup_{\alpha\in\mathscr{R}}|z_{\alpha}(s)-\tilde{z}_{\alpha}(s)|.$ (154)

It is demonstrated in Billingsley [6] that convergence in the Skorohod metric
implies convergence in the supremum norm, i.e as
$d^{\circ}_{\infty}(z,\tilde{z})\to 0$, it must be that

$\sup_{\alpha\in\mathscr{R}}\sup_{0\leq s\leq
t}\big{|}z_{\alpha}(s)-\tilde{z}_{\alpha}(s)\big{|}\to 0.$ (155)

The Lemma now follows from a standard application of Gronwall’s Inequality. ∎

## References

  * [1] Andrea Agazzi, Luisa Andreis, Robert I.A. Patterson, and D.R. Michiel Renger.  Large deviations for markov jump processes with uniformly diminishing rates.  Stochastic Processes and their Applications, 7 2022.
  * [2] Andrea Agazzi, Amir Dembo, and Jean Pierre Eckmann.  Large deviations theory for markov jump models of chemical reaction networks.  Annals of Applied Probability, 28:1821–1855, 2018.
  * [3] David F Anderson and Thomas G Kurtz.  Stochastic Analysis of Biochemical Systems.  Springer, 2015.
  * [4] Romain Azaïs and Aurélie Muller-Gueudin.  Optimal choice among a class of nonparametric estimators of the jump rate for piecewise-deterministic markov processes.  Electronic Journal of Statistics, 10:3648–3692, 2016.
  * [5] Howard C. Berg.  Random Walks in Biology.  Princeton University Press, 1993.
  * [6] Patrick Billingsley.  Convergence of Probability Measures.  1999\.
  * [7] Freddy Bouchet, Tobias Grafke, Tomás Tangarife, and Eric Vanden-Eijnden.  Large deviations in fast–slow systems.  Journal of Statistical Physics, 162:793–812, 2 2016.
  * [8] Paul Bressloff.  Stochastic Processes in Cell Biology. I.  Springer, 2 edition, 2021.
  * [9] Paul Bressloff and Olivier Faugeras.  On the hamiltonian structure of large deviations in stochastic hybrid systems.  Journal of Statistical Mechanics: Theory and Experiment, 2017.
  * [10] Paul C Bressloff.  Stochastic Processes in Cell Biology. II, volume II.  Springer, 2 edition, 2021.
  * [11] Paul C. Bressloff and Jay M. Newby.  Path integrals and large deviations in stochastic hybrid systems.  Physical Review E - Statistical, Nonlinear, and Soft Matter Physics, 89:1–15, 2014.
  * [12] P.C. Bressloff and J.N. Maclaurin.  Stochastic hybrid systems in cellular neuroscience.  Journal of Mathematical Neuroscience, 8, 2018.
  * [13] Amarjit Budhiraja and Paul Dupuis.  Analysis and Approximation of Rare Events, volume 94.  Springer, 2019.
  * [14] Amarjit Budhiraja and Eric Friedlander.  Diffusion approximations for controlled weakly interacting large finite state systems with simultaneous jumps.  Source: The Annals of Applied Probability, 28:204–249, 2018.
  * [15] Pengxing Cao, Graham Donovan, Martin Falcke, and James Sneyd.  A stochastic model of calcium puffs based on single-channel data.  Biophysical Journal, 105:1133–1142, 2013.
  * [16] S. Coombes, R. Hinch, and Y. Timofeeva.  Receptors, sparks and waves in a fire-diffuse-fire framework for calcium release.  Progress in Biophysics and Molecular Biology, 85:197–216, 2004\.
  * [17] S. Coombes and Y. Timofeeva.  Sparks and waves in a stochastic fire-diffuse-fire model of ca2+ release.  Physical Review E, 68:1–8, 2003.
  * [18] Amir Dembo and Ofer Zeitouni.  Large Deviations Techniques and Applications 2nd Edition.  Springer, 1998.
  * [19] Geneviève Dupont, Martin Falcke, Vivien Kirk, and James Sneyd.  Models of calcium signalling, volume 43.  Springer, 2016.
  * [20] Paul Dupuis, Kavita Ramanan, and Wei Wu.  Large deviation principle for finite-state mean field interacting particle systems.  Arxiv Preprint, 1 2016.
  * [21] Weinan E and Eric Vanden-Eijnden.  Transition-path theory and path-finding algorithms for the study of rare events.  Annual Review of Physical Chemistry, 61:391–420, 5 2010.
  * [22] A. Faggionato, D. Gabrielli, and M. Ribezzi Crivellari.  Non-equilibrium thermodynamics of piecewise deterministic markov processes.  Journal of Statistical Physics, 137:259–304, 2009.
  * [23] Victor Nicolai Friedhoff, Lukas Ramlow, Benjamin Lindner, and Martin Falcke.  Models of stochastic ca 2 + spiking: Established approaches and inspirations from models of neuronal spikes.  European Physical Journal: Special Topics, 230:2911–2928, 10 2021\.
  * [24] Tobias Grafke and Eric Vanden-Eijnden.  Non-equilibrium transitions in multiscale systems with a bifurcating slow manifold.  Journal of Statistical Mechanics: Theory and Experiment, 2017, 2017\.
  * [25] Martin Hasler, Vladimir Belykh, and Igor Belykh.  Dynamics of stochastically blinking systems. part i: Finite time properties.  SIAM Journal on Applied Dynamical Systems, 12:1007–1030, 2013.
  * [26] Qi He and G Yin.  Large deviations for multi-scale markovian switching systems with a small diffusion *.  Asymptotic Analysis, 2014.
  * [27] Alexandru Hening, Dang H. Nguyen, Nhu Nguyen, and Harrison Watts.  Random switching in an ecosystem with two prey and one predator.  SIAM Journal on Mathematical Analysis, 55:347–366, 2023.
  * [28] Robert Hinch and S. Jon Chapman.  Exponentially slow transitions on a markov chain: The frequency of calcium sparks.  European Journal of Applied Mathematics, 16:427–446, 2005.
  * [29] James Keener and James Sneyd.  Mathematical Physiology I: Cellular Physiology 2nd Edition.  Springer, 2008.
  * [30] James P. Keener.  Stochastic calcium oscillations.  Mathematical Medicine and Biology, 23:1–25, 2006.
  * [31] James P. Keener and Jay M. Newby.  Perturbation analysis of spontaneous action potential initiation by stochastic ion channels.  Physical Review E - Statistical, Nonlinear, and Soft Matter Physics, 84, 7 2011.
  * [32] Joel Keizer and Gregory D. Smith.  Spark-to-wave transition: Saltatory transmission of calcium waves in cardiac myocytes.  Biophysical Chemistry, 72:87–100, 1998.
  * [33] Joel Keizer, Gregory D. Smith, Silvina Ponce-Dawson, and John E. Pearson.  Saltatory propagation of ca2+ waves by ca2+ sparks.  Biophysical Journal, 75:595–600, 1998.
  * [34] Yuri Kifer.  Large deviations and adiabatic transitions for dynamical systems and Markov processes in fully coupled averaging.  2009\.
  * [35] Hyunjoong Kim, Yoichiro Mori, and Joshua B. Plotkin.  Optimality of intercellular signaling: direct transport versus diffusion.  Physical Review E, 4 2022.
  * [36] Richard C. Kraaij and Mikola C. Schlottke.  A large deviation principle for markovian slow-fast systems.  11 2020.
  * [37] Rohini Kumar and Lea Popovic.  Large deviations for multi-scale jump-diffusion processes.  Stochastic Processes and their Applications, 127:1297–1320, 2017\.
  * [38] Thomas G. Kurtz.  Limit theorems for sequences of jump markov processes approximating ordinary differential processes.  Journal of Applied Probability, 8:344–356, 1971.
  * [39] Sean D. Lawley, Janet A. Best, and Michael C. Reed.  Neurotransmitter concentrations in the presence of neural switching in one dimension.  Discrete and Continuous Dynamical Systems - Series B, 21:2255–2273, 2016.
  * [40] James MacLaurin and Jay M. Newby.  Extreme first passage times for populations of identical rare events.  Arxiv Preprint, 9 2023.
  * [41] Gregor Moenke, Martin Falcke, and Keven Thurley.  Hierarchic stochastic modelling applied to intracellular ca2+ signals.  PLoS ONE, 7, 12 2012.
  * [42] Etienne Pardoux and Brice Samegni-Kepgnou.  Large deviation principle for epidemic models.  Source: Journal of Applied Probability, 54:905–920, 2017.
  * [43] Robert I.A. Patterson and D. R.Michiel Renger.  Large deviations of jump process fluxes.  Mathematical Physics Analysis and Geometry, 22, 9 2019.
  * [44] Lukas Ramlow, Martin Falcke, and Benjamin Lindner.  An integrate-and-fire approach to ca2+ signaling. part i: Renewal model.  Biophysical Journal, 122:713–736, 2 2023.
  * [45] Lukas Ramlow, Martin Falcke, and Benjamin Lindner.  An integrate-and-fire approach to ca2+ signaling. part ii: Cumulative refractoriness.  Biophysical Journal, 122:4710–4729, 12 2023.
  * [46] Martin G. Riedler, Michèle Thieullen, and Gilles Wainrib.  Limit theorems for infinite-dimensional piecewise deterministic markov processes. applications to stochastic excitable membrane models.  Electronic Journal of Probability, 17, 2012.
  * [47] Ryszard Rudnicki and Marta Tyran-Kaminska.  Piecewise Deterministic Processes in Biological Models.  Springer Briefs in Applied Sciences and Technology, 2017.
  * [48] Sten Rüdiger.  Stochastic models of intracellular calcium signals.  Physics Reports, 534:39–87, 2014.
  * [49] Ivo Siekmann, Larry E. Wagner, David Yule, Edmund J. Crampin, and James Sneyd.  A kinetic model for type i and ii ip3r accounting for mode changes.  Biophysical Journal, 103:658–668, 2012.
  * [50] Alexander Skupin, Helmut Kettenmann, Ulrike Winkler, Maria Wartenberg, Heinrich Sauer, Stephen C. Tovey, Colin W. Taylor, and Martin Falcke.  How does intracellular ca2+ oscillate: By chance or by the clock?  Biophysical Journal, 94:2404–2411, 3 2008.
  * [51] Alan Weiss and Adam Schwartz.  Large Deviations for Performance Analysis: queues, communication and computing. Routledge, 1995.