Paul C Bressloff
email: bressloff@math.utah.edu

Olivier Faugeras
email: olivier.faugeras@inria.fr

On the Hamiltonian structure of large deviations in stochastic hybrid systems

Keywords: Hamiltonian, large deviations, WKB, stochastic hybrid systems

We present a new derivation of the classical action underlying a large deviation principle (LDP) for a stochastic hybrid system, which couples a piecewise deterministic dynamical system in R d with a time-homogeneous Markov chain on some discrete space Γ. We assume that the Markov chain on Γ is ergodic, and that the discrete dynamics is much faster than the piecewise deterministic dynamics (separation of time-scales). Using the Perron-Frobenius theorem and the calculus-of-variations, we show that the resulting Hamiltonian is given by the Perron eigenvalue of a |Γ|-dimensional linear equation. The corresponding linear operator depends on the transition rates of the Markov chain and the nonlinear functions of the piecewise deterministic system. We compare the Hamiltonian to one derived using WKB methods, and show that the latter is a reduction of the former. We also indicate how the analysis can be extended to a multiscale stochastic process, in which the continuous dynamics is described by a piecewise stochastic differential equations (SDE). Finally, we illustrate the theory by considering applications to conductance-based models of membrane voltage fluctuations in the presence of stochastic ion channels.

Introduction

