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 calculusof-variations, we show that the resulting action 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 multi-scale 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 a type of stochastic hybrid system (SHS) known as a piecewise deterministic Markov process (PDMP) [10] ‡. One important example at the single-cell level is the occurrence of 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] Newby | Spontaneous excitability in the Morris-Lecar model with ion channel noise[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 variable represents the activation state of the gene [23,[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 relaxation dynamics then t/ is the characteristic time-scale of the Markov chain 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). This is a consequence of the fact that the paths of escape are constant energy solutions that converge to a stable fixed point of the underlying deterministic system in the limit t → -∞, and the Hamiltonian is zero at fixed points. 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-time 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][START_REF] Bressloff | Path-integral methods for analyzing the effects of fluctuations in stochastic hybrid neural networks[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 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 classical action arising from 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]. Using the Perron-Frobenius theorem and the calculus-of-variations, we evaluate the LDP rate function in terms of the classical action, whose equations of motion along the most probable paths of escape are given by a Hamiltonian dynamical system. We show that one major difference between hybrid and non-hybrid stochastic processes is that in the former case the WKB Hamiltonian tends to be a reduced version of the LDP Hamiltonian. However, the corresponding Hamiltonian dynamical systems yield the same most probable paths of escape.

The structure of the paper is as follows. In section 2 we define a stochastic hybrid system corresponding to a one-dimensional PDMP, and specify our various mathematical assumptions. We then present our detailed derivation of the classical LDP action in section 3 (theorem 3.1 and its proof). Although we focus on a onedimensional PDMP, we also indicate how to extend our results to higher-dimensional PDMPs and multi-scale processes. The relationship between the classical action of large deviation theory and the quasipotential of WKB theory is developed in section 4. In particular, we show how the WKB quasipotential satisfies a Hamilton-Jacobi equation, whose associated Hamiltonian is consistent with the Hamiltonian obtained from the classical LDP action. 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]. Finally, we illustrate our analysis of LDPs for stochastic hybrid systems by considering 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 x ∈ Ω ⊂ R and a discrete variable n ∈ Γ ≡ {0, • • • , N } [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) ẋ = F n (x), (2.1)
where the vector field F n : Ω → R is a continuous function, locally Lipschitz. That is, for any 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 ∈ K. (2.2)
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), see section 3.5, 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. § 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].

Let us decompose the transition matrix of the Markov chain as W nm (x) = P nm (x)λ n (x), with m,m =n P nm (x) = 1 for all x. That is, for a given x, the jump times from state n are exponentially distributed with rate λ n (x) and P nm (x) is the probability distribution that when a jump occurs the new state is m for some m = n. 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 .
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 . 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ρ, (2.5)
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 (x) = W nm (x) -δ n,m k∈Γ W nk (x), (2.7)
such that m∈Γ A nm (x) = 0 ∀n ∈ Γ and ∀x ∈ Ω, 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 with the piecewise deterministic dynamics for a given n, whereas the second term represents jumps in the discrete state. Note that we have rescaled the transition matrix W (and hence 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. 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 action from a large deviation principle

We now turn to the major part of the paper, namely, a new derivation of the classical action corresponding to the abstract rate function of the LDP for a one-dimensional PDMP 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]. We begin by introducing some notation in order to state 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] (section 3.1). The associated rate function J T is expressed as a variational principle. In section 3.2 we summarize the main result of our paper in the form of theorem 3.1, which expresses the equivalence of the rate function J T with a classical action, whose associated Lagrangian is related to the Perron eigenvalue of a certain linear operator. The latter depends on the transition matrix W of the Markov chain and the nonlinear functions F n (x), n ∈ Γ of the piecewise deterministic system (2.1). The proof of theorem 3.1 is presented in section 3.3, where we explicitly derive the classical action by solving the variational principle using the Perron-Frobenius theorem and the calculus-of-variations. Extensions of the classical action to higher dimensional PDMPs and multi-scale processes are presented in sections 3.4 and 3.5, respectively.

3.1. Large deviation principle 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] 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 + (Γ) denote the space of all probability measures on Γ and introduce the product space

M + (Γ) [0,T ] . For each t ∈ [0, T ], let ψ(t) = (ψ 0 (t), . . . , ψ N (t)) ∈ M + (Γ) 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 ] 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 T ] with topology defined by the metric [START_REF] Faggionato | Averaging and large deviation principles for fully-coupled piecewise deterministic Markov processes and applications to molecular motors[END_REF]

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
ψ 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,
d({(x(t), ψ(t))} t∈[0,T ] , {( x(t), ψ(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 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 z∈(0,∞) Γ (n,m)∈Γ×Γ ψ 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. Let R N denote
the set of (N + 1) × (N + 1) positive matrices with diagonal zero. Define the mapping J : R N → R as

J (c) = sup z∈[0,∞) Γ J (c, z), J (c, z) = (n,m)∈Γ∆ c nm (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 J (c, •) is the unique solution of the set of equations

m∈Γ c nm z m z n = m∈Γ c mn 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 nonvanishing 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 [23,[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].)

Statement of main theorem

We now state the main result of our paper in the form of a theorem, which will then be proven in section 3.3. Theorem 3.1 Let λ = λ(x, µ) for fixed x, µ be 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.8)
where z(x, µ) = (z n (x, µ), n ∈ Γ) is the unique positive eigenvector of the linear equation under the normalization

n z n (x, µ) = 1. Similarly, let R(x, µ) = (R n (x, µ), n ∈ Γ) be the positive eigenvector of the adjoint equation m∈Γ A nm (x)R m (x, µ) + µF n (x)R n (x, µ) = λ(x, µ)R n (x, µ) (3.9) under the normalization m z m (x, µ)R m (x, µ) = 1.
The rate function J T ({x(t)} t∈[0,T ] ) defined in equation (3.4) can then be written in the form of the classical action

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

L(x, ẋ) = µ(x, ẋ) ẋ -λ(x(t), µ(x, ẋ)), (3.11)
and µ = µ(x, ẋ) is given by the unique solution of the invertible equation

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

ψ m (x, µ) = z m (x, µ)R m (x, µ). (3.13) 3.3. Proof of Theorem 3.1
The proof proceeds in two main steps.

[A] Evaluating the supremum of equation (3.5). We begin 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.14)
for some bounded vector q = (q n , n ∈ Γ). Here, for any

a, b ∈ R N +1 , [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 N +1 ) = λ(x, q) -c. (3.15)
For simplicity, we choose c = q N and take Q = (q 0 , q 1 , . . . , q N -1 , 0) so that λ(x, Q) = λ(x, q) -q N and z = z(x, Q) are solutions of the matrix equation

A(x)z(x, Q) + Q • z(x, Q) = λ(x, Q)z(x, Q). (3.16)
There are N independent variables, q n , n = 0, . . . , N -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=0,...,N -1

{|q n |}, (3.17)
it is clear that the new matrix A(x) + diag(Q) + κI N +1 is irreducible and positive.

According to the Perron-Frobenius theorem, it has a unique strictly positive eigenvector z(x, Q) with the normalization m z m (x, Q) = 1, and this eigenvector is also an eigenvalue of A + diag(Q) (but with a shifted eigenvalue). We now show that for a given Q and fixed x, the Perron eigenvector z(x, Q) is the solution to the variational problem (3.5) provided that the measure ψ takes the specific form

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

A (x)R(x, Q) + Q • R(x, Q) = λ(x, Q)R(x, Q), (3.19)
such that z R = 1. Equation (3.18) ensures that ψ is a strictly positive measure and that N m=0 ψ m = 1. We proceed by establishing that z(x, Q) satisfies equation (3.7) for c nm = ψ n (x, Q)W nm (x) and ψ n (x, Q) given by equation (3.18). The left-hand side of (3.7) becomes

N m=0 ψ n W nm (x) z m z n = ψ n z n N m=0 W nm (x)z m = ψ n z n λ(x, Q) -Q n + k W nk z n = ψ n λ(x, Q) -Q n + k W nk , (3.20)
whereas the right-hand side of (3.7) reads

z n N m=0 W mn (x) ψ m z m = z n N m=0 A mn (x) + δ mn N k=0 W nk ψ m z m .
Cancelling the factors of z n , setting ψ m /z m = R m and using equation (3.19) we recover (3.20). Hence, equation (3.7) holds. Finally, setting

ψ n = ψ n (x, Q) in equation (3.5)
we can evaluate the supremum to obtain

j(x, ψ) = N m,n=0 ψ n W nm (x) 1 - z m z n = N n=0 ψ n N m=0 W nm (x) - 1 z n N m=0 A nm (x) + δ nm N k=0 W nk (x) z m = N n=0 ψ n [Q n -λ(x, Q)] = N -1 n=0 q n ψ n -λ(x, Q), (3.21)
since N m=0 ψ m = 1 and Q N = 0. For ease of notation, we have suppressed the explicit dependence of z and ψ on x, Q.

In the given variational problem for fixed x, there are N independent variables ψ n , n = 0, . . . , N -1 with ψ N = 1 -N -1 n=0 ψ n . Similarly, there are N independent variables q n , n = 0, . . . , N -1. Therefore, equation (3.18) determines a mapping between the sets {q n , n = 0, . . . , N -1} and {ψ n , n = 0, . . . , N -1}. It remains to show that there exists a unique solution Q for each ψ ∈ M + (Γ), that is, the mapping is invertible. For then we can set q n = q n (x, ψ) in equation (3.21) to obtain the value of the supremum in equation (3.5) for any positive measure ψ ∈ M + (Γ). Differentiating equation (3.19) with respect to q m , m = 0, . . . , N -1, yields the inhomogeneous linear equation

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

∂λ(x, Q) ∂q m = R m (x, Q)z m (x, Q), m = 0, . . . , N -1. (3.23)
Since R and z are strictly positive, λ(x, Q) is a monotonically increasing function of the q m . Moreover, equations (3.14) and (3.19) 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 + (Γ) there exists a vector Q such that ψ n = R n (x, Q)z n (x, Q) for all n = 0, . . . , N -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.23) 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.24)
for all m, n = 0, . . . , N -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.25)
with n, m = m . On the other hand, differentiating equation (3.21) for j(x, ψ) with respect to ψ n gives ∂j(x, ψ)

∂ψ n = q n + N -1 j=0 ψ 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 = 0, . . . , N -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.16). Now suppose ψ is a nonnegative rather than a strictly positive 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.21) 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 N = 1. 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 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 → ±∞.

[B] Evaluating the infimum of equation (3.4). The next step in the proof is to substitute for j(x, ψ) in the rate function (3.4) using equation (3.21), which gives

J T ({x(t)} t∈[0,T ] ) = inf ψ: ẋ= N n=0 ψnFn(x) T 0 N -1 n=0 q n (t)ψ n (t) -λ(x(t), Q(t)) dt, (3.26)
with q n (t) = q n (x(t), ψ(t)) and N m=0 ψ 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.27)
where

S[x, µ, ψ] = T 0 N -1 n=0 q n (t)ψ n (t) -λ(x(t), Q(t)) (3.28) +µ(t) ẋ - N -1 n=0 [F n (x) -F N (x)]ψ n (t) -F N (x) dt,
and we have imposed the constraint

q m = µ[F m (x) -F N (x)], (3.30a
)

ψ m (x, µ) = z m (x, µ)R m (x, µ), (3.30b)
for all m = 0, . . . N -1, with R n (x, µ), z n (x, µ) the positive eigenvectors of the matrix equation

A (x)R(x, µ) + µF(x) • R(x, µ) = λ(x, µ)R(x, µ), (3.31)
and its adjoint, respectively. We have thus established Eqs. (3.8) and (3.9). Here

λ = λ(x, µ) ≡ F N (x) + λ(x, Q|q m = µ(F m -F N ), m = 0, . . . , N -1).
Substituting equation (3.30b) into (3.29a) shows that µ = µ(t) where µ(t) is the solution to the equation

ẋ(t) = N n=0 F n (x(t))ψ n (x(t), µ(t)), (3.32)
provided that the equation is invertible for the given trajectory {x(t)} t∈[0,T ] . Hence, evaluating the action at the infimum we have

J T ({x(t)} t∈[0,T ] ) = T 0 µ(t) N n=0 F n (x(t))ψ n (µ(t), x(t)) -λ(x(t), µ(t)) dt,
The final step is to show that equation (3.32) is invertible so that the function µ(t) = µ(x(t), ẋ(t)) exists and, hence, the rate function J T has the required Lagrangian form (3.10). From the inverse function theorem we require that