There are a growing number of problems in biology that involve the coupling between a piecewise deterministic dynamical system in R d and a time-homogeneous Markov chain on some discrete space Γ [START_REF] Bressloff | Path-integrals and large-deviations in stochastic hybrid systems[END_REF], resulting in an a type of stochastic hybrid system (SHS) known as a piecewise deterministic Markov process (PDMP) [START_REF] Davis | Piecewise-deterministic Markov processes: A general class of non-diffusion stochastic models[END_REF] ‡. One important example at the single-cell level is given by membrane voltage fluctuations in neurons due to the stochastic opening and closing of ion channels [START_REF] Fox | Emergent collective behavior in large numbers of globally coupled independent stochastic ion channels[END_REF][START_REF] Chow | Spontaneous action potentials due to channel fluctuations[END_REF][START_REF] Keener | Perturbation analysis of spontaneous action potential initiation by stochastic ion channels[END_REF][START_REF] Goldwyn | The what and where of adding channel noise to the Hodgkin-Huxley equations[END_REF][START_REF] Buckwar | An exact stochastic hybrid model of excitable membranes including spatio-temporal evolution[END_REF][START_REF] Pakdaman | Fluid limit theorems for stochastic hybrid systems with application to neuron models[END_REF][START_REF] Wainrib | Reduction of stochastic conductance-based neuron models with time-scales separation[END_REF][START_REF] Newby | Breakdown of fast-slow analysis in an excitable system with channel noise[END_REF][START_REF] Bressloff | Stochastic hybrid model of spontaneous dendritic NMDA spikes[END_REF][START_REF] Anderson | Stochastic representations of ion channel kinetics and exact stochastic simulation of neuronal dynamics[END_REF]. Here the discrete states of the ion channels evolve according to a continuous-time Markov process with voltage-dependent transition rates and, in-between discrete jumps in the ion channel states, the membrane voltage evolves according to a deterministic equation that depends on the current state of the ion channels. In the limit that the number of ion channels goes to infinity, one can apply the law of large numbers and recover classical Hodgkin-Huxley type equations. However, finite-size effects can result in the noise-induced spontaneous firing of a neuron due to channel fluctuations. Another major example is a gene regulatory network, where the continuous variable is the concentration of a protein product and the discrete variables represents the activation state of the gene [25,[START_REF] Assaf | Spectral theory of metastability and extinction in birth-death systems[END_REF][START_REF] Newby | Isolating intrinsic noise sources in a stochastic genetic switch[END_REF][START_REF] Newby | Metastable behavior in Markov processes with internal states[END_REF][START_REF] Newby | Bistable switching asymptotics for the self regulating gene[END_REF]. A third example concerns a stochastic formulation of synaptically-coupled neural networks that has a mathematical structure analogous to gene networks [START_REF] Bressloff | Metastability in a stochastic neural network modeled as a velocity jump Markov process[END_REF][START_REF] Bressloff | Path-integrals and large-deviations in stochastic hybrid systems[END_REF].

In many of the above examples, one finds that the transition rates between the discrete states n ∈ Γ are much faster than the relaxation rates of the piecewise deterministic dynamics for x ∈ R d . Thus there is a separation of time scales between the discrete and continuous processes, so that if t is the characteristic time-scale of the Markov chain then t is the characteristic time-scale of the relaxation dynamics for some small positive parameter . Assuming that the Markov chain is ergodic, in the limit → 0 one obtains a deterministic dynamical system in which one averages the piecewise dynamics with respect to the corresponding unique stationary measure. This then raises the important problem of characterizing how the law of the underlying stochastic process approaches this deterministic limit in the case of weak noise, 0 < 1. A rigorous mathematical approach to addressing the above issue is large deviation theory, which has been developed extensively within the context of stochastic differential equations (SDEs) [START_REF] Freidlin | Random perturbations of dynamical systems[END_REF][START_REF] Dembo | Large Deviations Techniques and Applications[END_REF][START_REF] Touchette | The large deviation approach to statistical mechanics[END_REF]. In particular, consider some random dynamical system in R d for which there exists a well defined probability density functional or law P [x] over the different sample trajectories {x(t)} T 0 in a given time interval [0, T ]. Here is a small parameter that characterizes the noise level, with x(t) given by the solution x * (t) of some ODE ẋ = F (x) in the limit → 0. A large deviation principle (LDP) for the random paths of the SDE over some time interval [0, T ] is

P [x] ∼ e -J T [x]/ , → 0,
where J T [x] is known as the rate function and J T [x * ] = 0. In the case of SDEs, the rate function can be interpreted as a classical action with corresponding Lagrangian ‡ In this paper we use the term SHS mainly within the restricted sense of a PDMP. However, we do briefly consider a more general class of SHS in section 3, in which the piecewise deterministic continuous process is replaced by a piecewise stochastic continuous process.

L [START_REF] Freidlin | Random perturbations of dynamical systems[END_REF],

J T [x] = T 0 L(x, ẋ)dt.
Such a Lagrangian formulation is more amenable to explicit calculations. In particular, it can be used to solve various first passage time problems associated with the escape from a fixed point attractor of the underlying deterministic system in the weak noise limit. This involves finding the most probable paths of escape, which minimize the action with respect to the set of all trajectories emanating from the fixed point.

Evaluating the action along a most probable path from the fixed point to another point x generates a corresponding quasipotential Φ(x). From classical variational analysis, it can be shown that the quasipotential satisfies a Hamilton-Jacobi equation H(x, ∂ x Φ) = 0, where H is the Hamiltonian obtained from L [START_REF] Freidlin | Random perturbations of dynamical systems[END_REF] via a Fenchel-Legendre transformation:

H(x, p) = sup y {py -L(x, y)}.

The optimal paths of escape correspond to solutions of Hamilton's equations on the zero energy surface (H = 0). Interestingly, the same Hamilton-Jacobi equation is obtained by considering a Wentzel-Kramers-Brillouin (WKB) approximation of the stationary state of the continuous process x(t) in the weak-noise limit [START_REF] Matkowsky | The exit problem for randomly perturbed dynamical systems[END_REF][START_REF] Naeh | A direct approach to the exit problem[END_REF][START_REF] Maier | Limiting exit location distribution in the stochastic exit problem[END_REF]. Analogous connections between large deviation theory and WKB methods have also been established for continuous Markov chains [START_REF] Hanggi | Bistable systems: master equation versus Fokker-Planck modeling[END_REF][START_REF] Knessl | An asymptotic theory of large deviations for Markov jump processes[END_REF][START_REF] Dykman | Large fluctuations and optimal paths in chemical kinetics[END_REF][START_REF] Vellela | Stochastic dynamics and non-equilibrium thermodynamics of a bistable chemical system: the Schlögl model revisited[END_REF]. More recently, rigorous large deviation theory has been applied to PDMPs [START_REF] Faggionato | Averaging and large deviation principles for fully-coupled piecewise deterministic Markov processes and applications to molecular motors[END_REF][START_REF] Faggionato | Non-equilibrium thermodynamics of piecewise deterministic Markov Processes[END_REF][START_REF] Kifer | Large deviations and adiabatic transitions for dynamical systems and Markov processes in fully coupled averaging[END_REF]. Independently of these developments in large deviation theory, a variety of techniques in applied mathematics and mathematical physics have been used to solve first passage time problems in biological applications of stochastic hybrid systems. These include WKB approximations and matched asymptotics [START_REF] Keener | Perturbation analysis of spontaneous action potential initiation by stochastic ion channels[END_REF][START_REF] Newby | An asymptotic analysis of the spatially inhomogeneous velocity-jump process[END_REF][START_REF] Newby | Isolating intrinsic noise sources in a stochastic genetic switch[END_REF][START_REF] Newby | Breakdown of fast-slow analysis in an excitable system with channel noise[END_REF][START_REF] Bressloff | Stochastic hybrid model of spontaneous dendritic NMDA spikes[END_REF], and path-integrals [START_REF] Bressloff | Path-integrals and large-deviations in stochastic hybrid systems[END_REF]. Although such approaches are less rigorous than large deviation theory, they are more amenable to explicit calculations. In particular, they allow one to calculate the prefactor in Arrhenius-like expressions for mean first passage times, rather than just the leading order exponential behavior governed by the quasipotential. A major aim aim of this paper is to make explicit the connection between large deviation theory and more applied approaches to stochastic hybrid systems, by highlighting the common underlying Hamiltonian structure. Consistent with this aim, we present a new derivation of the Hamiltonian formulation of an LDP for stochastic hybrid systems. We take as our starting point the LDP due to Faggionato et al. [START_REF] Faggionato | Averaging and large deviation principles for fully-coupled piecewise deterministic Markov processes and applications to molecular motors[END_REF][START_REF] Faggionato | Non-equilibrium thermodynamics of piecewise deterministic Markov Processes[END_REF], which is equivalent to the more abstract LDP of Kifer [START_REF] Kifer | Large deviations and adiabatic transitions for dynamical systems and Markov processes in fully coupled averaging[END_REF]. Using the Perron-Frobenius theorem and the calculus-of-variations, we evaluate the LDP rate function in terms of a classical action, whose equations of motion are given by a Hamiltonian dynamical system.

The structure of the paper is as follows. In section 2 we define a stochastic hybrid system and specify our various mathematical assumptions. We then present our detailed derivation of the LDP Hamiltonian in section 3. We thus establish that the Hamiltonian is given by the Perron eigenvalue of a matrix equation that depends on the transition rates of the Markov chain and the nonlinear functions of the piecewise deterministic system. This result is consistent with alternative formulations based either on more abstract probability theory [START_REF] Kifer | Large deviations and adiabatic transitions for dynamical systems and Markov processes in fully coupled averaging[END_REF] or formal path-integral methods [START_REF] Bressloff | Path-integrals and large-deviations in stochastic hybrid systems[END_REF][START_REF] Bressloff | Path-integral methods for analyzing the effects of fluctuations in stochastic hybrid neural networks[END_REF]. We also indicate how the analysis can be extended to a multi-scale stochastic process, in which the slow dynamics is given by a piecewise SDE. In section 4 we construct a WKB approximation of solutions to a differential Chapman-Kolmogorov equation, which describes the evolution of the probability density function for a stochastic hybrid system. In particular, we show that the resulting quasi-potential satisfies a Hamilton-Jacobi equation, whose associated Hamiltonian is a reduced version of the full LDP Hamiltonian. Finally, we illustrate our analysis of LDPs for stochastic hybrid systems by considered some conductance-based models of membrane voltage fluctuations in the presence of stochastic ion channels (section 5).

Stochastic hybrid systems

Consider a one-dimensional SHS corresponding to a PDMP with continuous variable [START_REF] Davis | Piecewise-deterministic Markov processes: A general class of non-diffusion stochastic models[END_REF][START_REF] Kifer | Large deviations and adiabatic transitions for dynamical systems and Markov processes in fully coupled averaging[END_REF]. (Note that one could extend the analysis to higher-dimensions, x ∈ R d . In this case Ω is taken to be a connected, bounded domain with a regular boundary ∂Ω, see section 3.4. It is also possible to have more than one discrete variable, but one can always relabel the discrete states so that they are effectively indexed by a single integer.) When the internal state is n, the system evolves according to the ordinary differential equation (ODE)

x ∈ Ω ⊂ R and a discrete variable n ∈ Γ ≡ {0, • • • , N -1}
ẋ = F n (x), (2.1)
where the vector field F n : R → R is a continuous function, locally Lipschitz. That is, given a compact subset K of Ω, there exists a positive constant K n such that

|F n (x) -F n (y)| ≤ K n |x -y|, ∀x, y ∈ Ω (2.2)
for some constant K n . We assume that the dynamics of x is confined to the domain Ω so that we have existence and uniqueness of a trajectory for each n. One final mild constraint is that the vector field does not have identical components anywhere in Ω, that is, for any x ∈ Ω, there exists at least one pair (n, m) ∈ Γ, n = m for which F m (x) = F n (x). For fixed x, the discrete stochastic variable evolves according to a homogeneous, continuous-time Markov chain with transition matrix W(x). We make the further assumption that the chain is irreducible for all x ∈ Σ, that is, for fixed x there is a non-zero probability of transitioning, possibly in more than one step, from any state to any other state of the Markov chain. This implies the existence of a unique invariant probability distribution on Γ for fixed x ∈ Σ, denoted by ρ(x, n), such that

n∈Γ ρ(x, n)W nm (x) = ρ(x, m), ∀m ∈ Γ, ∀x ∈ Ω. (2.
3)

The existence of the unique invariant measure is a consequence of the well known Perron-Frobenius Theorem §. The above stochastic model defines a one-dimensional PDMP. It is also possible to consider generalizations of the continuous process, in which the ODE (2.1) is replaced by a stochastic differential equation (SDE) or even a partial differential equation (PDE). In order to allow for such possibilities we will refer to all of these processes as examples of a stochastic hybrid system.

Let us decompose the transition matrix of the Markov chain as W nm (x) = P nm (x)λ m (x), § A finite-dimensional, real square matrix with positive entries has a unique largest real eigenvalue (the Perron eigenvalue) and the corresponding eigenvector has strictly positive components [START_REF] Grimmett | Probability and random processes 3rd[END_REF].

with n =m P nm (x) = 1 for all x. Hence λ m (x) determines the jump times from the state m whereas P nm (x) determines the probability distribution that when it jumps the new state is n for n = m. The hybrid evolution of the system with respect to x(t) and n(t) is then described as follows. Suppose the system starts at time zero in the state (x 0 , n 0 ). Call x 0 (t) the solution of (2.1) with n = n 0 such that x 0 (0) = x 0 . Let θ 1 be the random variable such that

P(θ 1 < t) = 1 -exp - t 0 λ n0 (x 0 (t ))dt .
for some ε, 0 < ε 1. Then in the random time interval s ∈ [0, θ 1 ) the state of the system is (x 0 (s), n 0 ). We draw a value of θ 1 from P(θ 1 < t), choose an internal state n 1 ∈ Γ with probability P n1n0 (x 0 (θ 1 )), and call x 1 (t) the solution of the following Cauchy problem on [θ 1 , ∞):

ẋ1 (t) = F n1 (x 1 (t)), t ≥ θ 1 x 1 (θ 1 ) = x 0 (θ 1 )
Iterating this procedure, we construct a sequence of increasing jumping times (θ k ) k≥0 (setting θ 0 = 0) and a corresponding sequence of internal states (n k ) k≥0 . The evolution (x(t), n(t)) is then defined as

(x(t), n(t)) = (x k (t), n k ) if θ k ≤ t < θ k+1 .
(2.4)

Note that the path x(t) is continuous and piecewise C 1 . In order to have a well-defined dynamics on [0, T ], it is necessary that almost surely the system makes a finite number of jumps in the time interval [0, T ]. This is guaranteed in our case.

Given the above iterative definition of the stochastic hybrid process, let X(t) and N (t) denote the stochastic continuous and discrete variables, respectively, at time t, t > 0, given the initial conditions X(0) = x 0 , N (0) = n 0 . Although the evolution of the continuous variable X(t) or the discrete variable N (t) is non-Markovian, it can be proven that the joint evolution (X(t), N (t)) is a strong Markov process [START_REF] Davis | Piecewise-deterministic Markov processes: A general class of non-diffusion stochastic models[END_REF]. Introduce the probability density ρ(x, n, t|x 0 , n 0 , 0) with P{X(t) ∈ (x, x + dx), N (t) = n|x 0 , n 0 ) = ρ(x, n, t|x 0 , n 0 , 0)dx. It follows that ρ evolves according to the forward differential Chapman-Kolmogorov (CK) equation [START_REF] Gardiner | Handbook of stochastic methods[END_REF][START_REF] Bressloff | Path-integrals and large-deviations in stochastic hybrid systems[END_REF] ∂ρ ∂t = -Lρ,

with the adjoint generator L (dropping the explicit dependence on initial conditions) defined according to

Lρ(x, n, t) = ∂F n (x)ρ(x, n, t) ∂x - 1 ε m∈Γ A nm (x)ρ(x, m, t), (2.6)
with

A nm = W nm -δ n,m , (2.7)
such that m∈Γ A nm (x) = 0 ∀m ∈ Γ and n∈Γ ρ(x, n)A nm (x) = 0. The first term on the right-hand side of equation (2.6) represents the probability flow associated by the piecewise deterministic dynamics for a given n, whereas the second term represents jumps in the discrete state n. Note that we have rescaled the matrix A by introducing the dimensionless parameter ε, ε > 0. This is motivated by the observation that in many biological applications there a separation of time-scales between the relaxation time for the dynamics of the continuous variable x and the rate of switching between the different discrete states n [START_REF] Bressloff | Path-integrals and large-deviations in stochastic hybrid systems[END_REF].

Let us now define the averaged vector field F : R → R by

F (x) = n∈Γ ρ(x, n)F n (x)
It can be shown [START_REF] Faggionato | Averaging and large deviation principles for fully-coupled piecewise deterministic Markov processes and applications to molecular motors[END_REF] that, given the assumptions on the matrix W, the functions ρ(x, n) on Ω belong to C 1 (R) for all n ∈ Γ and that this implies that F (x) is locally Lipschitz. Hence, for all t ∈ [0, T ], the Cauchy problem

ẋ(t) = F (x(t)) x(0) = x 0 (2.8)
has a unique solution for all n ∈ Γ. Intuitively speaking, one would expect the stochastic hybrid system (2.1) to reduce to the deterministic dynamical system (2.8) in the limit → 0. That is, for sufficiently small , the Markov chain undergoes many jumps over a small time interval ∆t during which ∆x ≈ 0, and thus the relative frequency of each discrete state n is approximately ρ(x, n). This can be made precise in terms of a Law of Large Numbers for stochastic hybrid systems proven in [START_REF] Faggionato | Averaging and large deviation principles for fully-coupled piecewise deterministic Markov processes and applications to molecular motors[END_REF].

Classical Hamiltonian from a large deviation principle

We now turn to the major part of the paper, namely, a new derivation of a classical Hamiltonian action corresponding to the abstract rate function of an LDP for stochastic hybrid systems introduced by Faggionato et al. [START_REF] Faggionato | Averaging and large deviation principles for fully-coupled piecewise deterministic Markov processes and applications to molecular motors[END_REF][START_REF] Faggionato | Non-equilibrium thermodynamics of piecewise deterministic Markov Processes[END_REF]. The advantage of our approach is that it avoids many of the technical difficulties associated with probabilistic approaches to large deviation theory [START_REF] Kifer | Large deviations and adiabatic transitions for dynamical systems and Markov processes in fully coupled averaging[END_REF].

In order to write down the LDP of Faggionato et al [START_REF] Faggionato | Averaging and large deviation principles for fully-coupled piecewise deterministic Markov processes and applications to molecular motors[END_REF][START_REF] Faggionato | Non-equilibrium thermodynamics of piecewise deterministic Markov Processes[END_REF], it is first necessary to introduce some notation. Let M + ([0, T ]) denote the space of probability measures on the interval [0, T ] and take the K-dimensional vector {ψ(t)} t∈[0,T ] to be an element of the product space M + ([0, T ]) Γ . In other words, for each t ∈ [0, T ],

ψ(t) = (ψ 1 (t), . . . , ψ K (t)) such that ψ n (t) ≥ 0, n∈Γ ψ n (t) = 1.
A particular realization of the stochastic process, {(x(t), n(t))} t∈[0,T ] , then lies in the product space

C([0, T ]) × M + ([0, T ]) Γ with ψ n (t) = 1 {n(t)=n} ≡ 1, if n(t) = n, 0, if n(t) = n (3.1)
and

x(t) = x 0 + T 0 n∈Γ ψ n (s)F n (x(s))ds. (3.2) Let Y x0 denote the subspace of C([0, T ]) × M + ([0, T ]) Γ for which equation (3.2) holds but ψ is now a general element of M + ([0, T ]) Γ .
Such a space contains both the set of trajectories of the stochastic hybrid system with ψ n (t) given by equation (3.1) and n(t) evolving according to the Markov chain, and the solution x * (t) of the averaged system (2.8) for which

ψ n = ψ * n with ψ * n (t) = ρ(n, x * (t)). It can be proven that Y x0 is a compact subspace of C([0, T ]) × M + ([0, T ]) Γ with topology defined by the metric [15] d({(x(t), n(t))} t∈[0,T ] , {( x(t), n(t))} t∈[0,T ] ) = sup t∈[0,T ] |x(t) -x(t)| + n∈Γ sup 0≤t≤T t 0 [ψ n (s) -ψ n (s)]ds .
Finally, we take P x0 to be the probability density functional or law of the set of trajectories {x(t)} t∈[0,T ] ∈ C([0, T ], Ω).

The following large deviation principle then holds [START_REF] Faggionato | Averaging and large deviation principles for fully-coupled piecewise deterministic Markov processes and applications to molecular motors[END_REF][START_REF] Faggionato | Non-equilibrium thermodynamics of piecewise deterministic Markov Processes[END_REF]: Given an element

{x(t)} t∈[0,T ] ∈ C([0, T ], Ω), P x0 {x(t)} t∈[0,T ] ∼ e -J T ({x(t)} t∈[0,T ] )/ . (3.3)
where the rate function

J [0,T ] : C([0, T ], Ω) → [0, ∞) is given by J T ({x(t)} t∈[0,T ] ) = inf {ψ(t)} t∈[0,T ] :{x(t),ψ(t)} t∈[0,T ] ∈Yx 0 T 0 j(x(t), ψ(t)) dt (3.4) and j(x, ψ) = sup ψ n W nm (x) 1 - z m z n . (3.5)
Here the symbol ∼ means asymptotic logarithmic equivalence in the limit → 0. We summarize a few useful properties of j(x, ψ) defined by equation (3.5). Since the double sum in equation (3.5) excludes diagonal terms, we introduce the set Γ ∆ = Γ × Γ\∆ where ∆ is the diagonal of Γ × Γ. Let c nm = ψ n W nm (x). Suppose that ψ is a strictly positive measure, ψ n > 0 for all n ∈ Γ. It then follows from the properties of the transition matrix W that the mapping c : Γ ∆ → [0, ∞) is irreducible in the sense that, for all n = m ∈ Γ, there exists a finite sequence

n 1 , n 2 , • • • , n k such that n 1 = n, n k = m and c nini+1 > 0 for i = 1, • • • , k -1. Define the mapping J : [0, ∞) Γ∆ → R as J (c) = sup z∈[0,∞) Γ Ĵ (c, z), Ĵ (c, z) = (n,m)∈Γ∆ c nn (1 -z m /z n ) (3.6)
The following lemma is proven in [START_REF] Faggionato | Averaging and large deviation principles for fully-coupled piecewise deterministic Markov processes and applications to molecular motors[END_REF][START_REF] Kifer | Large deviations and adiabatic transitions for dynamical systems and Markov processes in fully coupled averaging[END_REF] and follows largely from material found in [START_REF] Rockafellar | Convex analysis[END_REF].

Lemma 1: The function J is convex and continuous and takes its values in [0, ∞). Moreover, for each irreducible c, the supremum on [0, ∞) Γ of the function Ĵ (c, •) is a maximum and is the unique solution of the set of equations

m∈Γ c nm z m z n = m∈Γ c n n z n z m n ∈ Γ (3.7)
under the normalization n∈Γ z n = 1.

A key idea behind the above LDP is that a slow dynamical process coupled to the fast Markov chain on Γ rapidly samples the different discrete states of Γ according to some non-negative measure ψ. In the limit → 0, one has ψ → ρ, where ψ is the ergodic measure of the Markov chain. On the other hand, for small but nonzero , ψ is itself distributed according to a LDP, whereby one averages the different functions F n (x) over the measure ψ to determine the dynamics of the slow system. In most biological applications, one is not interested in the internal discrete state of the system, that is, one only observes the statistical behavior of the continuous variable x(t). For example, x(t) could represent the membrane voltage of a single neuron [START_REF] Newby | Breakdown of fast-slow analysis in an excitable system with channel noise[END_REF] or the synaptic current in a population of neurons [START_REF] Bressloff | Metastability in a stochastic neural network modeled as a velocity jump Markov process[END_REF]. Faggionato et al. [START_REF] Faggionato | Averaging and large deviation principles for fully-coupled piecewise deterministic Markov processes and applications to molecular motors[END_REF][START_REF] Faggionato | Non-equilibrium thermodynamics of piecewise deterministic Markov Processes[END_REF] explicitly calculated the rate function (3.4) for a restricted class of stochastic hybrid systems, whose stationary density is exactly solvable. One major constraint on this class of model is that the vector field of the piecewise deterministic system has non-vanishing components within a given confinement domain. However, this constraint does not hold for biological systems such as ion channels [START_REF] Keener | Perturbation analysis of spontaneous action potential initiation by stochastic ion channels[END_REF][START_REF] Newby | Breakdown of fast-slow analysis in an excitable system with channel noise[END_REF][START_REF] Bressloff | Stochastic hybrid model of spontaneous dendritic NMDA spikes[END_REF], gene networks [25,[START_REF] Newby | Isolating intrinsic noise sources in a stochastic genetic switch[END_REF], and neural networks [START_REF] Bressloff | Metastability in a stochastic neural network modeled as a velocity jump Markov process[END_REF] (Faggionato et al. were motivated by a model of molecular motors that is exactly solvable. Such a solvability condition also breaks down for molecular motors when local chemical signaling is taken into account [START_REF] Newby | Local synaptic signalling enhances the stochastic transport of motor-driven cargo in neurons[END_REF].)

Summary of main results

We will show (without the restrictions of Faggionato et al. [START_REF] Faggionato | Averaging and large deviation principles for fully-coupled piecewise deterministic Markov processes and applications to molecular motors[END_REF][START_REF] Faggionato | Non-equilibrium thermodynamics of piecewise deterministic Markov Processes[END_REF]) that the rate function J T ({x(t)} t∈[0,T ] ) of the non-Markov process {x(t)} t∈[0,T ] can be written in the form of an action

J T ({x(t)} t∈[0,T ] ) = T 0 L(x, ẋ)dt, (3.8)
with Lagrangian given by

L(x, ẋ) = µ(x, ẋ) ẋ -λ(x(t), µ(x, ẋ)), (3.9)
Here λ(x, µ) is the Perron eigenvalue of the linear equation

m∈Γ A nm (x)z m (x, µ) + µF n (x)z n (x, µ) = λ(x, µ)z n (x, µ), (3.10)
and µ = µ(x, ẋ) is the solution of the invertible equation

ẋ = m∈Γ ψ m (x, µ)F m (x), (3.11)
with

ψ m (x, µ) = z m (x, µ)R m (x, µ), (3.12)
where R m is the positive eigenvector of the adjoint equation (cf. equation (4.5))

m∈Γ A nm (x)R m (x, µ) + µF n (x)R n (x, µ) = λ(x, µ)R n (x, µ) (3.13) under the normalization m z m (x, µ)R m (x, µ) = 1.
Given the Lagrangian L, we can determine the Hamiltonian H according to the Fenchel-Legendre transformation Differentiating equation (3.10) with respect to µ gives

H(x, p) = sup y [(p -µ(x, y))y + λ(x, µ(x, y))] . ( 3
m∈Γ A nm (x) ∂z m (x, µ) ∂µ + [µF n (x) -λ(x, µ)) ∂z n (x, µ) ∂µ = ∂λ(x, µ) ∂µ -F n (x) z n (x, µ). (3.16)
Since the adjoint of the linear operator on the left-hand side has a one-dimensional null space spanned by R n , it follows from the normalization m z m R m = 1 that

∂λ(x, µ) ∂µ = m∈Γ ψ m (x, µ)F m (x) = y.
Equation (3.15) thus shows that p = µ, and we can identify p as the "conjugate momentum" of the Hamiltonian

H = λ(x, p), (3.17)
where λ(x, p) is the Perron eigenvalue of the linear equation (3.10) and its adjoint (3.13).

Evaluating the supremum

We now present the details of deriving the classical action (3.8) from the LDP given by equation (3.3). We first evaluate the supremum and then in section 3.3 we evaluate the infimum. We proceed by introducing the following Ansatz regarding the solution z = (z n , n ∈ Γ) of the variational problem (3.5) for fixed x and strictly positive measure ψ, namely, that it is an eigenvector of the following matrix equation:

A(x)z + q • z = λz (3.18)
for some bounded vector q = (q n , n ∈ Γ), Here, for any

a, b ∈ R K , [a • b] n ≡ [diag(a)b] n = a n b n .
Note that we are free to shift the vector q by a constant since, under the transformation q n → q n -c, the eigenvalue shifts by λ → λ -c and the eigenvector is unchanged. That is, for fixed x,

λ(x, q -c1 K ) = λ(x, q) -c. (3.19)
For simplicity, we choose c = q K and take Q = (q 1 , q 2 , . . . , q K-1 , 0) so that λ(x, Q) = λ(x, q) -q K and z = z(x, Q) are solutions of the matrix equation

A(x)z + Q • z = λz (3.20)
There are K -1 independent variables, q n , n = 1, . . . , K -1. One of the crucial features of the above Ansatz is that we can then ensure z n = z n (x, Q) is strictly positive by using the Perron-Frobenius theorem and taking λ = λ(x, Q) to be the Perron eigenvalue. Indeed, choosing κ such that

κ > max n=1,...,K-1 {|q n |}, (3.21)
it is clear that the new matrix A(x) + diag(Q) + κI K is irreducible and positive.

According to the Perron-Frobenius theorem, it has a unique strictly positive eigenvector z with the normalization m z m = 1, and this eigenvector is also an eigenvalue of A + diag(Q) (but with a shifted eigenvalue).

Applying equation (3.20) to the double sum in equation (3.5), which holds for strictly positive ψ, we find K m,n=1

ψ n W nm (x) 1 - z m z n = K n=1 ψ m K m=1 W nm (x) - 1 z n m (A nm (x) + δ nm ) z m = K m=1 ψ m [Q m -λ] = K-1 m=1 q m ψ m -λ(x, Q), since K m=1 ψ m = 1 and Q K = 0.
In order to ensure that we have found the true supremum with respect to z, we use equation (3.20) to show that the supremum satisfies equation (3.7). Its left-hand side is

K m=1 ψ n W nm (x) z m z n = ψ n z n K k=1 W nm (x)z m = ψ n z n (λz n -Q n z n ) + z n = ψ n (λ -Q n + 1),
whereas the right-hand side reads

z n K m=1 W mn (x) ψ m z m = z n K m=1 (A mn (x) + δ mn ) ψ m z m
Dividing the left and right-hand sides of the second equality by z n , we deduce that

ψ n = R n (x, Q)z n (x, Q) (3.22)
for each n = 1, . . . , K -1, where R = (R n , n ∈ Γ) is the corresponding unique strictly positive eigenvector (up to scalar multiplication) of the adjoint linear equation

A (x)R + Q • R = λR. (3.23)
such that z R = 1. Equation (3.22) ensures that ψ is a strictly positive measure and that

K m=1 ψ m = 1.
Assuming that such a solution exists, the corresponding function j(x, ψ) is given by

j(x, ψ) = K-1 n=1 q n ψ n -λ(x, Q), (3.24)
where we have dropped the constant term.

In the given variational problem there are K -1 independent variables ψ n , n = 1, . . . , K -1 with ψ K = 1-K-1 n=1 ψ n . Similarly, there are K -1 independent variables q n , n = 1, . . . , K -1. Therefore, equation (3.22) determines a mapping between the sets {q n , n = 1, . . . , K -1} and {ψ n , n = 1, . . . , K -1}. It remains to show that there exists a unique solution q for each ψ ∈ M + ([0, T ]) Γ , that is, the mapping is invertible. Differentiating equation (3.23) with respect to q m , m = 1, . . . , K -1, yields the inhomogeneous linear equation

L(x, Q) ∂R ∂q m ≡ A (x) + diag(Q) -λI K ∂R ∂q m = ∂λ ∂q m R -R m e m , (3.25)
where (e m ) n = δ m,n . Multiplying both sides of equation (3.25) on the left with z and using (3.20) we obtain

∂λ(x, Q) ∂q m = R m (x, Q)z m (x, Q), m = 1, . . . , K -1. (3.26)
Since R and z are strictly positive, λ(x, Q) is a monotonically increasing function of the q m . Moreover, equations (3.18) and (3.23) imply that, in the limit q l → ∞ with all other q m fixed, R l , z l → 1 and ∂λ/∂q l → 1. On the other hand, if q l → -∞ then R l , z l → 0 and the Perron eigenvalue becomes independent of q l with ∂λ/∂q l → 0. Hence, by continuity, for each ψ ∈ M + ([0, T ]) Γ there exists a vector Q such that

ψ n = R n (x, Q)z n (x, Q) for all n = 1, . . . , K -1.
For such a solution to be unique, the inverse function theorem implies that the Jacobian of the transformation must be invertible. Differentiating equation (3.26) with respect to q n shows that the Jacobian is equivalent to the Hessian of λ with respect Q, since

D mn (x, Q) ≡ ∂R m (x, Q)z m (x, Q) ∂q n = ∂ 2 λ(x, Q) ∂q m ∂q n (3.27)
for all m, n = 1, . . . , K -1. This also establishes that the Jacobian is a symmetric matrix with real eigenvalues. Invertibility follows from the convexity of the function J (c) defined by equation (3.6). That is, differentiating

j(x, ψ) = J (c), c nm = ψ n W nm (x),
with respect to ψ for fixed x gives

W nm (x)W mm (x) ∂ 2 J (c) ∂c nm ∂c mm = ∂ 2 j(x, ψ) ∂ψ n ∂ψ m (3.28)
with n, m = m . On the other hand, differentiating equation (3.24) for j(x, ψ) with respect to ψ n gives ∂j(x, ψ)

∂ψ n = q n + K-1 j=1 ψ j - ∂λ ∂q j ∂q j ∂ψ n = q n ,
and so

∂ 2 j(x, ψ) ∂ψ n ∂ψ m = ∂q n ∂ψ m = [D -1 ] nm .
Hence, convexity of J (c) together with irreducibility of the non-negative transition matrix W means that the Jacobian is invertible and positive definite.

In summary, we have shown that for a strictly positive measure ψ and fixed x, a unique solution q n = q n (x, ψ) exists for all n = 1, . . . , K -1 and we have solved the first variational problem by identifying z with the unique (up to scalar multiplication), strictly positive eigenfunction of the matrix equation (3.20). Now suppose ψ is a nonnegative rather than a strictly positive measure measure, that is, ψ m = 0 for at least one state m ∈ Γ. In this case c nm = ψ n W nm (x), n = m, is not irreducible and lemma 1 no longer applies. However, as proven by Faggionato et al. [START_REF] Faggionato | Averaging and large deviation principles for fully-coupled piecewise deterministic Markov processes and applications to molecular motors[END_REF], the function J (c) is continuous with respect to c. Hence, assuming that the form of the rate function (3.3) still holds (which isn't necessarily true), we can take a sequence of strictly positive measures ψ (l) on Γ such that ψ (l) n → ψ n for each n ∈ Γ. This implies that (for fixed

x) c[ψ (l) ] → c[ψ] and j(x, ψ (l) ) = J (c[ψ (l) ]) → J (c[ψ]) = j(x, ψ).
so that one can extend equation (3.24) to non-negative measures by taking q n (x, ψ) = lim l→∞ q n (x, ψ (l) ).

Example. We will illustrate the Perron eigenvalue solution to the supremum variational problem by considering an example for K = 2. Let us take the transition matrix to be

W = 1/2 1/3 1/2 2/3 Consider the eigenvalue equation Wz + q • z = λz,
where we have absorbed the diagonal terms k=1,2 W km into the definition of q m . The resulting characteristic equation is a quadratic in λ and the leading or Perron eigenvalue is given by

λ = q 1 + q 2 2 + 7 12 + 1 2 (q 1 -q 2 ) 2 -(q 1 -q 2 )/3 + 25/36.
It follows that

ψ 1 ≡ ∂λ ∂q 1 = 1 2 + f (q 1 -q 2 ) ψ 2 ≡ ∂λ ∂q 2 = 1 2 -f (q 1 -q 2 )
with

f (q) = 1 4
2q -1/3 q 2 -q/3 + 25/36 .

Note that ψ 1 + ψ 2 = 1 as required. The function f (q) is a monotonically increasing function of q with f (-∞) = -1/2 and f (∞) = 1/2. Thus, one can find a unique, finite value of q = q 1 -q 2 for all ψ 1 ∈ (0, 1), that is, for all strictly positive ψ. In the case of a non-negative ψ with ψ 1 = 0 or ψ 2 = 0, we have q → ±∞.

Evaluating the infimum

The next step is to substitute for j(x, ψ) in the rate function (3.4), which gives

J T ({x(t)} t∈[0,T ] ) = inf ψ: ẋ= K n=1 ψnFn(x) T 0 K-1 n=1 q n (t)ψ n (t) -λ(x(t), Q(t)) dt, (3.29)
with q n (t) = q n (x(t), ψ(t)) and K m=1 ψ m = 1. In order to solve this variational problem, we introduce a Lagrange multiplier µ(t) and set

J T ({x(t)} t∈[0,T ] ) = inf ψ,µ S[x, µ, ψ], (3.30)
where

S[x, µ, Ψ] = T 0 K-1 n=1 q n (t)ψ n (t) -λ(x(t), Q(t)) (3.31) +µ(t)( ẋ - K-1 n=1 [F n (x) -F K (x)]ψ n (t) -F K (x)) dt,
q m = µ[F m (x) -F K (x)], (3.33a)
ψ m (x, µ) = z m (x, µ)R m (x, µ), (3.33b)
for all m = 1, . . . K -1, with R n (x, µ), z n (x, µ) the positive eigenvectors of the matrix equations

A R + µF(x) • R = λR (3.34)
and its adjoint, respectively. Here

λ = λ(x, µ) ≡ λ(x, Q|q m = µ(F m -F K ), m = 1, . . . , K -1).
The corresponding function j(x, ψ) becomes

j(x, ψ) = K n=1 ψ n (x, µ)F n (x) -λ(x, µ), (3.35)
We have thus established Eqs. (3.10) and (3.13). The final step is to show that for each x, the equation

ẋ = K n=1 F n (x)ψ n (x, µ)
is invertible so that the function µ = µ(x, ẋ) exists, and the rate function has the required Lagrangian form (3.8). From the inverse function theorem we require that

L(x, µ) ∂R ∂µ = ∂λ ∂µ R -F • R. (3.36)
with L(x, µ) = L(x, p)| pm=µFm . Using the same arguments as previously we obtain

∂λ(x, µ) ∂µ = K n=1 F n (x)z n (x, µ)R n (x, µ) (3.37)
Hence, we require In the above construction we considered the case of one-dimensional piecewise deterministic dynamics by taking x ∈ R. However, it is relatively straightforward to derive the corresponding LDP for x ∈ R d . When the internal state is n, the system now evolves according to the ODE ẋ = F n (x), (3.38) where the vector field F n : Ω → R d is a continuous function, locally Lipschitz. That is, given a compact subset K of Ω, there exists a positive constant with Lagrangian given by L(x, ẋ) = µ(x, ẋ), ẋ -λ(x(t), µ(x, ẋ)), Here λ(x, µ) is the Perron eigenvalue of the linear equation

∂ 2 λ(x, µ) ∂µ 2 ≡ K-1 m,n=1 ∂ 2 λ(x, Q) ∂q m ∂q n qm=µ(Fm(x)-F K (x)) × (F m (x) -F K (x))(F n (x) -F K (x)) =
K n (K) such that |F n (x) -F n (y)| ≤ K n (K)|x -y|, ∀x, y ∈ Ω. ( 3
A(x)z(x, µ) + (F(x)µ) • z(x, µ) = λ(x, µ)z(x, µ). (3.41) Here, for any a, b ∈ R K , [a • b] n ≡ [diag(a)b] n = a n b n , n = 1, • • • , K, (3.42

Extension to a multi-scale process on R

So far we have assumed that the slow process is piecewise deterministic. However, one of the useful features of taking the Lagrangian LDP [START_REF] Faggionato | Averaging and large deviation principles for fully-coupled piecewise deterministic Markov processes and applications to molecular motors[END_REF][START_REF] Faggionato | Non-equilibrium thermodynamics of piecewise deterministic Markov Processes[END_REF] as our starting point is that it is relatively straightforward to extend our analysis to the case where the slow process is a piecewise SDE. First, recall that the key idea behind the Faggionato et al. Lagrangian construction is that the slow dynamical process coupled to the fast Markov chain on Γ rapidly samples the different discrete states of Γ according to some non-negative measure ψ. In order to extend this construction to a piecewise SDE, it is necessary to take account of the fact that there are now two levels of stochasticity.

That is, after averaging the transition rates of the drift and variance of the SDE with respect to a given measure ψ, the resulting system is still stochastic. Since the slow system operates in a weak noise regime, it follows that one can apply an LDP to the slow system for a given ψ. The LDP for the full system is then obtained by combining the rate function of the slow system with the supremum rate function for ψ. Here we sketch how to extend the analysis to a piecewise SDE. (Approaching large deviation theory for multi-scale stochastic processes in terms of solutions to an eigenvalue problem has also been considered by Feng and Kurtz [START_REF] Feng | Large Deviations for Stochastic Processes[END_REF].) Consider the piecewise Ito SDE

dX(t) = F n (X) + √ σ n (X)dW (t), (3.46)
where n ∈ Γ and W (t) is a Wiener process. The drift term F n (X) and diffusion term σ n (X) are both taken to be Lipschitz. When the SDE is coupled to the fast discrete process on Γ, the stochastic dynamics is described by a differential Chapman-Kolmogorov equation. That is, writing

P{X(t) ∈ (x, x + dx), n(t) = n|x 0 , n 0 ) = ρ(x, n, t)dx, we have (see also equation (2.6)) ∂ρ(x, n, t) ∂t = - ∂ ∂x [F n (x)ρ(x, n, t)] + 2 ∂ 2 ∂x 2 [σ 2 n (x)ρ(x, n, t)] + 1 m A nm (x)ρ(x, m, t). (3.47)
We define the measure space M + ([0, T ]) as before, but now modify the definition of the subspace Y x0 ⊂ C([0, T ]) × M + ([0, T ]) Γ by taking it to be the set of stochastic trajectories satisfying

dX(t) = K n=1 ψ n F n (X) + √ K n=1 ψ n (t)σ 2 n (X) dW (t) ≡ F (X, ψ) + √ σ 2 (X, ψ)dW (t) (3.48)
for ψ ∈ M + ([0, T ]) Γ . Such a space contains the set of trajectories of the SDE (3.46) with ψ n (t) given by equation (3.1) and n(t) evolving according to the Markov chain on Γ. Consider a particular realization of the Wiener process W (t) on [0, T ], which is independent of {x(t), ψ(t)} [0,T ] . For a given realization, one can write down an LDP along identical lines to the case of a piecewise deterministic system, see equation (3.3)). Assuming that we can then average with respect to the Wiener process, we obtain an additional contribution to the rate function so that

J T ({x(t)} t∈[0,T ] ) = inf {ψ(t)} t∈[0,T ] J T ({(x(t), ψ(t))} t∈[0,T ] ). (3.49) with J T ({(x(t), ψ(t))} t∈[0,T ] ) = T 0 j(x(t), ψ(t)) dt + A T ({(x(t), ψ(t))}) t∈[0,T ] ) (3.50)
for j(x, ψ) given by equation (3.5) and

A T ({x(t), ψ(t)} t∈[0,T ] ) = T 0 ( ẋ -F (x, ψ)) 2 2σ 2 (x, ψ) dt. (3.51)
Equation (3.51) is the well known action functional for a one-dimensional SDE in the case of a fixed, time-independent ψ [START_REF] Freidlin | Random perturbations of dynamical systems[END_REF], see section 2. One can evaluate the rate function (3.49) along similar lines to the case of a piecewise deterministic system (see appendix A.3) to obtain a classical action with Hamiltonian given by the Perron eigenvalue λ(x, p) of the linear equation

m A nm (x)R m (x, p) + (pF n (x) + p 2 σ 2 n /2)R n (x, p) = λ(x, p)R n (x, p). (3.52)
The above analysis can be extended to the LDP (3.49) for a multi-scale stochastic process. In particular, using equation (3.24) we have

J T ({x(t)} t∈[0,T ] ) = inf {ψ(t)} t∈[0,T ] S[x, ψ] (3.53)
where

S[x, ψ] = T 0 K-1 n=1 q n (t)ψ n (t) -λ(x(t), Q(t)) + ( ẋ - K n=1 ψ n F n (x)) 2 2 K-1 n=1 ψ n σ 2 n dt, (3.54)
with q n (t) = q n (x(t), ψ(t)), K m=1 ψ m = 1 and λ the Perron eigenvalue of equation (3.18). Taking the infimum with respect to ψ k , k = 1, . . . , K -1, gives

0 = δS δψ k = q k - ( ẋ - K n=1 ψ n F n (x))[F k (x) -F K (x)] K n=1 ψ n σ 2 n (3.55) - ( ẋ - K n=1 ψ n F n (x)) 2 [σ 2 k (x) -σ 2 K (x)] 2[ K n=1 ψ n σ 2 n ] 2 + K-1 n=1 ψ n - ∂λ ∂q n ∂q n ∂ψ k
Introducing the new variables

µ = K n=1 ψ n F n (x), σ 2 = K n=1 ψ n σ 2 n , p = ẋ -µ σ 2 , (3.56)
and noting that ∂λ/∂q m = ψ m , we have

q k = p[F k (x) -F K (x)] + p 2 2 [σ 2 k (x) -σ 2 K (x)] (3.57)
and

S[x, ψ] = T 0 p K n=1 F n ψ n + p 2 2 K n=1 ψ n σ 2 n -λ(x, p) + ( ẋ -µ) 2 2σ 2 dt = T 0 pµ + p 2 2 σ 2 -λ(x, p) + p 2 ( ẋ -µ) dt = T 0 [p ẋ -λ(x, p)] (3.58)
with λ the Perron eigenvalue for the linear equation (which is independent of µ)

n A mn R n + (pF m (x) + p 2 σ 2 m /2)R m = λR m . (3.59)
Finally µ and σ 2 are determined from the identities

µ = n R n (x, p)z n (x, p)F n (x), σ 2 = n R n (x, p)z n (x, p)σ 2 n (x). (3.60)

Classical Hamiltonian and the WKB approximation of the stationary state

As indicated in the introduction, for non-hybrid stochastic processes the quasipotential of WKB theory satisfies the same Hamilton-Jacobi equation as obtained from the variational principle of large deviation theory. It turns out that such a connection is more subtle in the case of an SHS. In order to show this, we consider a WKB approximation of the stationary state ρ ss of the CK equation (2.5) (assuming it exists), with

Lρ ss (x, n) ≡ ∂F n (x)ρ ss (x, n) ∂x - 1 ε m∈Γ A nm (x)ρ ss (x, m) = 0. ( 4.1)
The WKB approximation of ρ ss takes the form

ρ ss (x, n) ∼ η(x, n) exp - Φ(x) , (4.2)
Substituting into equation (4.1) yields m∈Γ

A nm (x) + Φ (x)δ n,m F m (x) η(x, m) = dF n (x)η(x, n) dx , (4.3)
where Φ = dΦ/dx. Introducing the asymptotic expansions η ∼ η (0) + η (1) + . . . and Φ ∼ Φ 0 + Φ 1 + . . ., the leading order equation is

m∈Γ A nm (x)η (0) (x, m) = -Φ 0 (x)F n (x)η (0) (x, n). (4.4)
As it stands, it is not clear that (4.4) has a solution for which η (0) (x, m) is positive for all x, m, and the relationship to the Hamiltonian structure of large deviation theory is not explicit. However, the structure of equation (4.4) is a reduced version of equation (3.13). This suggests introducing the family of eigenvalue equations

m∈Γ A nm (x)R m (x, p) + pF n (x)R n (x, p) = λ(x, p)R n (x, p), (4.5)
which are parameterized by the pair (x, p) with p an auxiliary variable and λ(x, p) the Perron eigenvalue. Comparison of equation (4.4) with (4.5) then shows that we can make the identifications Φ 0 (x) = p, η (0) (x, m) = R m (x, p) and λ(x, p) = 0. It immediately follows that η (0) is positive. Hence, the quasi-potential is the solution of the Hamilton-Jacobi equation

λ(x, Φ 0 (x)) = 0, (4.6)
This is equivalent to finding zero energy solutions of Hamilton's equations

ẋ = ∂λ(x, p) ∂p , ṗ = - ∂λ(x, p) ∂x , (4.7)
and identifying Φ 0 as the action along the resulting solution curve (x(t), p(t)):

Φ 0 (x) = x x p(x )dx . (4.8)
Note that here t is the parameter of a curve rather than physical time. One of the major applications of WKB methods is to solving escape problems for stochastic process in the weak-noise limit, both for non-hybrid systems [START_REF] Matkowsky | The exit problem for randomly perturbed dynamical systems[END_REF][START_REF] Naeh | A direct approach to the exit problem[END_REF][START_REF] Maier | Limiting exit location distribution in the stochastic exit problem[END_REF][START_REF] Hanggi | Bistable systems: master equation versus Fokker-Planck modeling[END_REF][START_REF] Knessl | An asymptotic theory of large deviations for Markov jump processes[END_REF][START_REF] Dykman | Large fluctuations and optimal paths in chemical kinetics[END_REF], and hybrid systems [START_REF] Keener | Perturbation analysis of spontaneous action potential initiation by stochastic ion channels[END_REF][START_REF] Newby | An asymptotic analysis of the spatially inhomogeneous velocity-jump process[END_REF][START_REF] Newby | Isolating intrinsic noise sources in a stochastic genetic switch[END_REF][START_REF] Newby | Breakdown of fast-slow analysis in an excitable system with channel noise[END_REF][START_REF] Bressloff | Stochastic hybrid model of spontaneous dendritic NMDA spikes[END_REF]. For example, suppose that the mean-field equation (2.8) is bistable with two stable fixed-points x ± separated by an unstable fixed point x 0 , see also Fig. 1. Given the quasi-potential Φ 0 , the mean first passage time τ to escape from x -, say, can be calculated by considering higher order terms in the WKB approximation, and using matched asymptotics to deal with the absorbing boundary at x 0 . One finds that τ takes the general Arrhenius form [START_REF] Keener | Perturbation analysis of spontaneous action potential initiation by stochastic ion channels[END_REF][START_REF] Newby | An asymptotic analysis of the spatially inhomogeneous velocity-jump process[END_REF] τ ∼ χ(x 0 , x -)

|Φ 0 (x 0 )|Φ 0 (x -) e [Φ0(x0)-Φ0(x-)]/ ,
where χ is an appropriate prefactor. Hence, the WKB method provides a powerful calculational tool. On the other hand, there is no a priori justification for interpreting the quasi-potential and its associated Hamiltonian in terms of an underlying variational principle for optimal paths in the space of stochastic trajectories. This becomes crucial when solving escape problems in higher dimensions, since a metastable state is now surrounded by a non-trivial boundary (rather than a single point) and one needs to determine the relative weighting of optimal paths crossing different points on the boundary. Establishing the connection between WKB analysis and large deviation theory provides such a variational principle.

Applications to stochastic ion channels

In this final section we illustrate the Hamiltonian structure of stochastic hybrid systems by considering a few explicit models taken from neuroscience.

Binary model

Before considering conductance-based models of membrane voltage fluctuations, we begin with the simple example of a binary stochastic hybrid process (two discrete states n = 0, 1), which was analyzed in some detail by Faggionato et al [START_REF] Faggionato | Non-equilibrium thermodynamics of piecewise deterministic Markov Processes[END_REF] using a different method. The latter authors exploited the fact that the model is exactly solvable, in the sense that the stationary density of the corresponding Chapman-Kolmogorov equation can be computed explicitly, and used a fluctuation-dissipation theorem to determine the Hamiltonian and quasipotential. Here we will obtain the same results more directly by calculating the Perron eigenvalue. One biological application of the binary model is to the bidirectional transport of a molecular motor along a onedimensional microtubular track, in which x represents the spatial location of the motor on the track and the two discrete states represent the motor moving either towards the + end or -end of the track. (In more complex models, the discrete space Γ represents multiple internal conformational states of the motor, each of which has an associated velocity on the track [START_REF] Bressloff | Metastability in a stochastic neural network modeled as a velocity jump Markov process[END_REF].) Suppose that the continuous variable evolves according to piecewise dynamics on some finite interval (a, b),

ẋ = F n (x), n = 0, 1, (5.1)
with F 0 , F 1 continuous and locally Lipschitz. Suppose that F 0 , F 1 are non-vanishing within the interval (a, b), and F n (a) ≥ 0, F n (b) ≤ 0 for n = 0, 1; the dynamics is then confined to (a, b). Denote the transition rates of the two-state Markov chain by ω ± (x) with

{n = 0} ω+(x) ω-(x) {n = 1}.
The stationary measure of the Markov chain is given by

ρ(x, 0) = ω -(x) ω -(x) + ω -(x) , ρ(x, 1) = ω + (x) ω -(x) + ω -(x)
.

(5.

2)

The linear equation (3.13) can be written as the two-dimensional system

-ω + (x) + pF 0 (x) ω -(x) ω + (x) -ω -(x) + pF 1 (x) R 0 R 1 = λ R 0 R 1 .
(5.

3)

The corresponding characteristic equation is

0 = λ 2 + λ[ω + (x) + ω -(x) -p(F 0 (x) + F 1 (x))] + (pF 1 (x) -ω -(x))(pF 0 (x) -ω + (x)) -ω -(x)ω + (x).
It follows that the Perron eigenvalue is given by

λ(x, p) = 1 2 Σ(x, p) + Σ(x, p) 2 -4γ(x, p) (5.4)
where

Σ(x, p) = p(F 0 (x) + F 1 (x)) -[ω + (x) + ω -(x)],
and

γ(x, p) = (pF 1 (x) -ω -(x))(pF 0 (x) -ω + (x)) -ω -(x)ω + (x).
A little algebra shows that D(x, p) ≡ Σ(x, p) 2 -4γ(x, p) = [p(F 0 -F 1 ) -(ω + -ω -)] 2 + ω + ω -> 0 so that as expected λ is real. Hence, from Hamilton's equations ẋ = ∂λ(x, p) ∂p

= F 0 (x) + F 1 (x) 2 + ∂D(x, p) ∂p 1 2 D(x, p) (5.5) = F 0 (x) + F 1 (x) 2 + F 0 (x) -F 1 (x) 2 p(F 0 -F 1 ) -(ω + -ω -) [p(F 0 -F 1 ) -(ω + -ω -)] 2 + ω + ω -
which is the same result as obtained in example 10.4 of Faggionato et al [START_REF] Faggionato | Non-equilibrium thermodynamics of piecewise deterministic Markov Processes[END_REF]. Moreover, writing

ẋ = F 0 (x)ψ 0 (x) + F 1 (x)ψ 1 (x),
we see that

ψ 0 (x) = 1 2 1 + p(F 0 -F 1 ) -(ω + -ω -) [p(F 0 -F 1 ) -(ω + -ω -)] 2 + ω + ω - , (5.6)
and

ψ 1 (x) = 1 2 1 - p(F 0 -F 1 ) -(ω + -ω -) [p(F 0 -F 1 ) -(ω + -ω -)] 2 + ω + ω - , (5.7)
so that ψ 0,1 (x) ≥ 0 with ψ 0 (x) + ψ 1 (x) = 1.

Stochastic Na + ion channels and the initiation of spontaneous action potentials

An important example of stochastic hybrid systems at the single-cell level concerns a conductance-based model of a neuron, in which the stochastic opening of membrane ion channels generates a stochastic ionic current that drives the membrane voltage. It is then possible that ion channel noise induces spontaneous action potentials (SAPs), which can have a large effect on a neuron's function [START_REF] Fox | Emergent collective behavior in large numbers of globally coupled independent stochastic ion channels[END_REF][START_REF] Chow | Spontaneous action potentials due to channel fluctuations[END_REF][START_REF] Keener | Perturbation analysis of spontaneous action potential initiation by stochastic ion channels[END_REF][START_REF] Goldwyn | The what and where of adding channel noise to the Hodgkin-Huxley equations[END_REF][START_REF] Buckwar | An exact stochastic hybrid model of excitable membranes including spatio-temporal evolution[END_REF][START_REF] Newby | Breakdown of fast-slow analysis in an excitable system with channel noise[END_REF][START_REF] Bressloff | Stochastic hybrid model of spontaneous dendritic NMDA spikes[END_REF]. If SAPs are too frequent, a neuron cannot reliably perform its computational role. Hence, ion channel noise imposes a fundamental limit on the density of neural tissue. Smaller neurons must function with fewer ion channels, making ion channel fluctuations more significant and more likely to cause a SAP. Here we will consider the simple case of a single type of ion channel, namely, a fast sodium (Na) channel, which was previously analyzed using WKB methods [START_REF] Keener | Perturbation analysis of spontaneous action potential initiation by stochastic ion channels[END_REF]. Let x(t) denote the membrane voltage of the neuron at time t and N be the fixed number of sodium channels. We assume that each channel can either be open (O) or closed (C), and can switch between each state according to the kinetic scheme

C α(x) β(x)
O.

(5.8)

with voltage-dependent transition rates. (A more detailed biophysical model would need to treat each ion channel as a cluster of subunits rather than a single unit. In other words, the Markov chain of events associated with opening and closing of an ion channel would involve transitions between more than two internal states,.) The stochastic membrane voltage is taken evolves according to the piecewise deterministic equation

dx dt = F n (x) ≡ n N f (x) -g(x), (5.9)
where n is the number of open ion channels at time t, and

f (x) = g Na (V Na -x), g(x) = -g L [V L -x] -I.
Here g Na is the maximal conductance of a sodium channel and V Na is the corresponding membrane reversal potential. Similarly, g L and V l are the effective maximal conductance and reversal potential of any other currents, which are assumed to be independent of the opening and closing of ion channels, and I is an external curent. The four quantities (g Na , g L , V Na , V L ) are taken to be constants. Since the right-hand side of (5.9) is negative for large x and positive for small x, it follows that the voltage x is confined to some interval Ω = [x L , x R ]. The function F n (x) is clearly continuous and locally Lipschitz.

In this example the space Γ of discrete states is the set of integers {n = 0, 1, . . . , N } and the Markov chain is given by a birth-death process:

n → ω+(n,x)/ n + 1, n → ω-(n,x)/ n -1 (5.10)
with transition rates

ω + (x, n) = α(x)(N -n), ω -(x, n) = β(x)n.
(5.11)

The small parameter reflects the fact that sodium channels open at a much faster rate than the relaxation dynamics of the voltage [START_REF] Keener | Perturbation analysis of spontaneous action potential initiation by stochastic ion channels[END_REF]. It follows that the matrix A(x) for fixed x is tridiagonal matrix with

A n-1,n (x) = ω + (x, n -1), A n+1,n (x) = ω -(x, n + 1) (5.12a) A nn (x) = -ω + (x, n) -ω -(n), (5.12b)
for n = 0, 1, . . . , N . It is straightforward to show that the Markov chain is ergodic with unique invariant measure (for fixed n) given by

ρ(x, n) = N ! (N -n)!n! a(x) n b(x) N -n , (5.13)
with

a(x) = α(x) α(x) + β(x) , b(x) = β(x) α(x) + β(x)
.

(5.14)

The above stochastic hybrid system satisfies all of the conditions specified in section 3. Hence, the law of large numbers implies that in the mean-field limit → 0, we obtain the deterministic kinetic equation

dx dt = F (x) ≡ a(x)f (x) -g(x) (5.15)
where

a(x) = n /N, n = N n=1 nρ(x, n),
and ρ is the stationary density (5.13). One of the features of the averaged model is that it can exhibit bistability for a range of physiologically reasonable parameter values. This is illustrated in Fig. 1, where we plot the deterministic potential U (x) = -dF /dx as a function of x. Here x -represents a resting state of the neuron, whereas x + represents an active state; noise-induced transitions from x -to x + can be interpreted in terms of the initiation of a spontaneous action potential. Elsewhere, WKB methods and matched asyptotics have been used to calculate the MFPT to escape from x -.

Here, we will focus on the quasi-potential and its relation to the Perron eigenvalue U(x)

x

Stochastic Morris-Lecar model

One of the major simplifications of the above model is to assume that the slow potassium channels are frozen. If one now incorporates the slow opening and closing of these channels, then the underlying deterministic system becomes excitable rather than bistable. That is, there is a single stable fixed point such that for small stimuli the voltage returns directly to rest, whereas for stronger stimuli the voltage makes a large detour before returning to rest, which corresponds to an action potential. There is no longer a well-defined, unique firing threshold. A simple deterministic model of neural excitability is the Morris-Lecar (ML) model [START_REF] Morris | Voltage oscillations in the barnacle giant muscle fiber[END_REF]: ẋ = a(x)f Na (x) + yf K (x) -g(x)

(5.23) ẏ = y ∞ (x) -y τ y (x) , (5.24)

where x is voltage and y represents the fraction of open K + channels. The dynamics of this system can be explored using phase-plane analysis as illustrated in Fig. 3. A slow/fast analysis of the deterministic system suggests that the initiation of an action potential occurs without any change in w, thus motivating the analysis of Keener and Newby [START_REF] Keener | Perturbation analysis of spontaneous action potential initiation by stochastic ion channels[END_REF]. However, it turns out that this adiabatic approximation breaks down when stochastic fluctuations in the opening and closing of K + channels are taken into account. This can be established by extending the WKB analysis outlined in section 4 to a stochastic version of the ML model [START_REF] Newby | Breakdown of fast-slow analysis in an excitable system with channel noise[END_REF]. Since one now has two continuous variables x and y in the deterministic limit, it follows that stochastic trajectories in the phase plane correspond to characteristic projections of an underlying Hamiltonian dynamical system. In general, it is difficult to solve FPT problems in more than one Green/blue curves represent an action potential trajectory in the limit of slow y.

dimension. In the case of a metastable state with a well-defined basin of attraction, one has to calculate the MFPT to cross the separatrices forming the boundary of the basin of attraction. There is an additional level of complexity for an excitable system, due to the fact that there is no well-defined deterministic separatrix. Interestingly, one finds that the stochastic ML model has an effective separatrix that any stochastic trajectory has to cross in order to generate a stochastic action potential [START_REF] Newby | Breakdown of fast-slow analysis in an excitable system with channel noise[END_REF], see also [START_REF] Khovanov | Noise-induced escape in an excitable system[END_REF].

Here we focus on calculating the Hamiltonian associated with the LDP for the stochastic Morris-Lecar model. Suppose (5.28b)

In contrast to the bistable sodium ion channel model of section 5.2, we cannot treat both the Na and K channel kinetics as fast, and therefore we cannot develop a variational problem by scaling all transition rates in terms of a small parameter and applying the analysis of section 5. In fact, rather than a piecewise deterministic system, we now have a multi-scale stochastic system, in which both fast and slow

. 14 )

 14 Evaluating the right-hand side yields the equation p -µ(x, y) ψ m (x, µ)F m (x).

and we have imposed the constraint K m=1 ψ m = 1 .

 1 The variational problem can now be expressed in terms of functional derivatives of Ss)[F m (x(s)) -F K (x(s))] = 0 (3.32b) for m = 1, . . . K -1. Combining with equations (3.22) and (3.26), we obtain the following solution to the variational problem in terms of µ:

  K m=1 ∂ψ m (x, µ) ∂µ F m (x) = 0 for all x ∈ Ω. Following along identical lines to the analysis of equation (3.25), we differentiate the linear equation (3.34) with respect to µ to give

  0. This holds since the Jacobian D of equation (3.27) is invertible and F m (x) = F K (x) for at least one m = K. Finally, from equation (3.21), we require µ to be bounded, that is, there exists a κ for which κ > µ max m=1,...,K {|F m (x)|} for all x ∈ Ω. Hence our derivation of the classical action (3.8) for the LDP rate function (3.4) is complete.

3. 4 .

 4 Extension to higher-dimensions, x ∈ R d , d > 1.

. 39 )

 39 The rate function of the LDP (3.3) becomesJ T ({x(t)} t∈[0,T ] ) =

  ) the K-dimensional vector F(x)µ is the product of the K ×d matrix F(x) whose K rows are the d-dimensional vectors F m (x), m = 1, • • • , K, with the d-dimensional vector µ, and µ = µ(x, ẋ) is the solution of the invertible equationẋ = m∈Γ ψ m (x, µ)F m (x), (3.43) with ψ m (x, µ) = z m (x, µ)R m (x, µ),(3.44) where z is the positive eigenvector of the adjoint equationA (x)R(x, µ) + (F(x)µ) • R(x, µ) = λ(x, µ)R(x, µ)(3.45) under the normalizations n z n = 1 and m z m (x, µ)R m (x, µ) = 1. The associated Hamiltonian system is the 2d-dimensional phase space (x, p) with Hamiltonian H(x, p) = λ(x, p) and Hamilton's equations ẋr = ∂λ ∂p r , ṗr = -∂λ ∂x r , r = 1, . . . , d.