L(x, µ) ∂R ∂µ = ∂λ ∂µ R -F • R, (3.33) with L(x, µ) = L(x, Q)| qm=µ[Fm-F N ],m=0,...,N -1 .
Using the same arguments as previously, we obtain

∂λ(x, µ) ∂µ = N n=0 F n (x)z n (x, µ)R n (x, µ). (3.34)
Hence, we require

∂ 2 λ(x, µ) ∂µ 2 ≡ N -1 m,n=0 ∂ 2 λ(x, Q) ∂q m ∂q n qm=µ(Fm(x)-F N (x)) × (F m (x) -F N (x))(F n (x) -F N (x)) = 0.
This holds since the Jacobian D of equation (3.24) is invertible and F m (x) = F N (x) for at least one m = N . Finally, from equation (3.17), we require µ to be bounded, that is, there exists a κ for which κ > µ max In theorem 3.1 we considered the case of one-dimensional piecewise deterministic dynamics by taking x ∈ R. However, it is relatively straightforward to derive the corresponding classical action for x ∈ R d . When the internal state is n, the system now evolves according to the ODE ẋ = F n (x), (3.35) 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 K n (K) such that

|F n (x) -F n (y)| ≤ K n (K)|x -y|, ∀x, y ∈ K. (3.36)
The rate function of the LDP (3.3) becomes

J T ({x(t)} t∈[0,T ] ) = T 0 L(x, ẋ)dt, (3.37)
where L is the Lagrangian L(x, ẋ) = µ(x, ẋ), ẋ -λ(x(t), µ(x, ẋ)), and λ(x, µ) is the Perron eigenvalue of the linear equation

A(x)z(x, µ) + (F(x)µ) • z(x, µ) = λ(x, µ)z(x, µ), (3.38)
where z(x, µ) is the positive eigenvector of the linear equation under the normalization

n z n (x, µ) = 1. Here [a • b] n = a n b n , n = 0, • • • , N (3.39)
for any a, b ∈ R N +1 and F(x) is the (N + 1) × d matrix whose N + 1 rows are the d-dimensional vectors F m (x), m = 0, • • • , N . Finally, the d-dimensional vector µ = µ(x, ẋ) is the solution of the invertible equation

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

ψ m (x, µ) = z m (x, µ)R m (x, µ), (3.41)
where R is the positive eigenvector of the adjoint equation

A (x)R(x, µ) + (F(x)µ) • R(x, µ) = λ(x, µ)R(x, µ) (3.42)
under the normalization m z m (x, µ)R m (x, µ) = 1.

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 infimum 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.43)
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 piecewise SDE is coupled to the fast discrete process on Γ, the stochastic dynamics is described by a differential Chapman-Kolmogorov equation of the form (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.44)
where P{X(t) ∈ (x, x + dx), n(t) = n|x 0 , n 0 ) = ρ(x, n, t)dx. We define the measure space M + (Γ) 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) = N n=0 ψ n (t)F n (X)dt + √ N n=0 ψ n (t)σ 2 n (X) dW (t) ≡ F (X, ψ)dt + √ σ(X, ψ)dW (t) (3.45) for ψ ∈ M + (Γ) [0,T ]
. Such a space includes the set of trajectories of the piecewise SDE (3.43) 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 of W , one can write down an LDP along identical lines to the case of a piecewise deterministic system, see equation (3.3), yielding

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.46)
where X(t) is the sample path generated by the particular realization of W (t), and j(X, ψ) is given by equation (3.21). If we now formally average with respect to the Wiener process we find that

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

S[x, ψ] = T 0 N -1 n=0 q n (t)ψ n (t) -λ(x(t), Q(t)) + ( ẋ - N n=0 ψ n F n (x)) 2 2 N n=0 ψ n σ 2 n (x) dt, (3.48)
with q n (t) = q n (x(t), ψ(t)), N m=0 ψ m (t) = 1 and λ the Perron eigenvalue of equation (3.16). Note that the final term on the right-hand side is the well-known Freidlin-Wentzell action for SDEs [START_REF] Freidlin | Random perturbations of dynamical systems[END_REF]. Taking the infimum with respect to ψ k , k = 0, . . . , N -1, gives

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

µ = N n=0 ψ n F n (x), σ 2 = N n=0 ψ n σ 2 n (x), ν = ẋ -µ σ 2 , (3.50)
we have

q k = ν[F k (x) -F N (x)] + ν 2 2 [σ 2 k (x) -σ 2 N (x)], (3.51)
and

ψ m = ∂λ/∂q m = R m (x, ν)z m (x, ν), (3.52)
with R n (x, ν), z n (x, ν) the unique positive eigenvectors of the matrix equation

n A mn (x)R n (x, ν) + (νF m (x) + ν 2 σ 2 m (x)/2)R m (x, ν) = λ(x, ν)R m (x, ν). (3.53)
and its adjoint. Moreover, equation (3.50) implies that ν = ν(x(t), ẋ(t)) ≡ ν(t) where ν(t) is the solution to the equation

ẋ(t) = ν(t) N n=0 σ 2 n (x(t))ψ n (x(t), ν(t)) + N n=0 F n (x(t))ψ n (x(t), ν(t)), (3.54)
provided that the equation is invertible for the given trajectory {x(t)} t∈[0,T ] . Finally, evaluating the action at the infimum yields the corresponding classical action:

J T ({x(t)} t∈[0,T ] ) = T 0 ν(t) N n=0 F n (x(t))ψ n (x(t), ν(t)) + ν(t) 2 2 N n=0 ψ n (x(t), ν(t))σ 2 n (x(t)) -λ(x(t), ν(t)) + ν(t) 2 2 dt = T 0 L(x, ẋ)dt. (3.55)

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 a Hamilton-Jacobi equation H(x, ∂ x Φ) = 0, where H is the Hamiltonian obtained from the classical action of large deviation theory. It turns out that such a connection is more subtle in the case of a PDMP, as we now show. We will focus on the one-dimensional case of theorem 3.1. Given the action (3.10), we can determine the Hamiltonian H from the Lagrangian L according to the Fenchel-Legendre transformation:

H(x, p) = sup y [(p -µ(x, y))y + λ(x, µ(x, y))] . (4.1)
Evaluating the right-hand side yields the equation