Figure 1 :Figure 2 :

 12 Figure 1: Plot of deterministic potential U (x) as a function of voltage x for different values of the external stimulus current I. Parameter values are VNa = 120 mV, VL = -62.3 mV, gNa = 4.4 mS/cm 2 , gL = 2.2 mS/cm 2 , and α(x) = β exp[(x -v1)/v2] with β = 0.8 s -1 , v1 = -1.2mV , v2 = 18mv.

Figure 3 :

 3 Figure 3: Deterministic phase plane dynamics. Thick curves show the nullclines: ẋ = 0 as grey and ẏ = 0 as black. Black stream lines represent deterministic trajectories.Green/blue curves represent an action potential trajectory in the limit of slow y.

  that at time t there are n = 0, 1, • • • , N open Na channels and m = 0, 1, • • • , M open K channels. The membrane voltage x then evolves asdx dt = F (x, m, n) ≡ n N f Na (x) + m M f K (x) + f L (x) + I. (5.25)We assume that each channel is either open or closed and switches between each state according toO αi(x) βi(x)C, i = Na, K.(5.26)The space Γ of discrete states is the set of integer pairs (n, m) with 0 ≤ n ≤ wN, 0 ≤ m ≤ M and the Markov chain is given by a birth-death process(n, m) → ω ± N a (n,x) (n ± 1, m), (n, m) → Na (n, x) = nβ Na (x), ω + Na (n, x) = (N -n)α Na (x),(5.28a)ω - K (m, x) = mβ K (x), ω + K (m, x) = (M -m)α K (x).