p -µ(x, y) + ∂λ ∂µ -y ∂µ ∂y = 0 (4.2) with y = m∈Γ ψ m (x, µ)F m (x).
Differentiating equation (3.8) with respect to µ gives

m∈Γ A nm (x) ∂z m (x, µ) ∂µ + [µF n (x) -λ(x, µ)) ∂z n (x, µ) ∂µ = ∂λ(x, µ) ∂µ -F n (x) z n (x, µ). (4.3)
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 (4.2) thus shows that p = µ, and we can identify p as the "conjugate momentum" of the Hamiltonian

H = λ(x, p), (4.4)
where λ(x, p) is the Perron eigenvalue of the linear equation (3.8) and its adjoint (3.9). In order to relate the Hamiltonian (4.4) with the quasipotential of WKB methods, consider the 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.5)
The WKB approximation of ρ ss takes the form

ρ ss (x, n) ∼ η(x, n) exp - Φ(x) . (4.6)
Substituting into equation (4.5) yields m∈Γ

A nm (x) + Φ (x)δ n,m F m (x) η(x, m) = dF n (x)η(x, n) dx , (4.7)
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.8)
As it stands, it is not clear that (4.8) 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.8) is a reduced version of equation (3.9). 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.9)
which are parameterized by the pair (x, p) with p an auxiliary variable and λ(x, p) the Perron eigenvalue. Comparison of equation (4.8) with (4.9) 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. The quasipotential is then the solution of the Hamilton-Jacobi equation

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

and identifying Φ 0 as the action along the resulting solution curve (x(t), p(t)):

Φ 0 (x) = T -∞ p(t) ẋ(t)dt, (4.12)
with x(T ) = x. 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 of section 5. Given the quasipotential Φ 0 , the mean first passage time τ to escape from x -via x 0 , say, can be calculated by considering higher order terms in the WKB approximation, and using matched asymptotics to deal with an imposed 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 quasipotential and its associated Hamiltonian in terms of an underlying variational principle for most probable 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), The stationary measure of the Markov chain is given by

ẋ = F n (x), n = 0, 1, (5.1
ρ(x, 0) = ω -(x) ω -(x) + ω -(x) , ρ(x, 1) = ω + (x) ω -(x) + ω -(x) . (5.2)
The linear equation (3.9) 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 current.

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 2. 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 asymptotics have been used to calculate the MFPT to escape from x -.

Here, we will focus on the quasipotential and its relation to the Perron eigenvalue Substituting the explicit expressions for A and F n (x) into equation (3.9), yields the following equation for the Perron eigenvalue λ and the right eigenvector R:

(N -n + 1)α(x)R n-1 -[λ + nβ(x) + (N -n)α(x)]R n + (n + 1)β(x)R n+1 = -p n N f (x) -g(x) R n . (5.16)
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α

Γ + Γβ(N -n) -λ -nβ -(N -n)α = -p n N f -g .
Collecting terms independent of n and terms linear in n yields the pair of equations Eliminating Γ from these equation gives

p = - N f (x) 1 Γ(x, p) + 1 (α(x) -β(x)Γ(x, p)) , (5.18)
p = 1 f (x) N β(x) λ(x, p) + N α(x) + pg(x) + 1 (λ(x, p) + pg(x)).
This yields a quadratic equation for λ of the form

λ 2 + σ(x)λ -h(x, p) = 0, (5.20)
with

σ(x) = (2g(x) -f (x)) + N (α(x) + β(x)), h(x, p) = p[-N β(x)g(x) + (N α(x) + pg(x))(f (x) -g(x))].
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

h(x, ∂ x Φ 0 ) = 0. (5.21)
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 :

Φ 0 = 0 and Φ 0 (x) = -N α(x)f (x) -(α(x) + β)g(x) g(x)(f (x) -g(x)) . (5.22)
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 .

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 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 that at time t there are n = 0, (n, m ± 1), (5.27) with transition rates ω - 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). (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 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

N m=0 ψ m = 1 .

 1 The variational problem can now be expressed in terms of functional derivatives of S:δS δµ(s) = ẋ(s) -s)[F m (x(s)) -F N (x(s))] = 0 (3.29b)for m = 0, . . . N -1. Combining with equations (3.18) and (3.23), we obtain the following solution to the variational problem in terms of µ:

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

  m=0,...,N {|F m (x)|} for all x ∈ Ω. Hence, we have obtained the classical action (3.10) from the LDP rate function (3.4) and the proof of Theorem 3.1 is complete.

3. 4 .

 4 Extension to higher dimensions, x ∈ R d , d > 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

Figure 1 .

 1 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 2 .

 2 Figure 2. Phase portrait of Hamilton's equations of motion for the ion channel model with Hamiltonian given by the Perron eigenvalue. (x and p are taken to be dimensionless.) The zero energy solution representing the maximum likelihood path of escape from x-is shown as the gray curve. (The corresponding path from x+ is not shown.) Same parameter values as Fig. 1 and I = 0.

26 )

 26 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 toThe 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) → ω ± K (m,x)

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.

[ 7 ]

 7 dX(t) = F n (X, Y )dt (5.29a)dY (t) = [Ω + (X, Y ) -Ω -(X, Y )]dt + σ 2 (X, Y )dW (t),(5.29b)withF n (x, y) = n N f Na (x) + yf K (x) + f L (x) + I, σ 2 (x, y) = Ω + (x, y) + Ω -(x, y). (5.30)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λ(x, p)R n (x, p) = m A nm (x)R m (x, p) + {p x F n (x) + p y [Ω + (x, y) -Ω -(x, y)]} R n (x, p)

Acknowledgements

We thank Jay Newby for drawing our attention to the work of Yuri Kifer [24] 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.

+ 1 2 p 2 y σ 2 (x, y)R n (x, p), (5.31) 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

where g(x, y) = yf K (x) + f L (x) and

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

and

). 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].


Alin Bostan
email: alin.bostan@inria.fr

Teresa Carlos D'andrea

Agnès Krick
email: krick@dm.uba.ar

Marcelo Szanto
email: aszanto@ncsu.edu

Valdettaro
email: mvaldettaro@gmail.com

C D'andrea
email: cdandrea@ub.edu

T Krick

A Szanto

M Valdettaro

Subresultants in multiple roots: an extremal case

Keywords: Subresultants, Hankel matrices, Ostrowski's determinant

come

Introduction