z∈(0,∞) Γ (n,m)∈Γ×Γ

Acknowledgements

We thank Jay Newby for drawing our attention to the work of Yuri Kifer [26] and helpful discussion. The work was conducted while PCB was visiting the NeuroMathComp group where he holds an INRIA International Chair. OF was partially supported by the European Union Seventh Framework Programme (FP7/2007-2013) under grant agreement no. 269921 (BrainScaleS), no. 318723 (Mathematics), and by the ERC advanced grant NerVi no. 227747.

Substituting the explicit expressions for A and F n (x) into equation (3.13), yields the following equation for the Perron eigenvalue λ and the right eigenvector R:

Consider the trial solution [START_REF] Bressloff | Stochastic hybrid model of spontaneous dendritic NMDA spikes[END_REF] R n (x, p) = Γ(x, p) n (N -n)!n! , (5.17) which yields the following equation relating Γ and λ: nα

Collecting terms independent of n and terms linear in n yields the pair of equations

and

Eliminating Γ from these equation gives

This yields a quadratic equation for λ of the form

Given the Perron eigenvalue, we can determine the quasipotential Φ(x) by solving the Hamiton-Jacobi equation λ(x, ∂ x Φ) = 0. Equation (5.20) then yields the reduced Hamilton-Jacobi equation

The latter is precisely the equation for the quasipotential previously derived using WKB methods [START_REF] Newby | Breakdown of fast-slow analysis in an excitable system with channel noise[END_REF]. It has the following pair of solutions for Φ 0 = ∂ x Φ 0 :

The trivial solution Φ 0 = constant occurs along deterministic trajectories, which converge to the fixed point, whereas the non-trivial solution for Φ 0 (x) occurs along the most likely escape trajectories. In Fig. 2 we show solutions to Hamilton's equations in the (x, p)-plane, highlighting the zero energy maximum likelihood curve linking x - and x 0 .

processes are intrinsically stochastic. Multi-scale stochastic processes also arise in models of gene regulatory networks [START_REF] Newby | Isolating intrinsic noise sources in a stochastic genetic switch[END_REF]. We will proceed along similar lines to Ref. [START_REF] Brooks | Quasicycles in the stochastic hybrid Morris-Lecar neural model[END_REF] by treating n(t) as a fast variable with α Na , β Na = O(1/ ), and treating y(t) = m(t)/M as a continuous (recovery) variable with M = 1/ . We can then derive a piecewise hybrid SDE by carrying out a system size expansion with respect to y.

Setting M Ω ± (x, y) = ω ± m (M y, x), we obtain the SDE [START_REF] Brooks | Quasicycles in the stochastic hybrid Morris-Lecar neural model[END_REF]