Let K be a field, and f = f m x m + • • • + f 0 and g = g n x n + • • • + g 0 be two polynomials in K[x] with f m = 0 and g n = 0. Set 0 ≤ d < min{m, n}. The order-d subresultant Sres d (f, g) is the polynomial in K[x] defined as Sres d (f, g) := det m+n-2d

f m • • • • • • f d+1-(n-d-1) x n-d-1 f . . . . . . . . . n-d f m . . . f d+1 f g n • • • • • • g d+1-(m-d-1) x m-d-1 g . . . . . . . . . m-d g n • • • g d+1 g , (1)
where, by convention, f ℓ = g ℓ = 0 for ℓ < 0. Although it is not immediately transparent from the definition, Sres d (f, g) is a polynomial of degree at most d, whose coefficients are equal to some minors of the Sylvester matrix of f and g. Subresultants were implicitly introduced by Jacobi [Jac1836] and explicitly by Sylvester [Syl1839, Syl1840], see [GL2003] for a comprehensive historical account1 .

For any finite subsets A = {α 1 , . . . , α m } and B = {β 1 , . . . , β n } of K, and for 0 ≤ p ≤ m, 0 ≤ q ≤ n, one can define after Sylvester [Syl1840b] the double sum expression: Syl p,q (A, B)(x) :=

A ′ ⊂A, B ′ ⊂B |A ′ |=p, |B ′ |=q R(A ′ , B ′ ) R(A\A ′ , B\B ′ ) R(A ′ , A\A ′ ) R(B ′ , B\B ′ ) R(x, A ′ ) R(x, B ′ ),
where R(Y, Z) := y∈Y z∈Z (y -z).

Sylvester stated in [Syl1840b], then proved in [Syl1853, Section II], the following connection between subresultants and double sums: assume that d = p + q, and suppose that f and g are the square-free polynomials

f = (x -α 1 ) • • • (x -α m ) and g = (x -β 1 ) • • • (x -β n ).
Then, d p Sres d (f, g) = (-1) p(m-d) Syl p,q (f, g).

This identity can be regarded as a generalization to subresultants of the famous Poisson formula [Poi1802] for the resultant of f and g:

Res(f, g) = m i=1 n j=1 (α i -β j ).
(2)

We note however that the Poisson formula also holds when f or g have multiple roots, since it does not involve denominators in terms of differences of roots in subsets of A or in subsets of B. To demonstrate the challenges in finding closed formulae for subresultants in the most general case, consider the instance when f

= (x -α 1 ) m 1 • • • (x -α r ) mr , g = (x -β 1 ) n 1 • • • (x -β s ) ns
with α i = α j , β k = β ℓ and d = 1. A (quite intricate) closed formula for Sres 1 (f, g) appears in [DKS2013, Th.2.7] and has the form:

Sres 1 (f, g) = r i=1 (-1) m-mi 1 ≤ j ≤ r j = i g(α j ) mj (α i -α j ) mj g(α i ) mi-1 (x -α i )• k 1 + • • • + k i + • • • • • • + k r+s = m i -1 1 ≤ j ≤ r j = i mj-1+kj kj (α i -α j ) kj 1≤ℓ≤s n ℓ -1+k r+ℓ k r+ℓ (α i -β ℓ ) k r+ℓ + min{1, m i -1} k 1 + • • • + k i + • • • • • • + k r+s = m i -2 1 ≤ j ≤ r j = i mj -1+kj kj (α i -α j ) kj 1≤ℓ≤s n ℓ -1+k r+ℓ k r+ℓ (α i -β ℓ ) k r+ℓ .
This is a nontrivial expression, and nothing similar has been found yet for subresultants of general orders. It is worth noticing, however, that determinantal formulations for subresultants of square-free polynomials readily generalize to the case of polynomials with multiple roots (see [DKS2013, Th.2.5]), so that the difficulty seems to lie in finding expanded expressions.

In this article we take a completely different approach and focus on an extremal case, which is when both f and g have only one multiple root each: we get explicit expressions for Sres d ((x-α) m , (x-β) n ) for all d < min{m, n}.

To do this, we set 0 < d < min{m, n} or d = min{m, n} when m = n and c = c(m, n, d) := m + n -2d -1. We introduce the d × (d + 1) integer Hankel matrix with binomial entries by

H(m, n, d) := c m -i -j 1≤i≤d 0≤j≤d =      c m-1 c m-2 . . . . . . c m-d-1 c m-2 . . . . . . c m-d-2 . . . . . . . . . . . . c m-d c m-d-1 . . . . . . c m-2d     
, where, by convention, c k = 0 for k < 0 and for k > c. Denote with q j (m, n, d) the j-th maximal minor of H(m, n, d) defined as the determinant of the square submatrix H j (m, n, d) of H(m, n, d) obtained by deleting its (j + 1)-th column, for 0 ≤ j ≤ d. By convention, q 0 (m, n, 0), the determinant of an empty matrix, equals 1.

Clearly, all q j (m, n, d) are integer numbers. To regard them as elements of the field K, we consider their class via the natural ring homomorphism Z → K which maps the integer 1 to the unit 1 K of K.

We now describe our main result, which provides a closed-form expression for the coefficients of the subresultant Sres d ((x-α) m , (x-β) n ) when expressed in the set of Bernstein polynomials (x -α) j (x -β) d-j , 0 ≤ j ≤ d .

Theorem 1.1. Let m, n, d ∈ N with 0 ≤ d < min{m, n}, and α, β ∈ K. Then,

Sres d ((x-α) m , (x-β) n ) = (-1) ( d 2 ) (α-β) (m-d)(n-d) d j=0 q j (m, n, d)(x-α) j (x-β) d-j .
Note that Theorem 1.1 is consistent with the Poisson formula (2) for d = 0.

Our second result completes the first one by providing explicit expressions for the values of the minors q j (m, n, d), 0 ≤ j ≤ d, as products of quotients of explicit factorials.

Theorem 1.2. Let m, n, d ∈ N with 0 < d < min{m, n}, and c = m + n -2d -1. Then,

q 0 (m, n, d) = (-1) ( d 2 ) d i=1 (i -1)! (c + i -1)! (m -i -1)!(n -i)! ,
and for 1 ≤ j ≤ d the following identities hold in Q:

q j (m, n, d) = d j n-d+j-1 j m-1 j q 0 (m, n, d).
The proof of Theorem 1.1 yields as a byproduct (see Proposition 3.3) a nice description of the d-th principal subresultant PSres

d ((x -α) m , (x -β) n ), that is, of the coefficient of x d in Sres d ((x -α) m , (x -β) n ): PSres d ((x -α) m , (x -β) n ) = (α -β) (m-d)(n-d) d i=1 (i -1)! (c + i)! (m -i)!(n -i)! . (3)
The product in (3) is an integer number whose prime factors are less than m+n-d. Thus, if α = β and if the characteristic of K is either zero or at least equal to m+n-d, the subresultant Sres d ((x-α) m , (x-β) n ) is a polynomial of degree exactly d. When char(K) is positive but smaller than m+n-d, this is generally not true (though exceptions exist, e.g., for m = 5, n = 3, d = 2 and char(K) = 3). The change of behavior might be very radical. For instance, there exist triples (m, n, d) for which the degree of Sres

d ((x-α) m , (x-β) n ) is less than d for any positive characteristic p < m + n -d. Such an example is (m, n, d) = (6, 8, 2). Another interesting example is when p = m + n -d -1: in that case, Sres d ((x-α) m , (x-β) n ) reduces to a constant in characteristic p.
In general, the degree of Sres d ((x -α) m , (x -β) n ) can be determined using Theorem 1.2. For example, in characteristic 5, the order-8 subresultant of (x -α) 11 and (x -β) 9 is a polynomial of degree 6 for all α = β.

We briefly sketch our proof strategy for these results. We start from the basic fact that if Sres d ((x-α) m , (x-β) n ) has degree exactly d, then any linear combination F •(x-α) m +G •(x-β) n , of degree bounded by d with deg(F ) < n-d and deg(G) < m-d, is a scalar multiple of the subresultant (Lemma 3.1). In Proposition 3.4 we show that d j=0 q j (m, n, d)(x -α) j (x -β) d-j can be expressed as such a linear combination, and determine the scalar multiple which is the ratio between this expression and the subresultant. Theorem 1.1 then follows by specializing the "generic case" (in characteristic zero) to fields of positive characteristic.

To prove Theorem 1.2, we proceed in two main steps. We first evaluate q 0 (m, n, d) in Lemma 2.3 by using a result due to Ostrowski (Lemma 2.1) for the determinant evaluation of a Hankel matrix involving binomial coefficients. Then, in Lemma 2.4 we reduce the computation of the remaining q j (m, n, d), 1 ≤ j ≤ d, to a hypergeometric identity due to Pfaff and Saalschütz (Lemma 2.2).

Sylvester's original motivation for deriving expressions in roots for subresultants was to understand how Sturm's method for computing the number of real roots of a polynomial in a given interval works formally; see [Syl1853], where Sylvester applies the theory of subresultants developed there to the case when g = f ′ , with f having simple roots. Furthermore, Sylvester's formulae opened the door to have great flexibility in the evaluation of resultants and subresultants (see the book of Jouanolou and Apéry [AJ2006] for several ingenious formulae for the simple roots case). The search of explicit expressions for subresultants of polynomials having multiple roots is an active area of research; see for instance [Hon1999, LP2003, DHKS2007, DHKS2009, RS2011, DKS2013]. The quest for such formulae has uncovered some interesting connections of subresultants to other well-known objects and thus several new applications were discovered.

One of these applications are closed expressions for various rational interpolation problems, including the Cauchy interpolation or the osculatory rational interpolation problem [BL2000, DKS2015]. The search for formulae in multiple roots also uncovered the close connection of subresultants to multivariate symmetric Lagrange interpolation [KSV2017]; generalizations to symmetric Hermite interpolation are the topic of ongoing research. These formulae in roots may be used to analyze the vanishing of the coefficients of the subresultants (see the discussion after Theorem 1.2), a question related to the understanding of the performance of the Euclidean Algorithm for polynomials over finite fields [MG1990]. Our closed formulae also led us think about accelerating the computation of the subresultants in our particular extremal case; this will be explained in detail in a forthcoming paper.

The paper is organized as follows: We first derive Theorem 1.2 in Section 2 thanks to Lemmas 2.3 and 2.4. Section 3 then introduces the aforementioned multiple of the subresultant and proves Theorem 1.1. search project MTM2013-40775-P, Teresa Krick and Marcelo Valdettaro are partially supported by ANPCyT PICT-2013-0294 and UBACyT 2014-2017-20020130100143BA, and Agnes Szanto was partially supported by NSF grant CCF-1217557.

Proof of Theorem 1.2

All along this section, we work over the rational numbers to compute the coefficients q j (m, n, d) which appear in the expression of the subresultant given in Theorem 1.1 over a field K of characteristic zero. As these numbers are integers, we can regard them as elements of any field K via the natural ring homomorphism Z → K which maps 1 Z → 1 K .

We start by recalling Ostrowski's determinant evaluation (Lemma 2.1) for Hankel matrices with binomial coefficients entries, and the Pfaff-Saalschütz identity (Lemma 2.2) for the evaluation at the point 1 of a special family of 3 F 2 hypergeometric functions.

Lemma 2.1 ([Ost1964]

). For ℓ, k ∈ N and a 0 , a 1 , . . . , a k ∈ N,

det ℓ a i -j 0≤i,j≤k = ℓ! k+1 k i=1 (ℓ + i) k+1-i 0≤i<i ′ ≤k (a i ′ -a i ) k i=0 a i ! k i=0 (ℓ + k -a i )! . Lemma 2.2 ([Pfa1797, Saa1890, And1996, And1997], [Sla1966, §2.3.1]).
Let x, y, z be indeterminates over Q. Then, for any k ∈ N, the following identity holds in Q(x, y, z):

k j=0 (x) j (y) j (-k) j (z) j (1 + x + y -z -k) j j! = (z -x) k (z -y) k (z) k (z -x -y) k .
Here (x) 0 := 1 and (x) j := x(x + 1) • • • (x + j -1) for j ≥ 1 denotes the j-th Pochhammer symbol of x.

By applying these two results, Theorem 1.2 follows straightforwardly from Lemmas 2.3 and 2.4 below. Lemma 2.3 computes q 0 (m, n, d) as a direct consequence of Ostrowski's determinant evaluation. Lemma 2.4 computes all q j (m, n, d) for j > 0, and is a consequence of a binomial identity (given in (6)) which is, in fact, the Pfaff-Saalschütz identity in disguise. Recall that we have set c = m + n -2d -1.