with

We can now determine the Hamiltonian of the associated LDP by combining our analysis in sections 3.4 and 3.5. That is, H(x, p) = λ(x, p) with x = (x, y), p = (p x , p y ), and λ is the Perron eigenvalue of the linear equation

where the matrix A is the same as the Na model of section 5.2. Equation (5.31) can be solved along similar lines to (5.16) using the Ansatz R n (x, p) = Γ(x, p) n /(n!(N -n)!).

Collecting terms linear in n gives Γ(x, p) = α Na (x) -1 N (p x g(x, w) + h(x, y, p y ) -λ(x, p)),

where g(x, y) = yf K (x) + f L (x) and h(x, y, p y ) = p y [Ω + (x, y) -Ω -(x, y)] + 1 2 p 2 y σ 2 (x, y).

(5.32)

On the other hand, collecting terms independent of n and substituting for Γ(x, p) gives the following quadratic equation for λ:

and

with y ∞ (x) = α K (x)/(α K (x) + β K (x)). Note, in particular, for escape problems we are interested in zero energy solutions of Hamilton's equations, which reduce to solutions of h 1 (x, p) = 0. Our derivation of the Hamiltonian based on an LDP is equivalent to one obtained previously using formal WKB methods [START_REF] Newby | Breakdown of fast-slow analysis in an excitable system with channel noise[END_REF].


Cyril Cohen
email: cyril.cohen@inria.fr

Damien Rouhling
email: damien.rouhling@inria.fr

A renement-based approach to large scale reection for algebra

Large scale reection tactics are often implemented with ad-hoc data-structures and in a way which is specic to the problematic. This makes it hard to add improvements and to implement variations without writing an extensive theory of the specic data-structures involved. We suggest to replace the core of such tactics with procedures that are proven correct using CoqEAL renement framework, and to build a modular methodology around it. This renement framework addresses the problem of duplication by promoting the use of one extensive proof-oriented library together with one or several more ecient implementations, with a reduced amount of proofs, but destined to computation and proven correct with regard to the proof-oriented library. We show on the example of the ring tactic of Coq that this gain in exibility opens the door to dierent improvements. This is a presentation to trigger discussion about ideas for a prototype based on the existing but improved CoqEAL renement framework, described in [5] and [3].

Introduction

One of the interests of proof-assistants is to make proofs more reliable through systematic checking.

They require the user to provide the steps of a proof and they automate the verication of their validity.

However, having to write every single step makes the proof process quite tedious. Fortunately, most proof-assistants embed automatic procedures that help the user to build proofs. Some of them rely on a technique called reection ( [START_REF] Boutin | Using Reection to Build Ecient and Certied Decision Procedures[END_REF]).

The purpose of reection is to replace proof steps with computation, thus automating and shortening the proof. The proof-assistant Coq ( [START_REF] Tahar | The Coq proof assistant reference manual[END_REF]) makes an extensive use of computation for such a purpose. For instance, the following term is typable in Coq: eq_refl : 0 + n = n. Coq accepts the proof by reexivity because it is able to use the denition of + to compute the left handside into n. Reection performs more complex reasonings thanks to a translation (called metaication in [START_REF] Boutin | Using Reection to Build Ecient and Certied Decision Procedures[END_REF]) of the goal and the description of rules of computation inside a logical language.