Lemma 2.3. Let d, m, n ∈ N with 0 < d < min{m, n}. Then, q 0 (m, n, d) = (-1) ( d 2 ) d i=1 (i -1)!(c + i -1)! (m -i -1)!(n -i)! .
Proof. Using Lemma 2.1 with k = d -1 and a i = m -i -2 for 0 ≤ i ≤ d -1, and ℓ = c we get

q 0 (m, n, d) = det c m -i -j 1≤i,j≤d = det c m -i -j -2 0≤i,j≤d-1 = c! d d-1 i=1 (c + i) d-i 0≤i<i ′ ≤d-1 (i -i ′ ) d-1 i=0 (m -i -2)! d-1 i=0 (c + d -1 -(m -i -2))! = d i=1 c! i-1 j=1 (c + j) • (-1) ( d 2 ) d i=1 (i -1)! d i=1 (m -i -1)! d i=1 (n -d + i -1)! .
The statement follows by rearranging terms. ✷ Lemma 2.4. Let j, d, m, n ∈ N with 0 < j ≤ d < min{m, n}. Then,

q j (m, n, d) = d j n-d+j-1 j m-1 j q 0 (m, n, d).
Proof. Observe that the matrix H has full rank d since by Lemma 2.3, its minor q 0 (m, n, d) is non-zero. Therefore, an elementary linear algebra argument shows that the kernel of the induced linear map H: Q d+1 → Q d has dimension 1, and is generated by the (non-zero) vector q(m, n, d) := (q 0 (m, n, d), -q 1 (m, n, d), . . . , (-1) d q d (m, n, d)).

Set k j (m, n, d) := Therefore, to prove (4) it is enough to check the following identities

d j=0 m + n -2d -1 m -j -i (-1) j k j (m, n, d) = 0 for 1 ≤ i ≤ d. ( 5
)
We actually prove that a more general identity holds for any i ∈ N:

d j=0 m + n -2d -1 m -j -i (-1) j d j n-d+j-1 j m-1 j = i-1 d m+n-d-1 m-i m-1 d , (6)
The expressions in (5) are then recovered by specializing i to 1, . . . , d.

The equalities in (6) follow from the Pfaff-Saalschütz identity described in Lemma 2.2. Since both sides of (6) are polynomials in n (of degree at most m -i) it is enough to verify them for an infinite number of values n. We will show that they hold for n ≥ 2d.

By observing that (a + j -1)! = (a -1)!(a) j , a+j-1 j = (a) j j! , (a -j)! = (-1) j a! (-a) j and a j = (-1) j (-a) j j! , we deduce that the left-hand side of (6) is equal, for n ≥ 2d, to

(m + n -2d -1)! (m -i)!(n -2d + i -1)! • d j=0 (n -d) j (-(m -i)) j (-d) j (n -2d + i) j (-(m -1)) j j! .
To simplify the latter sum, we now apply Lemma 2.2 for k = d, and for x, y, z specialized respectively to n -d, -(m -i), n -2d + i, and get

d j=0 (n -d) j (-(m -i)) j (-d) j (n -2d + i) j (-(m -1)) j j! = (i -d) d (m + n -2d) d (n -2d + i) d (m -d) d = i-1 d (m + n -2d) • • • (m + n -d -1) m-1 d (n -2d + i) • • • (n -d + i -1) ,
from which (6) follows immediately. ✷ Therefore, by the definition of the principal subresultant,

PSres d ((x -α) m , (x -β) n ) = = det m+n-2d 1 • • • 0 • • • 0 . . . . . . . . . n-d . . . . . . . . . 1 0 • • • 0 1 • • • n d (α -β) n-d • • • n m-1 (α -β) n-(m-1) . . . . . . . . . m-d 1 . . . n 2d-m+1 (α -β) n-(2d-m+1) . . . n d (α -β) n-d = det m-d n d (α -β) n-d . . . n m-1 (α -β) n-(m-1) n d-1 (α -β) n-(d-1) . . . n m-2 (α -β) n-(m-2) m-d . . . . . . n 2d-m+1 (α -β) n-(2d-m+1) . . . n d (α -β) n-d = (α -β) (m-d)(n-d) det n d -i + j 1≤i,j≤m-d = (α -β) (m-d)(n-d) d i=1 (i -1)! (c + i)! (m -i)!(n -i)! .
The third equality above follows from the "weighted" homogeneities of the determinant. Indeed, by multiplying the i-th row in the second matrix above by (α -

β) i-1 , 1 ≤ i ≤ m -d, the whole determinant gets multiplied by (α - β) 1+•••+(m-d-1) = (α -β) ( m-d 2 )
, but now for each j = 1, . . . , m -d, column j has the same term (α -β) n-d+j-1 that can be factored out, obtaining (α -

β) (n-d+m-d-1)+•••+(n-d+0) = (α-β) (m-d)(n-d)+( m-d
2 ) and one can then clear out the spurious (α -β) ( m-d

2 ) ), and the equality can be derived from Lemma 2.1 with ℓ = n, k = m -d -1 and a j = d + j + 1 for 0 ≤ i ≤ m -d -1.

✷

We now show how to express a scalar multiple of the polynomial expression

d j=0 q j (m, n, d)(x -α) j (x -β) d-j as a polynomial combination F • (x -α) m + G • (x -β) n , with
h d (α, β, m, n) = F • (x -α) m + G • (x -β) n .
Proof. Set f := (x -α) m and g := (x -β) n , and write

(α -β) c = (α -x + x -β) c = c k=0 (-1) k c k (x -α) k (x -β) c-k . Fix 0 ≤ j ≤ d. Then, (α -β) c (x -α) j (x -β) d-j = c k=0 (-1) k c k (x -α) k+j (x -β) c-k+d-j .
For k + j ≥ m the corresponding terms in the right-hand side are polynomial multiples of f , with coefficient F j of degree bounded by (k + j)

+ (c -k + d -j) -m = n -d -1.
Similarly, for c -k + d -j ≥ n, the corresponding terms are multiples of g, with coefficient G j of degree bounded by (k + j)