Reection is used to design veried decision procedures such as the ring ( [START_REF] Mahboubi | Proving Equalities in a Commutative Ring Done Right in Coq[END_REF]) and field ([4]) Coq tactics. However, these decision procedures are often designed in a monolithic fashion: they rely on ad-hoc data-structures to which specic transformations are applied. To implement variations and/or improvements, one has to dive into the core of such procedures and change the structures/proofs. This may require the development of an extensive theory of the involved structures.

In this paper we propose a more modular methodology for reection, on the particular case of the ring tactic. This methodology relies on CoqEAL renement framework and makes a clear distinction between the proofs of soundness for decision procedures and the computations they will perform. Section 2 presents briey CoqEAL. Then, Section 3 explains how we built a tactic around it. Finally, Section 4 compares our approach to the historical one.

The CoqEAL renement framework

CoqEAL 1 is built on top of the Mathematical Components library and the SSReflect extension of Coq tactic language ( [START_REF] Gonthier | A Small Scale Reection Extension for the Coq system[END_REF]). It contains ecient data-structures and some optimized algorithms in polynomial and matrix computations, and it provides a framework for changing the representation of objects.

Renement is a term usually used to describe a step-by-step approach in the verication of a program. One starts with an abstract representation of the program and its specication and then renes the program to more concrete representations, nally reaching an actual implementation. Each step is proven correct with respect to the previous one. What we mean by renement here ( [START_REF] Dénès | A Renement-Based Approach to Computational Algebra in Coq[END_REF]) is only a veried step for changing the representation of some data-type or the algorithm that is used, typically going from a proof-oriented version of your object to a computation-oriented one.

CoqEAL's version of renement is of interest to us because of its modularity and its automation ( [START_REF] Cohen | Renements for free![END_REF]). Indeed, thanks to the use of parametricity ( [START_REF] Wadler | Theorems for Free![END_REF]), a reformulation of Reynolds' abstraction theorem ( [START_REF] Reynolds | Types, Abstraction and Parametric Polymorphism[END_REF]), it is possible to get for free soundness results on the goal representation of your object if you have proven them on its initial representation. We extended the work of Cohen et al. ( [START_REF] Cohen | Renements for free![END_REF]) with a better automation of the use of parametricity, thanks to Keller and Lasson's plugin ( [START_REF] Keller | Parametricity in an Impredicative Sort[END_REF]) generating the needed instantiations of the parametricity theorem.

Moreover, we also use Coq type class mechanism ( [START_REF] Sozeau | First-Class Type Classes[END_REF]) in CoqEAL to keep a data-base of renement theorems. We use a type class refines, where refines P C stands for the computationoriented object C is a renement of the proof-oriented object P. Together with particular instances dening rules to guide the inference, this gives us a logical program computing renements. Thanks to this logical program, expressions are decomposed into smaller parts for which renement theorems are needed. This means that each specic operation can be considered independently from the others, thus simplifying the proof process. For instance, if we want to get a renement for the polynomial expression X + Y -(1 * Y), it is sucient to give respective renements X', Y', 1', +', -' and *' for X, Y, 1, +,and *.

We will use renement to normalize polynomials by computation. It is thus interesting to be able to get back proof-oriented polynomials after normalization. For this, we use a specic renement spec of the identity function, together with the following lemma:

Lemma refines_spec P C : refines P C -> P = spec C.

When given P (e.g. X + Y -(1 * Y)), rewriting with this lemma will trigger type class inference, thus inferring C (e.g. X' +' Y' -' (1' *' Y')) together with a proof of refines P C, and replace P with spec C. Computation will then result in a proof-oriented object P' replacing P (e.g. X). We wrote a tactic, coqeal_simpl, which goes through these two steps, renement and computation, using simpl as strategy of evaluation (see Figure 1). The prototype we implemented performs polycation by rst computing an abstract syntax tree (AST) together with a variable map, in order to interpret the variables into terms, and then turning it into a Mathematical Components polynomial. This rst step of polycation, called metaication ( [START_REF] Boutin | Using Reection to Build Ecient and Certied Decision Procedures[END_REF]), is performed by recursion on the head of the ring expression.

The tactic recognizes terms generated by the following grammar e ::= 0 j 1 j e 1 + e 2 j e 1 £ e 2 j e 1 e 2 j e j e n j c ; where n ranges over natural numbers and c over integers, and terms that do not t this grammar are metaied as variables. For The second step of polycation computes a Mathematical Components polynomial from the AST (using ast_to_poly) evaluated (using eval_poly) with the variable map to nish polycation.

In our example, we get the polynomial X + Y -(1 * Y) (with adapted notations). This approach is the fastest implementation of polycation we found. Our structure of AST also comes with an interpretation function to original expressions (ast_to_expr), so the soundness of polycation is justied by the equality of the original expression and the evaluation of the resulting polynomial:

Lemma polyficationP env ast : ast_to_expr env ast = eval_poly env (ast_to_poly ast).

In the second step of our tactic, we simply call CoqEAL simplication tactic on the polynomial obtained by polycation, resulting in the previous example in the polynomial X. The last step, depolycation, is implemented for now as a set of rewriting rules to push Mathematical Components polynomial evaluation through the polynomials and get back an expression in the original ring, for instance a in our example. The fundamental steps are the construction of the AST and its evaluation as a sparse Horner normal form. Since, these two datatypes are ad-hoc, any modication of metaication, hence of the AST, would require a new proof of soundness. Adapting proofs is not much of an issue, but doing so on computation-oriented data-structures implies writing the theory of these structures, which signicantly increases the amount of proofs.

Our approach has the following benets. First, one can reuse the results of libraries such as Mathematical Components to prove soundness. This reduces the amount of lemmas to be proven.

Moreover, such libraries usually contain structures that are well-suited for proofs, which means that you don't have to bother with the subtleties of your specic data-structures. Then, as we saw in Section 2, in this framework the soundness of each computable operation can be proven independently from the others. Thus, not only soundness is easier to prove but adding/removing/changing an operation has less impact in terms of proofs than with the former methodology.

As a consequence, our proof of soundness of polyfication (implemented via metaication) is very short (less than 300 lines), and computations came for free, because a renement was already implemented in CoqEAL.

Presently, this implementation is incomplete. To be as powerful as the ring tactic, our prototype lacks some renements for coqeal_simpl, like the power function and some coercions, as well as the possibility to declare translations for ring-specic constants or to use other rings than the set of integers as the ring of coecients for the polynomials. Also, because of the structure of iterated univariate polynomials, the lemma polyficationP holds only for commutative rings.

Concerning eciency, we cannot yet compare to the existing because CoqEAL still suers some problems of eciency in the renement of nested data-structures (polynomials over polynomials over : : : over integers). Because of this, the bottleneck of our current prototype is the renement of polynomial operations which looks exponential in the number of indeterminates, whereas it could be made linear using a better proof-search strategy, which implementation is work in progress. As of today, depolycation is instantaneous on our small examples. But we expect that for bigger terms it could become slower. We suggest to replace this step either by the same process as polycation (going through an AST to replace rewriting by computation), or by a renement of polynomial evaluation.

Future work

To improve the proof power of our tactic, the rst step is to catch up with the ring tactic by adding the aforementioned missing parts and using a better mathematical representation of polynomials, for instance multinomials 3 ([1]), provided there is a renement for them 4 . Then, a slight modication of polycation would allow us to handle morphisms, making it possible to translate the expression f(x + y) -f(y) into X + Y -Y with the variable map [f(x); f(y)] and thus to simplify it into f(x). Another possible improvement would be the use of Gröbner bases 5 to reason modulo equations. The current ring tactic already deals with hypotheses of the form m = p where m is a monomial and p a polynomial (after metaication), but thanks to Gröbner bases m could be any polynomial. This would require a modication of polycation, and the algorithms on Gröbner bases would be proven thanks to the renement framework.

Finally, we could bring exibility to the reduction strategy. We could easily replace simpl in coqeal_simpl with compute, vm_compute or native_compute, with some locking. Moreover, instead of using the coqeal_simpl reduction strategy, we suggest the user could dene his own transformation, like root nding or factoring, which would make it possible to go from an equation of the form

Conclusion

We refactorized a widespread reection methodology ( [START_REF] Boutin | Using Reection to Build Ecient and Certied Decision Procedures[END_REF]) in order to introduce renement in the core of reection-based tactics. We think that they could benet in dierent ways from the use of renement. The major perk of this framework is the possibility to reuse well-developped libraries to prove the soundness of a reection-based procedure while keeping the eciency of structures that are specic to the problematic. Thanks to the automation of renement ( [START_REF] Cohen | Renements for free![END_REF]), that we have improved, the proof process is also simplied.

In the case of our procedure for reasoning using ring axioms, we turned a syntactic step, metaication, into the semantic translation we call polycation. This results in a better distinction between the translation of the object, where only soundness of the encoding has to be proven, and the computation step that follows, which is proven sound independently. While working on polycation, we also learnt the importance of replacing systematic rewriting with computation (i.e. in our case by going through an AST), using rewriting only in the proofs of theorems about functions performing the computation (as in the proof of polyficationP). As illustrated by the example of the ring tactic, we believe that this access to existing libraries and the gain in modularity in the proof process should make it easier to implement and prove variations of decision procedures. We plan to study the possibility to extend this methodology to other decision procedures (such as field, lra, etc...).

Figure 1 :Figure 2 :

 12 Figure 1: The coqeal_simpl tactic

  example, the metaication of a + b -(1 * b) gives the variable map [a; b] with the AST Plus (Var 0) (Plus (Var 1) (Opp (Mul (Const 1) (Var 1)))).

Figure 3 :

 3 Figure 3: The ring tactic

x 2 + 2 c

 22 2bx + c = 0 to the conjunction of b 2 c > 0 and x = b ¦ p b . One can even imagine plugging here an external tool which produces a proof witness.

Acknowledgements

We thank the anonymous reviewers for their useful feedback. 3 https://github.com/math-comp/multinomials 4 Work in progress by Pierre Roux and Erik Martin-Dorel: https://sourcesup.renater.fr/plugins/scmgit/ cgip-bin/gitweb.cgi?p=validsdp.git;a=blob_plain;f=theories/multipoly.v 5 https://github.com/thery/grobner/


Robert George
email: rgeorge@onera.fr

Mikaël Orain

Laurent Zimmer

Caractérisation d'un claquage induit par laser en mono et double-impulsion

Dans les foyers aéronautiques, l'allumage est actuellement réalisé à l'aide de bougies à arc, une méthode robuste mais qui possède de nombreux inconvénients. Une alternative serait de remplacer la bougie par un système d'allumage par laser. Afin de mieux comprendre les mécanismes intervenant dans le claquage laser et l'allumage du carburant, nous avons mis en place un montage multi-diagnostic pour l'étude du claquage dans l'air et dans un mélange kérosène/air. Les différents paramètres étudiés sont l'évolution spatio-temporelle du plasma et de l'onde de choc qu'il initie, le taux d'absorption de l'énergie incidente, ainsi que la spectroscopie d'émission. Une configuration en double-impulsion a également été étudiée.

Introduction

Contexte

Dans les foyers aéronautiques, l'allumage est actuellement réalisé à l'aide de bougies à arc. Cette méthode est techniquement robuste mais possède néanmoins de nombreux inconvénients. Pour l'allumage ou le ré-allumage en altitude, une alternative serait de remplacer la bougie par un système d'allumage par laser basée sur la focalisation d'un faisceau laser pour initier la combustion. Ce procédé possède plusieurs avantages par rapport à la technologie existante utilisant une bougie électrique, notamment le positionnement précis et aisément déplaçable du dépôt d'énergie, le déclenchement précis temporellement, un dépôt d'énergie plus efficacement réparti tant spatialement qu'en termes de transfert d'énergie au mélange carburant / air, un délai d'inflammation fortement réduit, une énergie requise pour l'allumage moindre, etc [START_REF] Phuoc | Laser-induced spark ignition fundamental and applications[END_REF]. Toutefois, ce procédé est très contraint par la technologie laser existante, peu compétitive en termes de coût, robustesse et compacité par rapport à la technologie de la bougie à arc.

Objectifs de l'étude

Afin de mieux comprendre les mécanismes intervenant dans le claquage laser et l'allumage du carburant, il est nécessaire de caractériser proprement la création et l'évolution du plasma, notamment en ce qui concerne le dépôt d'énergie. Nous avons donc mis en place un montage multi-diagnostic pour l'étude du claquage dans l'air et dans un mélange kérosène/air. Les configurations simple et double impulsion ont été étudiées.

Procédure expérimentale

Création du plasma

Le plasma est créé par focalisation d'un faisceau issu d'un laser Nd:YAG (Figure 1). Celui-ci produit des impulsions d'une durée égale à 6 ns (FWHM) à une longueur d'onde de 532 nm. Son régime de fonctionnement est à 10 Hz et l'énergie par impulsion peut aller jusqu'à plusieurs centaines de milli-joules. La focalisation s'effectue grâce à une lentille plan-convexe de 100 mm de focale.

Le claquage a été étudié dans deux configurations différentes :

• l'air du laboratoire • un pré-mélange air/kérosène gazeux chauffé à 180°C Pour la configuration double-impulsion, un second faisceau laser produit par le même laser permet d'obtenir des caractéristiques identiques entre les deux impulsions. Le décalage inter-impulsion est réglé finement grâce à un générateur de délai.

Diagnostics optiques

De multiples diagnostics optiques ont été implémentés afin de caractériser le phénomène de claquage (Figure 1) :

• L'énergie laser, modifiée en ajustant le délai entre le déclenchement des lampes flash et du Q-switch, est mesurée par deux photo-détecteurs placés en amont et en aval de la position du claquage. Ceci nous permet de mesurer l'énergie incidente et l'énergie transmise pour ainsi en déduire la proportion d'énergie absorbée par le plasma.

• Une caméra ICCD (512 x 512 pixels) a été utilisée pour visualiser l'émission du plasma.

Cette dernière a été placée perpendiculairement à l'axe laser. La combinaison d'un filtre interférentiel centré à 312 nm avec une largeur à mi-hauteur de 10 nm et d'un filtre passehaut (WG320, Schott) a permis de collecter l'émission du plasma dans une gamme spectrale très fine, ainsi que de ne pas saturer notre système de détection.

• Un montage de strioscopie a permis la visualisation de l'onde de choc et la mesure de sa propagation. L'éclairement est fourni par une source laser pulsée à 690 nm, fournissant des impulsions de 200 ns. La détection a été effectuée par une caméra commerciale possédant une matrice CCD de 1024 x 0280 pixels.

• Des spectres d'émission ont été enregistrés grâce à l'utilisation d'un spectromètre couplé à une caméra ICCD (512 x 512 pixels). Ces spectres pourront ensuite nous permettre de déterminer différentes caractéristiques du plasma (espèces présentes, richesse, température, etc.).

Pour l'enregistrement des spectres et images obtenues par les caméras ICCD, le gain ainsi que le temps d'exposition ont été ajustés pour obtenir un bon rapport signal-sur-bruit. Un générateur de délai a été utilisé pour déclencher les impulsions laser ainsi que les acquisitions de notre système de visualisation, de nos photo-détecteurs et de notre spectromètre. Nous avons ainsi pu synchroniser nos diagnostics ce qui nous a permis de mesurer différents paramètres pour un même évènement. Des mesures ont été effectuées entre 0 et 10 µs après le claquage laser.

Résultats

Mesures du dépôt d'énergie dans le plasma

A cause des processus physiques qui interviennent dans la création d'un plasma par claquage laser, il existe un seuil en dessous duquel aucun plasma n'est créé. Dans le cas de notre montage, on observe la création d'un plasma pour des énergies supérieures à 19 mJ, ce qui correspond au seuil de claquage de notre montage (Figure 2). On observe un taux d'absorption croissant de 20% jusqu'à 85% pour des énergies laser supérieures à 150 mJ. On remarquera une dispersion importante sur l'énergie absorbée qui est due aux processus stochastiques du claquage. Dans la configuration en double-impulsion, les images de l'émission du plasma (Figure 6) nous ont permis de mettre en évidence une particularité dans l'absorption de l'énergie de la seconde impulsion par le plasma. En effet, la seconde impulsion n'est pas directement absorbée par le premier, mais créé deux plasmas à ses extrémités. La Figure 7.a montre la propagation d'une onde de choc initiée par un claquage provenant d'une impulsion de 170 mJ. A partir de ces images, un algorithme de reconnaissance de cercles a permis de déterminer les évolutions temporelles du rayon de l'onde de choc (Figure 7.b) et de sa vitesse. Le rayon augmente très rapidement aux temps courts (t < 1 µs), puis croît quasilinéairement au-delà de 5 µs. On observe également le détachement de l'onde de choc par rapport au plasma aux environs de t= 1,5 µs.

Simulation de la propagation

La modélisation de la propagation d'une onde de choc a fait l'objet de nombreux travaux. Le premier modèle, et le plus connu, est le modèle de Taylor [START_REF] Taylor | The formation of a blast wave by a very intense explosion[END_REF]. Cependant, ce modèle suppose un choc dit « fort » (M>2), ce qui n'est le cas dans notre situation qu'aux tout premiers instants. Jones a ensuite développé un modèle pour les chocs « intermédiaires » [START_REF] Jones | Intermediate strength blast waves[END_REF] que nous avons pu appliquer à notre situation. Ce modèle analytique prévoit une évolution du rayon de l'onde de choc selon l'équation :

(1)

= 0,454.

.

1 + 1,542. . . - 1
Où E 0 est l'énergie de l'onde de choc (J), c 0 la vitesse du son dans le milieu (m/s), P 0 la pression du milieu (Pa) et t le temps (s). En comparant nos mesures aux simulations, il nous a donc été possible de déterminer l'énergie E 0 de l'onde de choc pour différentes énergies incidentes (Tableau 1).

Energie incidente (mJ)

33,5 +/-0,3 63,9 +/-0,5 Tableau 1 : Mesures de l'énergie dissipée par l'onde de choc pour plusieurs énergies incidentes

Spectroscopie d'émission du plasma

La spectroscopie d'émission est un diagnostic optique consistant à créer un plasma par claquage laser dans un environnement que l'on souhaite caractériser, et à en étudier l'émission pour quantifier divers paramètres (espèces présentes, concentrations, etc.). Dans le cadre de notre étude, la spectroscopie d'émission peut non seulement nous aider à analyser le plasma, mais également le milieu dans lequel l'allumage aura lieu.

Evolution temporelle du spectre d'émission

Les mesures spectroscopiques nous ont permis de distinguer deux régimes d'émission : une émission ionique majoritaire jusqu'à 400 ns après le claquage, suivie d'une émission atomique majoritaire. Dans l'air, l'émission atomique provient uniquement de l'oxygène (777, 822 et 845 nm) et de l'azote (744 et 870 nm). En présence d'un mélange air / kérosène gazeux (Figure 6), on observe également les raies de l'hydrogène atomique (486 et 656 nm). Pour une espèce atomique donnée dont la distribution de population des niveaux d'énergies suit une distribution de Boltzmann (hypothèse de l'équilibre thermodynamique local), il est possible de déterminer la température électronique du milieu à partir du rapport de l'émission de deux raies à l'aide de l'équation :

(2)

) ln( 1 2 2 2 2 1 1 1 2 1 I g A I g A k E E T B ν ν - =
Où k B est la constante de Boltzmann, E n l'énergie du niveau considéré (cm -1 ), A n sa probabilité de transition (s -1 ), g n son poids statistique, ν n sa fréquence de raie (Hz) et I n l'intensité de la raie correspondante (avec n=1 ou 2).

Le choix du couple de raies pour la mesure de température s'est porté sur les raies de l'oxygène à 777 et 795 nm. Celui-ci nous offre un bon compromis avec de fortes intensités de raies, ainsi qu'une différence de leurs niveaux d'énergie E n -E 2 nous permettant d'obtenir une bonne précision sur la détermination de la température. On note que la température du plasma 200 ns après l'impulsion laser est similaire (~19000K) quelle que soit l'énergie laser incidente. En revanche, on note une décroissance plus marquée pour des énergies incidentes décroissantes (Figure 9). Lorsque le claquage laser se produit dans un milieu comportant des hydrocarbures, l'émission des raies atomiques de l'hydrogène est directement liée à la concentration du milieu en hydrocarbure, et donc à la richesse locale au point où se créé le plasma. Une étude spectroscopique peut donc permettre la mesure de cette richesse, lorsque l'intensité de la raie est hydrogène est normalisée par la mesure d'une seconde raie [START_REF] Ferioli | Measurements of Hydrocarbons using LIBS[END_REF].

Conclusion

De multiples diagnostics optiques ont été mis en place afin de caractériser le phénomène de claquage laser, le dépôt d'énergie dans le milieu ainsi que l'évolution et la propagation du plasma résultant ainsi que celle de l'onde de choc créée. Nous avons ainsi pu quantifier l'énergie absorbée par le plasma en fonction des paramètres laser, ainsi que la part dissipée dans l'onde de choc. L'étude de la spectroscopie d'émission du plasma fût utilisée pour mesurer l'évolution temporelle de la température, et pourrait permettre des mesures de richesse en présence d'un mélange combustible.

La comparaison entre simple et double-impulsion a également permis de mettre en évidence un dépôt d'énergie et une expansion du plasma bien plus singulière avec cette seconde configuration, le dépôt d'énergie de la seconde impulsion étant répartie aux extrémités du premier plasma.

Figure 1 :

 1 Figure 1 : Schéma du banc de caractérisation

Figure 2 :

 2 Figure 2 : Taux d'absorption du faisceau laser en fonction de l'énergie incidente La Figure 3 montre que, dans le cas d'un claquage en double-impulsion, le taux d'absorption de la deuxième impulsion dépend aussi bien de l'énergie de cette impulsion, mais également du délai entre la première et la deuxième impulsion. En effet, pour des énergies de la seconde impulsion faibles (29 mJ) et intermédiaires (91 mJ), on observe l'apparition d'un domaine temporel durant lequel l'absorption chute fortement. Néanmoins, pour des énergies de la seconde impulsion croissantes, cet effet devient de moins en moins prépondérant, jusqu'à ce que l'absorption reste quasi-inchangée pour des énergies élevées).

Figure 3 :

 3 Figure 3 : Taux d'absorption de la 2 ème impulsion laser en fonction de son énergie laser incidente et du délai inter-impulsions 3.2 Visualisation du plasma La Figure 4 montre l'évolution spatiale du plasma à différents instants après le tir laser. Les images correspondent à une énergie incidente de 171 mJ. Le plasma est de forme allongée dans la direction du faisceau laser, sa forme et sa répartition spatiale d'intensité varient au cours du temps.La taille du plasma commence par augmenter, puis décroît 3 µs après le tir laser (figure5). On constate que la longueur est maximale vers 2 µs pour les trois énergies laser considérées. Comme on pouvait le supposer, la taille du plasma augmente également avec l'énergie laser incidente.

Figure 4 :Figure 5 :

 45 Figure 4 : Evolution spatio-temporelle de l'émission du plasma en simple impulsion pour une énergie incidente de 170 mJ

Figure 6 : 3 . 3

 633 Figure 6 : Evolution spatio-temporelle de l'émission du plasma en double impulsion (E 1 = 30 mJ ; E 2 = 140 mJ ; délai inter-impulsions = 300 ns) 3.3 Visualisation de l'onde de choc par strioscopie 3.3.1 Propagation de l'onde de choc La création et la propagation d'une onde de choc et l'un des modes de déperdition par lequel l'énergie absorbée par le plasma est dissipée [2]. Cette énergie mécanique ne peut cependant pas être considérée comme une perte pure dans le cadre d'un allumage car l'onde de choc impacte l'environnement immédiat, et peut affecter la propagation du noyau de flamme. Malgré cela, il reste intéressant de déterminer quelle proportion de l'énergie initiale est dissipée par ce processus.

Figure 8 :

 8 Figure 8 : Spectre d'émission du plasma crée par une impulsion laser de 170 mJ dans un mélange air/kérosène de richesse 0,612 3.4.2 Mesures de température par l'étude du spectre d'émission

Figure 9 :

 9 Figure 9 : Evolution temporelle de la température du plasma pour plusieurs énergies laser incidentes

Figure 10 :

 10 Figure 10 : Ratio de l'intensité des raies H α et O 777 pour différentes valeurs de richesse Sur la Figure 10, on peut observer une quasi-linéarité dans la relation entre richesse et ratio H/O jusqu'à une richesse de 0,7. Les points au-delà n'ont pas été considérés car le claquage était devenu instable.

Références