+ (c -k + d -j) -n = m -d -1. The remaining terms satisfy k + j < m, i.e. k < m -j and c -k + d -j < n, i.e. k > m -j -d -1. Therefore (α -β) c (x -α) j (x -β) d-j = F j f + G j g + m-j-1 k=m-j-d (-1) k c k (x -α) k+j (x -β) c-k+d-j = F j f + G j g + d i=1 (-1) m-i-j c m -i -j (x -α) m-i (x -β) n-d+i-1 .
Multiplying each of these equations by q j (m, n, d) for 0 ≤ j ≤ d and adding them up, we get

h d (α, β, m, n) = (α -β) c d j=0 q j (m, n, d)(x -α) j (x -β) d-j = F f + G g + d j=0 d i=1 (-1) m-i-j c m -i -j q j (m, n, d)(x -α) m-i (x -β) n-d+i-1 , with F := d j=0 q j (m, n, d)F j and G := d j=0 q j (m, n, d)G j . It turns out that d j=0 d i=1 (-1) m-i-j c m -i -j q j (m, n, d) (x -α) m-i (x -β) n-d+i-1 = d i=1 (-1) m-i (x -α) m-i (x -β) n-d+i-1 d j=0 (-1) j c m -i -j q j (m, n, d) = 0,
since, as observed in the proof of Lemma 2.4, (q 0 (m, n, d), -q 1 (m, n, d), . . . , (-1) d q d (m, n, d)) generates ker H.

Therefore h d (α, β, m, n) = F • (x -α) m + G • (x -β) n with deg(F ) < n -d and deg(G) < m -d. ✷
We now compute explicitly the d-th coefficient of h d (α, β, m, n), which also implies in particular that it has degree exactly d when α = β and char(K) = 0 or char(K) ≥ m + n -d. Proposition 3.5. Let d, m, n ∈ N with 0 < d < min{m, n}, and α, β ∈ K. Then,

coeff x d h d (α, β, m, n) = (-1) ( d 2 ) (α -β) c d i=1 (i -1)! (c + i)! (m -i)!(n -i)! . For d = 0 we have h 0 (α, β, m, n) = (α -β) c . Proof. It is clear that coeff x d (h d (α, β, m, n)) = (α-β) c d j=0 q j (m, n, d).
The d = 0 case follows from our convention that q 0 (m, n, 0) = 1. We now show that d j=1 q j (m, n, d) = q 0 (m + 1, n, d), which proves the statement by Lemma 2.3. Observe that To prove these equalities, we use the identities shown in Propositions 3.3 and 3.5. The final identity for µ also holds when d = 0. Plugging the expression of h d given in (7) in the identity (8), we deduce Theorem 1.1 in this case.

In the general case, we use the fact that Theorem 1.1 holds for (x -u α ) m and (x -u β ) n in K ⊃ Q(u α , u β ), where u α , u β are indeterminates over Q.

As subresultants are defined via the determinant (1), and in this case they actually belong to Z[u α , u β ][x], the expression (1.1) holds after specializing u α → α, u β → β, and the standard ring homomorphism Z → K. This concludes the proof of Theorem 1.1. ✷

  j ≤ d. It suffices then to show that k(m, n, d) := (k 0 (m, n, d), -k 1 (m, n, d), . . . , (-1) d k d (m, n, d)) ∈ ker H, (4) so then we would have k(m, n, d) = λ q(m, n, d) with λ = 1/q 0 (m, n, d), as k 0 (m, n, d) = 1.

F

  and G satisfying the hypothesis of Lemma 3.1. For this, we define for 0 ≤ d < min{m, n},h d (α, β, m, n) := (α -β) c d j=0 q j (m, n, d)(x -α) j (x -β) d-j . (7) Note that h d (α, β, m, n) ∈ K[x]has degree bounded by d. Proposition 3.4. Let d, m, n ∈ N with 0 ≤ d < min{m, n}, and α, β ∈ K. There exist F , G ∈ K[x] with deg(F ) < n -d, deg(G) < m -d such that

.

  For 0 ≤ j ≤ d let C(j) denote the (j + 1)-th column of the matrix above. We perform the following operations: C(j) + C(j -1) → C(j) for j = d, . . . , 0. By using the identity c k-1 + c k = c+1 k ,

=

  q 0 (m + 1, n, d).✷We are now ready to prove Theorem 1.1.Proof of Theorem 1.1. When α = β, both sides of the expression in Theorem 1.1 vanish. Assume now that α = β, and that char(K) = 0, or char(K) ≥ m + n -d. Thanks to Propositions 3.3 and 3.5, both Sres d ((x -α) m , (x -β) n ) and h d (α, β, m, n) are non-zero polynomials of degree exactly d. Recall that we have set c = m + n -2d -1. Proposition 3.4 and Lemma 3.1 with µ = 1/λ then imply thatSres d ((x -α) m , (x -β) n ) = µ • h d (α, β, m, n) d ((x -α) m , (x -β) n ) coeff x d (h d (α, β, m, n)) = (α -β) (m-d)(n-d) d i=1 (i-1)! (c+i)! (m-i)!(n-i)! (-1) ( d 2 ) (α -β) c d i=1 (i-1)! (c+i)! (m-i)!(n-i)! = (-1) ( d 2 ) (α -β) (m-d)(n-d)-c .

The Sylvester matrix was defined in[Syl1840], and the order-d subresultant was introduced in[Syl1839, Syl1840] under the name of "prime derivative of the d-degree".

Acknowledgements

This project started when the last four authors met at the FoCM Conference in Montevideo in December 2014 and at the University of Buenos Aires in September 2015. We are grateful to Professor Richard A. Brualdi, who suggested a joint collaboration with the first author. Carlos D'Andrea was partially supported by ANPCyT PICT-2013-0294, and the MINECO re-

Proof of Theorem 1.1

To prove Theorem 1.1, we make use of the following well-known result, which follows for instance from Lemmas 7.7.4 and 7.7.6 in [Mis1993].

Following Lemma 3.1, we first prove that Sres d ((x -α) m , (x -β) n ) has indeed degree d when α = β and char(K) = 0 or char(K) ≥ m + n -d, in other words, that its principal subresultant PSres d ((x -α) m , (x -β) n ) is non-zero. We start by recalling a well-known result, which is used in the proof.

Proposition 3.3. Let d, m, n ∈ N with 0 < d < min{m, n}, and α, β ∈ K. Then,

In particular, if α = β and char