Modeling of Chemical Reaction Systems with Detailed Balance Using Gradient Structures

We consider various modeling levels for spatially homogeneous chemical reaction systems, namely the chemical master equation, the chemical Langevin dynamics, and the reaction-rate equation. Throughout we restrict our study to the case where the microscopic system satisfies the detailed-balance condition. The latter allows us to enrich the systems with a gradient structure, i.e. the evolution is given by a gradient-flow equation. We present the arising links between the associated gradient structures that are driven by the relative entropy of the detailed-balance steady state. The limit of large volumes is studied in the sense of evolutionary \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma $$\end{document}Γ-convergence of gradient flows. Moreover, we use the gradient structures to derive hybrid models for coupling different modeling levels.


Introduction
In this work we discuss different models for chemical reactions taking place in a container of volume V . Throughout we assume that the spatial extent of the container and the position of the chemical species are irrelevant, which means that we are looking at a well-stirred system. We assume that the system is composed of I different species named X 1 to X I , which may represent different molecules, e.g., X 1 = H 2 , X 2 = O 2 , and X 3 = H 2 O. We assume that these I species undergo R different reactions of mass-action type: where the vectors α r , β r ∈ N I 0 contain the stoichiometric coefficients, and k r fw , k r bw > 0 are the forward and backward reaction rates, see Sect. 2. The reaction 2H 2 + O 2 − − 2H 2 O would lead to the vectors α = (2, 1, 0) and β = (0, 0, 2).
Denoting by c = (c 1 , . . . , c I ) ∈ C := [0, ∞[ I the vector of nonnegative densities, the simplest model is the macroscopic reaction-rate equation (RRE), which is a system of ODEs on the state space C: Here the monomials c α r := I i=1 c α r i I indicate that the probability for the right number of particles for the r th reaction to meet is given by a simple product of the corresponding densities, i.e., we assume that the positions of the particles are independent.
A truly microscopic model can be obtained as a stochastic process. Here we count the number of particles N V i (t) for each species X i and consider the random vector N V (t) = (N V 1 (t), . . . , N V I (t)) ∈ N := N I 0 . A forward or backward reaction of type r is modeled as an instantaneous event where the particle numbers jump from N V (t)+α r to N V (t)+β r or vice versa. The corresponding jump rates in a volume of size V > 0 are given by k r fw B α r V (N V (t)) and k r bw B β r V (N V (t)) respectively; see (3.1) for the definition of B α V (n) . Here we study the vector of probabilities that describes the probability distribution of the random variable N V (t). The time evolution of u(t) is given by the chemical master equation (CME), i.e., the Kolmogorov forward equation associated with the continuous time Markov chain above. This is a countable linear system of ODEs:u where B V is an (unbounded) linear operator on 1 (N ), see Sect. 3, where also existence and uniqueness of solutions is discussed. We refer to [39] for a short introduction to the CME and to [23] for a justification. The basis of this work is the observation from [40] that (RRE) can be interpreted as a gradient flow if the reaction system satisfies the detailed-balance condition, i.e., there exists a positive equilibrium c * = (c * i ) i=1,...,I ∈ ]0, ∞[ I such that κ r * := k r fw c α r * = k r bw c β r * for r = 1, . . . , R.
Defining the Boltzmann entropy E and the Onsager operator K via with λ B (z) = z log z − z + 1 and If (RRE) satisfies the detailed-balance condition, then (CME) does so with an equilibrium distribution w V ∈ P(N ) that is explicitly given as a product of one-dimensional Poisson distributions with mean c * i V , namely (cf. Theorem 3.1), Consequently, we are also able to interpret (CME) as a gradient flow induced by a gradient system (P(N ), E V , K V ), see (3.7). Here E V (u) is again the Boltzmann entropy with respect to w V , but now divided by the volume V : (1.2)

Large-volume approximations using gradient structures
A major challenge in modeling chemical reactions is the question of understanding the transition from small-volume effects to the macroscopic behavior in large volumes. The first breakthrough was obtained in [30][31][32][33] by connecting the particle numbers N(t) ∈ N to the concentrations c ∈ C and showing that almost surely for all t > 0, (1.3) where t → c(t) is the solution of (RRE) with c(0) = c 0 . This result may be interpreted as a justification for the RRE in terms of the Markovian model. In [45,46] a dynamic large deviation principle is applied to 1 V N V (·), which leads to a rate functional that generates a gradient structure (C, E, cosh ); see Sect. 2.5. Recent large deviation results for chemical reaction networks can be found in [1,2].
In this paper we study the limit V → ∞ for the gradient system (P(N ), E V , K V ), and hence for (CME), in the sense of evolutionary -convergence for gradient systems, as introduced in [54,58] and further developed in [13,41]. For this purpose we use a suitable embedding ι V : P(N ) → P(C) (Sect. 4) and obtain the coarse grained gradient system (P(C), E, K ) with

E(c) (dc) and K ( )ξ (c) = −div c (c)K(c)∇ c ξ(c) .
In particular, the coarse grained gradient flow equation is the Liouville equatioṅ (t, c) = div c (t, c)R(c) , associated with (RRE); here we used that ξ = D E = E and R = −KD c E. Thus, in this scaling a pure transport equation remains, while all diffusion disappears, as can be seen in the factor 1/V before the middle sum in (1.2). In particular, our result is consistent with Kurtz' result (1.3): by assuming (0) = δ c 0 ∈ P(C) we obtain (t) = δ c(t) . While Kurtz works directly on the Markovian random variables, we work at the level of their distributions: u ∈ P(N ) CME:u = −K V (u)DE V (u) Our convergence result for the gradient systems (P(N ), E V , K V ) to the limiting gradient system (P(C), E, K ) can be seen as a concrete example of the EDP convergence of gradient systems as discussed in [13,35]. Another example treating the convergence of "Markovian discretizations" towards a Fokker-Planck equation is studied in [11]; see also [16,18,56] for applications to interacting particle systems.
In addition to the extreme cases V finite and V → ∞ it is also important to study the case of intermediate V , where 1 V N V (t) already behaves continuously but still shows some fluctuations of standard deviation 1/ √ V , see [61] for a numerical approach to treat the hierarchy via a suitable hybrid method. In [34] it is shown that the random vector t → X V (t) ∈ C obtained by solving the stochastic differential equation with independent Brownian vectors B fw (t), B bw (t) ∈ R R , and fw (X) = κ r X α r (see [34,Eq. (1.7)]) yields an improved approximation because 1 This model is a so-called diffusion approximation, which in the reaction context also is termed 'chemical Langevin dynamics'. In [24,Eq. (23)] and [61, Eq. (7)] the stochastic differential Eq. (1.4) is called chemical Langevin equation (CLE).
The associated Kolmogorov forward equation takes the forṁ (1.5) Here the diffusion matrix K CLE can be written in the explicit form that is different from K(c), because in the former the arithmetic mean while in the latter the logarithmic mean is taken. One drawback of the chemical Langevin Eq. (1.5) is that it cannot be written as gradient flow of the relative entropy, as the Einstein relation for the drift flux and the diffusion flux is not satisfied. Therefore we propose other approximations that stay inside the theory of gradient flows and seem to work sufficiently well if the concentrations are not too large or small. Our simplest approximation is given by the gradient system (P(C), E V , K ) with which leads to the linear Fokker-Planck equatioṅ In Sect. 5 we show that by systematically deriving higher-order corrections to E V and K we can recover the asymptotically correct diffusion matrix K CLE while keeping the gradient structure, but have to accept several additional terms, or switch over to the notion of asymptotic gradient flow structures in the sense of [6].

Hybrid modeling using gradient structures
A major advantage of the gradient flow description is that the different structures can be combined to obtain hybrid models, in which the set of chemical species is divided into subclasses which may be treated differently depending on the desired or needed accuracy. Our approach is based on the idea of model reduction for gradient structures. The idea is to approximate a complicated gradient structure (X, via an embedding mapping x = (y). Staying within the class of gradient systems has the advantage that the most important features of the original system can be preserved. In particular, decay of the driving functional along the approximate flow holds automatically. By contrast, such crucial features could get lost in a direct approach based on the evolution equation itself. In Sect. 6 we shall deal with three examples for hybrid models where it is essential to keep V as a large but finite parameter. First, we shall consider a hybrid model in which an RRE is coupled to a Fokker-Planck equation. Here the set of species is divided into two classes: C = C s ×C m . Some of them will be described stochastically (s), while others are described macroscopically (m). This leads to a gradient flow structure on the hybrid state space Y = P(C s )×C m . The resulting gradient flow equation turns out to be a meanfield equation, in which the density of the component c s satisfies a linear equation which is nonlinearly coupled to an ODE for the component c m .
We also study the coupling of an RRE for macroscopic variables to a CME for n microscopic variables. This leads to a hybrid system on P(N n 0 )×C m . Finally we analyze a mixed CME / Fokker-Planck model with state space P(N), in which the underlying space N := {0, 1, . . . , N −1} ∪ [N /V , ∞[ contains a mixture of discrete and continuous components.
The present work concentrates solely on the analytical underpinnings of hybrid modeling for CME; for numerical approaches to CME and to spatio-temporal CME we refer to [4,12,15,[27][28][29]48,61]. Geodesic convexity properties for the gradient structures appearing in this paper are studied in [38]. Notational conventions Throughout the paper we will consistently use the following notation to distinguish the different modeling levels.

Liouville equation
The LE is denoted by (P(C), E, K ): state and state space = ρ dc ∈ P(C), steady state δ c * , dual variable ξ energy functional E( ) = C E(c) d (c), Onsager operator K ( ) = −div (·)K(·)∇ Fokker-Planck equation The FPE is denoted by (P(C), E V , K ): state and state space ∈ P(C) := P(C), steady state W V , dual variable ξ Hybrid systems are denoted by "mathfrak" letters: for coupling CME and RRE (P(N V ,N ), E V ,N , K V ,N ) for merging discrete and continuous modeling for one species.
The space of all signed Borel measures of bounded variation on C is denoted by M (C).

Reaction Rate Equations
We denote by c = (c 1 , . . . , c I ) ∈ C := [0, ∞[ I the concentrations of I different chemical species X 1 , . . . , X I reacting according to the mass action law, i.e., the reactions where R is the number of possible reactions, α r , β r ∈ N I 0 are the vectors of the stoichiometric coefficients, and k r fw , k r bw > 0 are the forward and backward reactionrates. In general these rates may depend on c, but for simplicity we keep them as constants in this work. A typical example is the splitting of water into hydrogen and oxygen, namely The corresponding reaction-rate equations (RRE) are given via the ODE systeṁ where c α = c α 1 1 · · · c α I I , see [17,19,26].

Stoichiometry, Conservation, and Decomposition of the State Space
The stoichiometric subspace S ⊂ R I and its orthogonal complement S ⊥ are defined via For each ξ ∈ S ⊥ the function C ξ (c) = ξ · c defines a first integral, which easily follows from ξ · R(c) ≡ 0. These conservation laws often go under the name conservation of atomic species, see [17]. Suppose now that S ⊥ is a non-trivial subspace of R I . We shall argue that the RRE induces a decomposition of the state space C = [0, ∞[ I into affine invariant subsets.
(If S ⊥ = {0}, the only invariant set is C itself.) Choosing a basis { m k ∈ R I | k = 1, . . . , m W } of S ⊥ we define the matrix Q ∈ R m W ×I , which has the rows m k ∈ R I . By construction we have Q[S] = {0}, and we conclude that the solutions c of (2.2) conserve Qc as follows: By construction every affine conserved quantity is of the form ξ · c + q for some ξ ∈ S ⊥ and q ∈ R. This allows us to decompose the full state space C = [0, ∞[ I into the invariant, affine subsets (c 0 +S) ∩ C for c 0 ∈ C. Using the notation we define, for all q ∈ Q, the sets Then, q 1 = q 2 implies I(q 1 ) ∩ I(q 2 ) = ∅, and we have C = q∈Q I(q). Let us note that this decomposition does not depend on the choice of the orthonormal basis which determines the matrix Q, although the set I(q) does depend on Q. Note also that we can always write I(q) = (c+S) ∩ C for some arbitrary c ∈ C satisfying Qc = q.

Detailed Balance and the Wegscheider Matrix
We say that the above reaction system fulfills the condition of detailed balance if there exists a positive equilibrium density vector c * ∈]0, ∞[ I such that all reactions are simultaneously in equilibrium, i.e., This condition implies that R(c * ) = 0, but we emphasize that this condition is stronger in general cases. The condition of detailed balance is also called the condition of microscopic reversibility, see [17, p. 45] or [10] for a general discussion of these concepts. We are looking for a characterization of detailed balance. Let W ∈ Z R×I be the matrix which has the row vectors γ r := α r −β r ∈ Z I , r = 1, . . . , R. We call W the Wegscheider matrix because of the pioneering work in [60]. We then have S = Ran W T and S ⊥ = Ker W, which explains the abbreviation m W := dim S ⊥ = dim Ker W. Since c * is strictly positive, we can take the logarithm of the polynomial conditions (2.6) and find the equivalent linear system These conditions on the reaction coefficients k r fw and k r bw are called Wegscheider conditions (see, e.g., [25,57,59,60]). By choosing a basis of Ker W T and exponentiation they can be rewritten as polynomial conditions without referring to the equilibrium state c * .
Since dim(Ran W) = dim(Ran W T ) = dim S by standard linear algebra, the number of Wegscheider conditions can be expressed as Hence, if the number R of reactions is smaller than the number I of species, the Wegscheider conditions can usually be satisfied easily. Local existence for solutions starting in the interior of C is trivial, as R is a polynomial vector field. Since the relative entropy E is a coercive Liapunov functional, the solutions cannot blow up and stay inside a region B R (0) ∩ C for some R > 0. Moreover, solutions cannot leave this region via the boundary ∂ C, since the vector field is either tangential to ∂ C or points inwards. Indeed, if c j (t 0 ) = 0 for some j, theṅ because each term in the sum is nonpositive: If α r j = β r j or min{α r j , β r j } > 0, then the term is 0. Thus, we are left with the cases (α r j , β r j ) ∈ {(n, 0), (0, n)} for some positive n. In the first case c j (t 0 ) = 0 implies c α r (t 0 ) = 0 and the result follows, and the second case is similar.

The Reaction-Rate Equations as a Gradient System
We show that a RRE satisfying the detailed-balance condition can be generated by a gradient system (C, E, K). Here, the state space C := [0, ∞[ I contains all possible concentration vectors c. The driving functional is the relative entropy E : C → [0, ∞[ and the Onsager matrix K is chosen suitably (recall that λ B (z) = z log z − z + 1 ≥ 0): where the logarithmic-mean function is given via 10) The following result shows that a RRE (2.2) satisfying the detailed-balance condition (2.6) is indeed generated by the gradient system (C, E, K). This was first established in [62, Sect. VII] to derive entropy bounds for hyperbolic conservation laws in reactive flows and was rederived in [40]   Proof Multiplying DE(c) = (log(c i /c * i )) i=1,...,I by α r −β r ∈ R I we obtain which is the denominator of where we used the detailed-balance condition (2.6) in DB = . Thus, the assertion is established.
Summarizing the above derivations, we have rewritten the RRE in thermodynamic forṁ (2.12) which is also called the Onsager principle [51,52]. The latter states that the rate (flux) of a macroscopic variable is given as the product of a symmetric positive definite matrix K and the thermodynamic driving force −μ, see e.g. [10, Ch. X, § 4]. The symmetry K = K is related to microscopic reversibility, i.e., detailed balance, see also [44,45]. Subsequently, we refer to K as the Onsager operator or matrix.
Here we clearly see the advantage of using the Onsager operator K to write the RRE as a gradient system,as opposed to working with the Riemannian tensor: we do not have to take care of the fact that K is not invertible except if S = R I .

Continuous Time Markov Chains as a Gradient System
The forward equation for a reversible CTMC on a discrete space {1, 2, . . . , I } is a special case of the RRE considered above. In this case all reactions are of the form and the reaction rates k i j fw (resp. k i j bw ) are interpreted as the transition rates from i to j (resp. from j to i). The reaction-rate equation is the linear system of ODEṡ (2.13) and the detailed-balance condition for the equilibrium state c * takes the form (2.14) Using this condition, the RRE can be written coordinate-wise aṡ Here we used the notational convention that κ i j * := κ ji * for j < i. The relative entropy E is as above and the Onsager matrix takes the form (2.15) where e i ∈ R I denotes the i-th unit vector. We then have the gradient structure (C, E, K M ), namelyċ This gradient flow structure has been found in the independent works [37] (which deals with Markov chains exclusively) and [40] (in the setting of reaction-diffusion systems, in which Markov chains are implicitly contained). The related work [7] deals with discretizations of Fokker-Plank equations.
In fact, for the construction of gradient structures for Markov chainsċ = Ac we do not need the summation rule for logarithms. Hence, following [37,42] there are more general gradient structures. Choosing a strictly convex function φ : / (α r −β r )·φ c c * , but this quantity can be negative in general. As a consequence, the corresponding Onsager matrix would not be positive definite. This cannot happen for Markov chains (i.e., when α = e i and β = e j ), by virtue of the convexity of φ.
In the following we will mainly concentrate on the gradient structure (C, E, K M ) with the logarithmic entropy, as it is the only one that connects with the RRE.
In the case of RRE, the monomial terms c α can only be generated by the logarithmic Hence, we stick to the relative entropy E defined in (2.9), i.e., φ(z) = λ B (z). However, we may replace the linear Onsager principlė c = −K(c)DE(c) by the more general nonlinear formċ = ∂ ζ * c, −DE(c) .
The case ψ r (ζ ) = 1 2 ζ 2 leads to the quadratic dissipation potential in (2.9), i.e., the functions L r are given in terms of the logarithmic mean. In [3] the choices ψ r (±ζ ) = e ζ − 1 − ζ is used. Based on a derivation via the large deviation principle (see [44][45][46]) a special role is played by the choice of a "cosh-type" function ψ r : Here C * is normalized such that C * (ζ ) = 1 2 ζ 2 + O(ζ 4 ). Hence, the dual dissipation potential takes the form * It is shown in [43,Prop. 4.1] that this generalized gradient structure is distinguished as the only tilt-invariant gradient structure for CTMCs. The tilt-invariance even extends to nonlinear RRE, see [47].

Modeling Discrete Particle Numbers via CME
The chemical master equation (CME) is a CTMC that is defined on the set N = N I 0 where n = (n 1 , . . . , n I ) ∈ N is the vector of particle numbers, see [39] for an introduction. This means that n i ∈ N 0 denotes the number of particles of species X i in a sufficiently big volume, whose size is denoted by V > 0. The modeling assumes that all particles move randomly in this big volume (well-stirred tank reactor) so that they can meet independently. The dynamics is formulated in terms of the probabilities u n (t) = probability that at time t there are n i particles of species X i for i = 1, . . . , I .
All the R reaction pairs may happen independently of each other according to the number of the available atoms needed for the reactions and the reaction coefficients k r fw ≥ 0 and k r bw ≥ 0, respectively. Moreover, the jump intensities k r fw B α r V (n) from n + α r to n + β r and k r bw B β r V (n) from n + β r to n + α r also depend on the volume V , as n i denotes the absolute particle number, while for the reaction the densities c i = n i /V matter. The specific form of B α V (n) (cf. [32,39]) reads To avoid clumsy notation we defined where the factor V indicates that the number of reactions is proportional to the volume of the container, if the densities are kept constant.
The CME associated with the RRE (2.2) is the Kolmogorov forward equation for the The r th forward reaction from n+α r to n+β r can only happen (i.e., We emphasize that the RRE as well as the CME are uniquely specified if the reaction network (2.1), the reaction rates k r fw and k r bw , and the volume V > 0 are given. Hence, there are obviously close relations between both models, in particular for V 1, see [4,5,23,32,61]. So far, we have not used the detailed-balance condition, i.e., we can even allow for k r bw = 0 in the above considerations. In all cases, the Kolmogorov forward equation is an infinitedimensional linear ODE as in Sect. 2.4. The following result shows that the detailed-balance condition is inherited from the RRE to the CME, and moreover a simple equilibrium w V can be given explicitly as a product distribution of individual Poisson distributions, namely This result can also be retrieved from [5] by combining Theorems 4.1 and 4.5 there, where it is shown that the weaker "complex-balance condition" is sufficient to guarantee that the Poisson distribution w V is an equilibrium for CME.
For completeness we give a short and independent proof of the fundamental result that for RRE with detailed balance the associated CME satisfies detailed balance again.
satisfies the detailed-balance condition for the CME (3.2), namely Proof For each reaction we obtain the relation Analogously we obtain the same result for k bw B Using the detailed-balance coefficients ν n,r V we can rewrite the operator B r V from (3.2) in a symmetrically balanced form as where e (m) is the unit vector, i.e., e (m) n = δ n−m . It is important to realize that in general the steady state for the detailed-balance condition is highly non-unique, because of the discrete versions Indeed, choosing n arbitrary and defining w = (w n ) ∈ P(N ) via w n = 1 Z w V n for n ∈ I(n) and w n = 0 elsewhere, we obtain another equilibrium for the CME (3.2). Defining convex combination we obtain a rich family of steady states.
The following counterexamples show that the above result, which is central to our work, cannot be expected for systems not satisfying the detailed-balance condition. Building the CME according to (3.2) based on the two reaction pairs we obtaiṅ For the case a = 2 and b = 1, where the detailed-balance condition holds with and it is easy to check that w V = (e −V V n /n!) n∈N 0 is a steady state. However, for a = 7 and b = 1 the detailed-balance condition fails with c (1) = 7/2 > c * = 2 > c (2) = 1. The CME readṡ An explicit calculation shows that the Poisson distribution w V based on c * = 2, i.e., w V n = e −2V (2V ) n /n!, is not a steady state. Indeed, inserting w V into the right-hand side of the last equation we find (for n ≥ 1)

Example 3.3 (Microscopic versus macroscopic detailed balance)
We may also consider a RRE that looks macroscopically as being in detailed balance, but is generated by a microscopic model that is not in detailed balance. The two reactions ∅ 2 X and 2X 1 ∅ produce the RREċ = 2(1−c 2 ) that has the equilibrium c * = 1. However the CME readṡ Note that the reversible reaction pair 2X 1 − − 1 ∅ yields the same RRE, and its associated CME satisfies the detailed-balance condition.

Existence and Uniqueness of Solutions of CME
In this part we establish well-posedness for the CME. We do this by combining classical results from the theory of Markov chains with abstract semigroup theory. For fixed n 0 ∈ N we construct a special Green's function p t (n 0 , ·). General Markov chain theory (e.g., [36,Ch. 2]) implies that there exist a unique minimal solution [0, ∞[×N (t, n) → p t (n 0 , n) to the backward equatioṅ associated with the CME with initial condition p 0 (n 0 , n) = δ n 0 (n). This minimal solution is non-negative and satisfies p t (n 0 , n) ≥ 0 and n∈N p t (n 0 , n) ≤ 1, but for general CTMC it can happen that the latter inequality is strict, which means that the corresponding Markov chain explodes in finite time. We will show that explosion does not happen for CME with detailed balance. For the functional analytic existence and uniqueness result we use the sequence spaces as well as the weighted spaces with the corresponding norms and the usual modification for p = ∞. Now, we consider the transition semigroup (P t ) t≥0 defined by which we shall study by induction over the number R of reactions using the Trotter-Kato formula, where the detailed-balance condition guarantees that each subsystem is a contraction semigroup on L 2 (N , w V ).

Theorem 3.4 Assume that the detailed balance condition (2.6) holds. Then, the semigroup
Moreover, the semigroup is selfadjoint on L 2 (N , w V ) and Markovian, i.e., P t 1 = 1 for all t ≥ 0.
A related existence result for the Markov semigroup of the CME was established in [22], which however does not apply to the case of reversible RRE, because of the restrictions on the growth of the transition rates.
Proof All of the above statements follow from the general theory of continuous time Markov chains, except for the Markovianity. To show the latter, we first consider the case of a single reaction, thus R = 1. Each of the irreducible components of the state space N is then onedimensional (see also [38]), and the Markov chain is a birth-death chain on a countable (possibly finite) set.
If there exist two components of α − β with opposite sign, then each of the irreducible components of the state space N is finite. Therefore it is clear that the Markov chain does not explode in finite time. Suppose now that all components of α − β have equal sign, say α i − β i ≥ 0 for all i = 1, . . . , I , and at least one component is strictly positive. Then each of the infinite irreducible components of N is of the form for some n (0) ∈ N , and the restricted Markov process is a birth-death process with birth rate b k and death rate d k given by [53,Thm. 11] gives a characterization of non-explosion for birth-death chains; it asserts that the chain is non-explosive if and only if k≥ j≥0 In our setting we have Since the summands tend to ∞ as k → ∞, we infer that the latter sum is infinite; hence the Markov chain is non-explosive, or equivalently P t 1 = 1 (see [36,Thm. 2.33]). The case of multiple reactions follows by induction on the number of reactions R. Indeed, for R ⊆ {1, . . . , R}, let (P R t ) t≥0 denote the semigroup corresponding to the reactions r ∈ R. Then the Trotter product formula for contraction semigroups on L 2 (N , w V ) (see e.g., [9]) asserts that Note that we can apply this formula, since the detailed-balance conditions hold for all reactions simultaneously, hence all of the semigroups are contractive on the same space L 2 (N , w V ). We also observe that the class of finitely supported functions is a core for each of the generators. The Markovianity of P {1,...,R+1} thus follows from the Trotter formula and the Markovianity of P {1,...,R} and P {R+1} .

Remark 3.5
The mere existence of a probability distribution satisfying the detailed-balance equations is not sufficient to guarantee non-explosion of a continuous time Markov chain. It might happen that the chain jumps infinitely often in a finite time interval, see [50,Sec. 3.5] for an example. The previous result shows that this phenomenon does not occur in CME satisfying the detailed-balance condition.
It remains to transfer the results from L 1 (N , This definition of B is consistent with the explicit formula for B given above. Since Q generates a C 0 -semigroup of contractions on L 1 (N , w V ), it follows that B generates a C 0semigroup (P t ) t≥0 of contractions on 1 (N ). Furthermore, since P t preserves positivity and P t 1 = 1, it follows that P(N ) is invariant under the semigroup generated by B.
As an immediate consequence we obtain global well-posedness for the CME in P(N ).

Gradient Structures for CME
Since the CME is the forward equation associated with a reversible CTMC, we can formulate it as a gradient flow in view of Proposition 2.3. Indeed, for a strictly convex function φ : where is defined after (2.16), ν n,r V is given in Theorem 3.1, and e (m) denotes the m-th unit vector in 1 (N ). The following result is then a special case of Proposition 2.3.
..,I , then the associated CME has the gradient structure In the following we will mainly be concerned with the case that E φ V is the logarithmic entropy, where φ is the Boltzmann function λ B (z) = z log z − z + 1. In that case we obtain where the logarithmic mean (a, b) is defined in (2.10). The above definitions do not only restrict to the entropy function φ = λ B , but also introduce a normalization with respect to the volume V . Hence, E V can be seen as an entropy per unit volume. The corresponding scaling of K V was chosen such that the evolution equationu . For later purposes we also provide the cosh-type gradient structure for CME, whose relevance and usefulness is discussed in [21,43,44,46]. Recall the definition of C * in (2.18) and note the special scaling via the volume V in (3.8) below, which is needed because * cosh,V (u, ·) is not scaling invariant. Proposition 3.8 (Cosh-type gradient structure for CME) If the RRE (2.2) satisfies the detailed-balance condition (2.6) for a positive steady state c * = (c * i ) i=1,...,I , then the associated CME has the gradient structure  0)) and consider probability measures (t, ·) ∈ P(C) that are obtained by transporting 0 with t , namely

Liouville and Fokker-Planck Equations
It is now easy to see that t → t ∈ P(C) satisfies the Liouville equation in the sense of distributions. We will regard (4.1) as an evolution equation in the space P(C). We will not always notationally distinguish between an absolutely continuous probability measure and its density, but if we want to distinguish them we will write (dc) = ρ(c) dc with ρ ∈ L 1 (C). The goal of this section is to give a rigorous connection between the CME for V → ∞ and the Liouville equation in terms of the associated gradient structures.

The Liouville Equation as a Gradient System
We show that the gradient structureċ = −R(c) = −K(c)DE(c) for the RRE, which was discussed in Sect. 2.3, induces a natural gradient structure for the Liouville equation. Consider the "Otto-Wasserstein-type" Onsager operator K ( ) that acts on functions ξ : C → R via where div and ∇ are taken with respect to c ∈ R I . We also consider the affine potential energy functional E : P(C) → [0, +∞] defined by In the next result we identify the formal gradient structure for the Liouville equation.
Proof Let ∈ P(C) and let σ ∈ M (C) be a signed measure of finite total variation such that σ (C) = 0 and + hσ ∈ P(C) for |h| sufficiently small. Then we have The gradient flow equation˙ = −K ( )DE( ) is thus given by the Liouville equation (4.1).

Passing to the Limit from CME to Liouville
In this section we shall demonstrate that the gradient flow structure for the CME converges in a suitable sense to the gradient structure for the Liouville equation if V → ∞. More precisely, we will show that after a suitable V -dependent embedding of P(N ) into P(C) the proper scalings of the functionals E V and * V : (u, μ) → 1 2 μ · K V (u)μ converge in the sense of -convergence to the corresponding structures for the Liouville equation given by the gradient system (P(C), E, K ), see Sect. 4.3 to 4.5. Following the approach in [41,54,58], and in particular [35], we are then able to establish the convergence for V → ∞ of solutions The main tool for proving this evolutionary -convergence for gradient systems is the so-called energy-dissipation principle, cf. [41,Sect. 3.3], which states that u V solves the CME if and only if for all T > 0 the following energy-dissipation estimate holds: where we use the quadratic dissipation potential V and its Legendre dual * V defined via * where the second form uses Theorem 3.1 and is especially useful to perform the limit V → ∞, see the proof of Proposition 4.6.
We refer to [8,41] for this equivalence and general methods for proving such results on evolutionary -convergence. In [11] a similar approach was used to establish the convergence of CTMC to a Fokker-Planck equation. However, there the convergence of a parabolic equation is established, where upper and lower bounds of the density can be used. Here, the importance is that our limit measures (t) may not have densities; indeed, because we want to recover the Kurtz result (1.3) we are interested in the "deterministic case" (t) = δ c(t) . So our analysis has to be more careful in dealing with general limit measures. For this, we use the dualization approach introduced in [35 for suitably chosen recovery functions t → μ V (t). In order to compare probability measures on different spaces N and C, we consider a suitable embedding ι V : P(N ) → P(C). Here ι V (u) is simply obtained by assigning the mass of u at n ∈ N uniformly to the cube More explicitly, ι V (u) is given by where 1 1 A denotes the indicator function with 1 1 A (b) = 1 for b ∈ A and 0 otherwise. The corresponding dual operation acting on functions ξ ∈ C b (C) is given by The final convergence result will be formulated in Theorem 4.7, which will be a direct consequence of the following three estimates where the dual dissipation potential * Lio is defined via * We will see in Sect. 4.6 that the limsup estimate for the dual potential * V in Sect. 4.5 provides a weak form of a liminf estimate for the primal potential V .
A fundamental fact of the chosen gradient structures of the underlying Markov processes is that all the three terms in the energy-dissipation principle define convex functionals, which is of considerable help in proving the desired liminf estimates. Note that the convergence ι V (u V ) * is rather weak. However, we can use that the coefficients of the transition rates defining the CME are quite regular, so that the other parts in the integral converge in a much better sense. Moreover, the functionals → E( ) and → * Lio ( , DE( )) are in fact linear in .

0-Limit of the Relative Entropies
We also define X V := ι V (P(N )) ⊂ P(C) and W V = ι V (w V ) ∈ X V and consider the functionals These definitions are chosen such that Finally we define a natural inverse of ι V , namely such that P V := ι V • κ V is a projection from P(C) onto X V ⊂ P(C).
To understand the limit of E V for V → ∞ we will use the representation In Lemma 4.2 below we will show that E V converges pointwise to E as defined in (2.9).

Lemma 4.2 (Pointwise bound for E V )
For all c * > 0 there exist K * > 0 and V * > 0 such that for all V ≥ V * the following bounds hold: Proof We decompose the error via with n defined by c ∈ A V n . For the second term we use the convexity of λ B and the estimate log z ≤ 1 + λ B (z). Hence, we have Summing this inequality over i = 1, . . . , I we obtain the upper bound For the opposite direction we use (a) that λ B decreases on [0, 1] and the convexity of if V ≥ V * 1 := max max{1/c * i , ec * i } i = 1, ..., I , where V c * i ≥ 1 and V ≥ ec * i are needed in (a) and (c), respectively. Together with (4.13) this controls the second error term in (4.12), viz. where For controlling the first error term in (4.12) we use (4.9) and obtain the identity with k n from (4.10). Because of 2πk n ≥ 1 we obtain, for all V ≥ 1, the lower bound For the upper bound we use 2πk 0 = 1 and 2πk n ≤ 8n for n ≥ 1. Hence for n i ≥ 1 we obtain, using again the estimate log z ≤ 1 + λ B (z), Summation over i = 1, . . . , I yields, for all c ∈ A V n and V ≥ V * 3 := 8e max{c * 1 , . . . , c * I }, the upper bound Together with the lower estimate we control the first error term in (4.12) via Adding the estimates for first and the second error term (4.12) we obtain the desired estimate (4.11) with the choices K * = K 2 + K 4 and V * = max{V * 2 , V * 4 }.
For consistency of notation we remark that E V ( ) can be rewritten as provided that this integral exists. The limit functional E is given by where we use that E is a continuous and non-negative function, so that E can be defined everywhere but attains the value +∞ if does not decay suitably at infinity. We will use the following semi-continuity result.
The following result gives the -convergence of E V to E with respect to the sequential weak* convergence as well as the equi-coercivity. (a) Compactness / equi-coercivity: (c) Limsup estimate / recovery sequence: where we may take Proof Obviously it is sufficient to show the lower bound (a) and the liminf estimate (b) for the smaller functional E V , and for = ρ dc with ρ ∈ L 1 (C) (resp. V = ρ V dc with ρ V ∈ L 1 (C)). We use the elementary convexity estimate ∀ r ≥ 0, a, w > 0 : wλ B (r /w) = r log(r /w) − r + w ≥ r log(a/w) − a + w.
We choose r (c) = (c), w(c) = W V (c), and a(c) = e −|c| 1 = I i=1 e −c i > 0. Note that a ∈ L ∞ (C) ∩ P(C) and W V /a is bounded from above, for any fixed V . Hence, c → log(a(c)/W V (c)) = −|c| 1 + V E V (c) is bounded from below, and we can integrate the above estimate to obtain the lower bound Since there exists a constant K 1 > 0 such that |c| 1 ≤ K 1 1+E(c) and since E V satisfies the lower bound in (4.11), we obtain the lower bound To show part (c) we use the indicated recovery sequence and the upper bounds for E V from (4.11). For a given ∈ P(C) we define V = ι V (κ V ( )). For an arbitrary continuous and bounded test function ψ we define the piecewise constant approximation ψ V via averaging over A V n . We obtain

This immediately implies (4.18) in part (a) with
where the convergence follows via Lebesgue's dominated convergence from the pointwise convergence ψ V → ψ and the uniform boundedness of ψ V . Thus, we conclude V * .

To show convergence of E V ( V ) it suffices to prove the upper bound lim sup
For this we use the bound ρ V (c) ≤ V I = 1/vol(A V n ) and the fact that ρ V and E V are constant on the same cubes to obtain where now only the measure is left. The first term tends to 0 for V → ∞, and the second can be estimated from above using the upper estimate in (4.11), which yields This implies the desired upper bound for V → ∞, and the proof is complete.

A Liminf Estimate for the Dual Dissipation Functional
Here we provide the liminf estimate for the dual dissipation potential * We observe that the latter term is linear in while the former term is convex in u V . Indeed, introducing the convex function G(a, b) = (a−b)(log a − log b) for a, b > 0 and noting the relation (a, b)(log a − log b) To establish the linear lower bound we use the elementary, affine lower bound This estimate follows easily by convexity, G(a, b) ≥ G(e ω , 1) + DG(e ω , 1) · (a−e ω , b−1), and 1-homogeneity giving G(e ω , 1) = DG(e ω , 1) · (e ω , 1). Note that equality holds in (4.23) if ω = log(a/b). Moreover, we have g(ω) + g(−ω) = 2 − e ω − e −ω ≤ 0, so a careful choice of ω depending on n will be necessary to obtain a good lower bound with a positive leading term.

Proposition 4.5 We have the liminf estimate
Proof The special forms of K(c), E(c), and * Lio in (4.21) give the formula * Since * V and * Lio are defined as sums over r = 1, . . . , R of nonnegative terms, it suffices to show the result for each r separately, where we suppress the index r .
Inserting (4.23) into (4.22) yields, with ω n ∈ R to be fixed afterwards, * where we used the detailed-balance conditions from Theorem 3.1 for the last identity. Rearranging the sum and recalling that A δ (n) = 0 for n / ∈ N we find * We now choose ω n = log A α V (n)/A β V (n) for n ∈ N and ω n = 0 otherwise and find, for all n with n ≥ α or n ≥ β, the relation The idea is now that as 1 V n → c > 0 we have the convergences To be more precise we define, for all ε ∈ ]0, 1[, the functions G ε (a, b) = −ε + min{(1−ε)G(a, b), 1/ε}, which converge monotonely to G(a, b) for ε 0. A lengthy calculation using the explicit structure of A δ (n) shows that for all ε > 0 there exists V ε 1 such that Hence, using the definition of ι V we find the lower bound * Since H ε is lower semi-continuous and bounded, this implies the liminf estimate Because ε > 0 was arbitrary we can use the monotone convergence H ε (c)

A Liminf Estimate for the Dissipation Functional
In the evolutionary -convergence method of [41,54,58] it is standard to provide a liminf estimate for the primal dissipation potential V which in our case is defined via the Legendre transform However, as our theory relies on the dualization V (u, v) ≥ n∈N u n ξ n − * V (u, ξ ) it will be sufficient to have the following limsup estimate for * V , which crucially relies on the concavity of the map (a, b) → (a, b).

Proposition 4.6 Consider any pair ( , ξ ) ∈ P(C)×C 1 c (C) and set ξ
. Then, for every family (u V ) V >1 we have the limsup estimate

c).(4.25)
Proof As in the proof of Proposition 4.5 we can exploit that * V is a sum of non-negative terms over r = 1, . . . , R. Hence, it is sufficient to show the desired limsup estimate for each reaction individually. For notational simplicity we drop the reaction index r . Defining and the functions a V , b V , and M ξ V are given Using where a ∞ (c) = c α /c α * , b ∞ (c) = c β /c β * , and γ = α−β. Using (r , t) ≤ 1 2 (r +t), the uniform boundedness of a V and b V on C R , and that V is a probability measure, we see that in the limsup of * V (u V , ξ V ) we can replace M ξ V by (α−β) · ∇ξ 2 without changing the limsup in the left-hand side of (4.25). Next we consider the functionals F : Note that f (r , t) ≥ 1 2 (r +t). Moreover, f is convex and positively homogeneous of degree 1. Thus, F is weak* lower semi-continuous on M (C R )×M (C R ), cf. [20,Thm. 6.57]. Now using the convergences Thus, in the view of the identity and observing that the first term on the right-hand side is weak * continuous, the limsup for This is the desired result for one reaction, and the full result follows by summation over

Convergence of Solutions
Here we provide the general convergence result as V → ∞ for the appropriately embedded Then, for all t > 0, we have the convergence For the proof we use the energy-dissipation principle for V ≥ 1 and pass to the limit in each of the terms. If u V is a solution of the CME, then for all T > 0 we have Following the ideas in [11] for the passage from a Markov chain to the Fokker-Planck equation or the general methods in evolutionary -convergence, we want to pass to the limit in each of the four terms. As a general fact, it will be sufficient to obtain liminf estimates on the left-hand side, since by a chain-rule argument an estimate with "≤" instead of equality can be turned back into an equality. Moreover, by the assumptions of the theorem we see that the right-hand side converges to the desired limit. However, it is rather delicate to pass to the limit in the integral because the potential V is only implicitly defined and we expect the limit to be given in terms of the Benamou-Brenier formula for the Wasserstein distance induced by the metric on (C, K). A major difficulty is even to obtain a suitable equi-continuity for the solutions u V to be able to extract a subsequence converging at all times. In particular, it is unclear how to pass to the limit in ι V (u V (t)) by a direct argument.
Hence, following [35], we estimate the primal dissipation potential V from below using the definition in terms of the Legendre transform of * V . Using additionally an integration by parts we have where u, η := n∈N u n η n . With this argument we can replace the energy-dissipation principle (4.28) by the estimate (4.29) which holds for all differentiable η. In this equation we are then able to pass to the limit V → ∞, when choosing η = η V = ι * V (ξ ) for a smooth function ξ . At the end we are then able to calculate the supremum over all ξ by using the especially simple quadratic structure in ξ , which mirrors the fact that the Liouville equation is a simple transport equation.

Proof of Theorem 4.7
Step 1: Embedding and uniform a priori bounds We now consider the family u V : [0, T ] → P(N ) and embed it into P(C) via ι V from (4.6). As in [11] we show an equi-continuity in a 1-Wasserstein distance, but introduce an additional weight accounting for our unbounded domain C. We define the maximal order p of all reactions via p := max{ |α r | 1 , |β r | 1 | r = 1, . . . , R }.
For μ ∈ M (C) and for 0 , 1 ∈ P(C) we set Using the definition of the Markov generators Q r V in terms of the coefficients B δ r V (n) , see (3.3), it is easy to derive the uniform estimate ι V (u V (t)) 1W ≤ C 1W independently of the initial conditions and V ≥ 1 (one simply needs u V n ≡ 1). Hence, we obtain the uniform Lipschitz bound (4.30) Step 2: Extraction of a subsequence The subset of P(C) defined by the boundedness of the above first moment is a compact subset of the metric space (P(C), d 1W ). Indeed, using Prokhorov's theorem one finds that this set is weak * sequentially compact. Since d 1W is dominated by the bounded Lipschitz metric (which metrizes weak * convergence), the compactness of (P(C), d 1W ) follows. Hence, we can apply the abstract Arzelà-Ascoli theorem in (P(C), d 1W ) to extract a subsequence V k → ∞ and a limit function At first, in place of (4.31a) one obtains d 1W ι V (u V (t)), (t) → 0. To derive (4.31a), we use the bound (4.30) together with the fact that any bounded continuous function can be uniformly approximated on compact sets by (multiples of) functions in F. Similarly, (4.31d) follows from (4.31b). In particular, combining (4.31d) and the assumption ι V (u V (0)) * 0 we conclude (0) = 0 . Finally, (4.31c) follows via (4.31a) from Theorem 4.4: Step 3: Limit passage in (4.29) Combining (4.31a) for t = T and Theorem 4.4 (cf. (4.19)), the first term satisfies the liminf estimate lim inf V →∞ E V (u V (T )) ≥ E( (T )). For the last term we use the assumption For the third term we employ Proposition 4.5 for each t ∈ [0, T ] based on (4.31a). Using Fatou's lemma we conclude the liminf estimate Thus, it remains to pass to the limit in J V (u V , η). For this we choose an arbitrary ξ ∈ C 1 c ([0, T ]×C) and define ξ V (t) = ι * V (ξ(t)), cf. (4.7). With this choice we can apply Proposition 4.6 for all t ∈ [0, T ] based on (4.31a). Now, Fatou's lemma yields In summary, we conclude that the limit function : [0, T ] → P(C) satisfies Step 4: Energy balance By inserting ξ ≡ 0 in (4.32) we obtain the upper bound We want to show energy balance, i.e., equality when the factor 2 is omitted. For this purpose, we observe that the measures (t, ·) ∈ P(C) decay at infinity such that (4.31c) holds. Hence, we may also use ξ(t, c) = λE(c) as testfunctions in (4.32). Writing shortly e(t) Maximizing with respect to λ leads to Taking the limit δ 0 and replacing ϕ by −ϕ, we obtain the desired result (4.26). With this, Theorem 4.7 is established.

Approximation via Fokker-Planck Equations
In the above section we have seen that the Liouville equation is the proper limit of the CME for V → ∞. However, for finite but large V it can still be advantageous to replace the discrete CME by a continuous PDE with V as a large parameter. In this range the stochastic modeling is done by the so-called Langevin dynamics, see [24,34,61], which is based on a stochastic perturbation of the reaction-rate equation (RRE), see (1.4). At the level of probability distributions the corresponding model is the associated Fokker-Planck equation (FPE). We will discuss two different gradient flow approximations: in the first we simply add a suitable "entropic term" to the driving functional, but keep the dissipation fixed (cf. Sect. 5.2), while in the second we expand E V and K V such that all terms of order 1/V are correct (cf. Sect. 5.3).

Improved Approximation of the Relative Entropy
We interpret the sum in the definition of E V as a Riemann sum and replace it by a corresponding integral. The main point of the improvement is that we keep the entropy term 1 V u n log u n in the definition of E V (u), which is in contrast to the limit E obtained in Theorem 4.4. Working with absolutely continuous probability measures (dc) = ρ(c) dc with ρ ∈ L 1 (C), we can define the V -dependent entropy by where the equilibrium density W V ∈ L 1 (C) has to be chosen suitably. A first simple However, a better and more refined W V is obtained using the next order of expansion in Stirling's formula (4.10) as well. For this we use the approximation k n ≈ n + 1/6, i.e., log(n!) = n log n − n + 1 2 log 2π(n+ 1 6 ) + O(1/n 2 ) for n → ∞. Hence, taking the limits V , |n| → ∞ such that n V → c, we obtain ) for E. We now take a probability measure = ρdc ∈ P(C) and a discrete approximation u ≈ κ V ( ) ∈ P(N ), where κ V : P(C) → P(N ) is the natural projection defined in (4.8). Then the Riemann-sum approximation results in The probability density W V is then defined by normalizing W (V , ·, c * ). We thus set Z(V , c * ) := ∞ 0 W(V , c, c * ) dc and This yields the expansion In summary, for E V defined via (5.1) and (5.2) we have and

Simple Fokker-Planck Approximation
Here we keep the V -independent Onsager operator K ( ) : ξ → −div K∇ξ of the Liouville equation and obtain the V -dependent continuous gradient system (P(C), ..,I . We expect that this FPE is a good approximation to the CME for all sufficiently large V . In particular, (5.4) has the steady state ρ = W V , which is close to the discrete steady state w V ∈ P(N ) using the embedding as above. In contrast, the only steady states of the Liouville equation (4.1) are concentrated on the equilibria ofċ = −R(c). Of course, the FPE still respects the invariant sets I(q), because the mobility K of the Onsager operator K is the same as for the Liouville equation. In particular, ρ = W V is the unique equilibrium density if and only if K has full rank, i.e., I(q) = C for all q ∈ Q.
The simpler choice W V (c) = 1 Z (V ) e −V E(c) for the equilibrium yields the relative entropy The flow equation˙ = −K ( )D E V ( ) induced by the gradient system (P(C), which is the same as (5.4) but with A V ≡ 0. The simplified equation will be used below as well, since W V has a simpler explicit form. We believe that this approximation is suitable for many purposes. However, it does not produce the correct diffusion as derived in [34,Eq. (1.7)]. This diffusion correction is used to improve the RREċ = −R(c) by replacing it by a stochastic differential equation called the chemical Langevin equations (CLE) in [24,61], see (1.4). The associated Fokker-Planck equation takes the forṁ where K CLE (c) ∈ R I ×I is given in (1.6) and differs from K as the logarithmic mean (a, b) between a = c α r /c α r * and b = c β r /c β r * is replaced by the arithmetic mean 1 2 (a+b). Obviously, (5.6) does not have a gradient structure with respect to K CLE , because there is no function c → E(c) such that R(c) = K CLE (c)∇ E(c).

Fokker-Planck Equation with Higher-Order Terms
To derive a proper expansion for the term of order 1/V in the evolution equation, we work with the V -dependent entropy E V defined in Sect. 5.1. Up to an irrelevant V -dependent constant, this functional approximates E V from (3.7) up to order 1/V 2 .
Similarly, we need to derive a suitable expansion for the dissipation potential, which can be done for each reaction independently. The discrete dual dissipation potential is given by (4.5b), namely * For a smooth function ξ : C → R we use the second-order accurate midpoint approximation where we used symmetric difference quotients to obtain second order accuracy. Moreover, for a smooth and sufficiently fast decaying = ρdc ∈ P(C) we define the associated discrete u ∈ P(N ) via u = κ V ( ) = ι * V , which yields V I u n = ρ 1 and similarly for V I u n+β . Hence, for the arguments of we find the expansion For all smooth functions f , g : C → R with compact support in int(C), the trapezoidal rule for Riemann integrals gives Hence, for smooth ρ and ξ we find the expansion * where the correction term ϒ takes the explicit form Here we used the relation ∂ a (a, b) = (a,b) a a− (a,b) giving  a∂ a (a, b)+b∂ b (a, b) = (a, b) and   a∂ a (a, b)−b∂ b (a, b) = (a, b) a+b Now we are in the position to calculate the first-order correction to the Liouville equation from the approximate entropy E V (cf. (5.3)) and the dual dissipation potential * V , namelẏ where the coefficients are given by It is interesting to see the cancellation in the term b 1 , where ϒ 1 did not have a sign, but after multiplication with (α−β)·∇ E(c) it becomes positive and increases the logarithmic mean (k fw c α , k bw c β ) to the arithmetic mean 1 2 (k fw c α +k bw c β ). Moreover, the coefficient b 0 consists of two terms, the first of which corresponds (up to order 1/V 2 ) to the correction A V in (5.4) arising from the improvement of E V , while the second term arises from improving the dissipation potential * V , namely via ϒ 0 . Putting these derivations together, summing over r = 1, . . . , R different reactions, and dropping all terms of order 1/V 2 , we find the following approximative Fokker-Planck equation:ρ , and K CLE is given in (1.6). The big disadvantage of equation (5.7) is that it is generally no longer a gradient system. However, it may be considered as an equation with an asymptotic gradient flow structure in the sense of [6]. To find the simplest true gradient system that is compatible with the Fokker-Planck equation (5.7), we have to find a true dual dissipation potential * V that is non-negative and coincides with * V from above to lowest order. To keep the notation light, we again explain the construction for the case of one reaction only and set 0 (c) = (k fw c α , k bw c β ). Our simplest choice is * where the higher-order corrections ϒ 2 (c) and ϒ 3 (c) need to be chosen such that * V (ρ, ξ ) is still coercive. Choosing θ 1 , θ 2 ∈ ]0, 1[ with θ 1 < θ 2 , we may require for ϒ 2 (c) and ϒ 3 (c) hold for the following choices (or any bigger ones) Of course, we fix the energy functional to be the improved entropy functional E V from (5.3), and the gradient system (P(C), E V , * V ) has the associated gradient-flow equatioṅ Because ∇ γ E V 1 is of order 1/V , we see that this equation involves terms up to order 1/V 4 , namely through a V 0 and through a V 2 /V 2 . Clearly, our gradient-flow equation (5.8) is much more complicated than those generated by the asymptotic gradient-flow structures in the sense of [6], where higher order terms are simply dropped.
There is also the question of well-posedness for equation (5.8). To have parabolicity of the leading terms we need that the mapping p → 1 V a V 1 p + 1 V 2 a V 2 p 2 + 1 V 3 a V 3 p 3 is monotone, which amounts to asking that a V 1 + 2 a V 2 q + 3 a V 3 q 2 ≥ 0 for all q ∈ R. This can be always be achieved by making ϒ 2 very big while keeping ϒ 3 constant, since ϒ 2 only enters once via a V 1 .

Comparison of Models
To appreciate the positive and negative aspects of the different approximations of the CME, we treat the simplest example, namely the linear RRE on C = [0, ∞[: Obviously, we have the explicit solution c(t) = 1 + (c(0)−1 e −t . The associated CME for u = P(N 0 ) is given bẏ Moreover, it can be easily checked that for any solution t → c(t) of the RRE (5.9) the following formula provides an explicit solution of the CME (5.10): Note that this is expression is compatible with the ODEs (5.11) for the moments, since for these Poisson distributions we have e(t) = c(t) and v(t) = c(t)/V . The Liouville equation and the simple Fokker-Planck equation read The Fokker-Planck equation for the chemical Langevin equation (cf. (5.6)) takes the form To compare the solutions of (FP) and (FP CLE ) with the true solutions of the CME (5.10), we assume that the solutions can be approximated by Gaußians. In general, for multidimensional Fokker-Planck equations of the formρ = 1 with a(t) ∈ R d and A(t) ∈ R d×d spd leads to the necessary conditionṡ see [55] for rigorous results of this type. Applying these formulas to (FP CLE ) we obtaiṅ a = 1 − a andȦ = −2 A + 1 + a, (5.13) hence the ODEs for a and A/V coincide with those for e and v in (5.11).
A similar argument indicates that solutions to (FP) are well approximated by Gaußians with mean a V and variance A V satisfyinġ On the one hand, this clearly indicates that (FP CLE ) provides a better approximation to the CME for t ∈ [0, T ]. By formally passing to the limit V → ∞ in (5.14), we see that the ODE for a V is asymptotically correct. This is not the case for the ODE for A V , since the arithmetic mean in (5.13) is replaced by the logarithmic mean in (5.14). However, the error of (1, c) compared to 1 2 (1+c) is less than 10 % for c ∈ [1/3, 3] and it converges to 0 for c → 1, i.e., in the limit t → ∞. Equations (5.13) are consistent with Kurtz' central limit theorem, which asserts that the normalized process 1 V N V (t) has fluctuations around c(t) of order 1/ √ V , and the rescaled process converges to a Gaußian process t → V (t) with covariance matrix A satisfyingȦ(t) = −DR(c(t))A(t) − A(t)DR(c(t)) T + 2 K CLE (c(t)), see, e.g., [34,Eq. (1.9)].
On the other hand, (FP) makes a better prediction for the equilibrium distribution that is attained for t → ∞. For (FP CLE ) we have the unique steady state Thus, E grows only like c, such that ρ eq,CLE V decays exponentially only. In contrast, the equilibrium ρ eq,FP V = Z −1 V e −V E(c) of (FP) produces the correct super-exponential decay of the stationary Poisson distribution equation for the CME (5.10).

Approximation via Cosh-Type Gradient Structure
The derivation of a gradient structure (4.3) for the Liouville equation (4.1) can be repeated very similarly by starting from the cosh-type gradient structure introduced in [44], see Proposition 3.8. We do not give the details here but provide the result only.
Starting from the cosh-type dual dissipation potential * cosh,V defined in (3.8) instead of the quadratic dual potential * V defined in (4.5) we obtain the counterparts to Propositions 4.5 and 4.6 but now with * Without any need to justify the approximation procedure in the sense of Sect. 4.2 we easily obtain the following result.
..,I , and where γ r = α r −β r . Using Using the identities √ ab (C * ) log(b/a) = b − a and √ ab (C * ) log(b/a) = (a+b)/2 the FP equation has the expansioṅ where K CLE is exactly the same as obtained in (1.6) by a completely different approach.

Hybrid Models
We show in this section how the different gradient structures for RRE, for CME, and for the FPE can be combined to obtain hybrid models, which are combinations of several models depending on the desired accuracy. The importance here is to use the proper rescalings in terms of the volume V to make the different descriptions compatible. We do not consider a full theory, but highlight first the general strategy of model reduction for gradient systems in Sect. 6.1 and then illustrate this by a simple example in Sect. 6.2. A nontrivial case of a rigorous coarse graining in this spirit is given in [43,47], where a linear RRE with a small parameter ε is considered. The elimination of the fast relaxations in the time scale ε leads to a coarse-grained gradient system. In Sect. 6.3 we discuss the general coupling of the FPE to a RRE and the similar coupling of the CME to a RRE, both leading to so-called mean-field equations, where a linear equation for a probability density is nonlinearly coupled to an ODE. Finally, we discuss the mixed discrete and continuous description, where the CME is used for small numbers of particles and the FPE is used for larger numbers.

Coarse Graining for Gradient Systems
If a gradient system (X, E X , X ) is more complicated than what is needed, one is interested in approximating the system by a simpler model that still contains the most important features. We explain how this can be done while keeping the gradient structure.
We assume that the relevant states x ∈ X can be described by states y ∈ Y and that there is a reconstruction mapping x = (y), i.e., (Y) is a subset (or submanifold) of X. We now pull back the gradient structure (X, E X , X ) to an approximative gradient structure (Y, E Y , Y ). The natural approach is to restrict the energy functional and the (primal) dissipation potential as follows: E Y (y) = E X ( (y)) and Y (y,ẏ) := X ( (y), D (y)ẏ). (6.1) The solutions y : [0, T ] → Y of the coarse-grained gradient system (Y, E Y , Y ) will provide good approximations x : t → (y(t)) ∈ X of the true solutions of the full GS (X, E X , X ), if the set (Y) approximates a flow-invariant subset of X.
In reaction systems, the primal dissipation potential X is usually not known explicitly. Hence, it is desirable to have a method for reducing the dual dissipation potential * X directly to * Y , in the case where A = D (y) : Y → X is injective but its adjoint mapping A * : X * → Y * has a large kernel. The following exact result will be the motivation for our modeling approximations in the subsequent subsections.
Proof For the proof we use the saddle-point theory in [14,Ch. VI.2].
This yields (6.2), since the right-hand side is clearly infinite as well.
Fix now η ∈ Ran(A * ) and define the Lagrangian function L : For notational convenience we set Classical duality theory yields the trivial inequality P ≥ D. Clearly, L(·, ξ) is convex and lower semicontinuous, whereas L(x, ·) is concave and upper semicontinuous, since the boundedness of A * implies that ξ ∈ X * A * ξ = η is closed. Using η ∈ Ran(A * ), we find ξ η ∈ X * with A * ξ η = η, so that our assumptions guarantee the coercivity of x → L(x, ξ η ) ∈ R. Hence, we can apply [14, Chap. VI, Prop. 2.3], which shows that there is no duality gap: We relate P and D with the two sides in our desired formula (6.2). On the one hand, where in the last step we have substituted ξ = ξ η +ζ with A * ξ η = η and introduced δ 0 ( η) = 0 for η = 0 and ∞ otherwise. Thus, we conclude h(x) = (Ay) − η, y for x = Ay and h(x) = ∞ for x / ∈ Ran(A).
Thus, taking the minimum over all of X is the same as taking it over Ran(A), namely On the other hand, the definition of g(ξ ) = inf x∈X L(x, ξ) immediately gives g(ξ ) = − * (ξ ) − χ * (ξ ). Hence, we arrive at As a result, formula (6.2) follows from P = D.
In our applications below (as well as in most others) the explicit minimization in (6.2) is too complicated to be executed. However, as the coarse-graining mapping through is usually only an approximation, it may suffice to approximate the minimizers suitably. In general, one has to find an approximation ξ = M(y, η) ∈ X * and sets * Y (y, η) = * X ( (y), M(y, η)) or 1 2 η, K Y (y)η = L(y)η, K X ( (y))L(y)η , (6.4) where L(y) : Y * → X * is the linear version of M. Of course, when constructing M or L one should keep (6.2) in mind to preserve all interesting properties inherited by the coarsegraining process.
Remark 6.2 (Coarse-graining for fast-slow reaction systems) Particular cases of rigorous coarse graining are discussed in [43] and [47], where linear and nonlinear fast-slow RRE are considered, respectively. However, the approach is somehow opposite: if ε > 0 is the time scale of the fast reactions, then in the limit ε → 0 a linear constraint is imposed on DE(c) defining a slow manifold c = (y), and c is constraint to lie on this linear or nonlinear manifold. The reduced or coarse-grained gradient system is then given by E Y (y) = E( (c)) and * Y (y, η) = * slow ( (y),

A Simple Example: from CME to RRE
We apply the above idea with (X, E X , * X ) being (P(N ), In the simple exampleċ = 1−c treated in Sect. 5.4 the image of V defines an exactly invariant submanifold, but this is no longer true for nonlinear equations or systems. Nevertheless our construction provides the surprising identity with the old E defined in (2.9) which is independent of V .
To reduce the dual dissipation potential X defined via K V we use the derivative Thus, the adjoint operator D V (c) * maps μ = (μ n ) to ζ = (ζ i ) i=1,...,I via In general, one is not able to solve the minimization problem (6.2) that produces * Y from * X , so instead we construct a linear mapping ζ → μ = M V (c)ζ that approximates the minimizer for V → ∞ and satisfies ζ = D V (c) * M V (c)ζ . Indeed, we search for μ in the linear form μ a n = a · n for n ∈ N and obtain where we used the identities n∈N e −V |c| 1 (V c) n n! n i = V c i and n∈N e −V |c| 1 (V c) n n! n i n j = V 2 c i c j + δ i j V c i . Thus, we choose the simple operator M V of the form For inserting μ = M V (c)ζ and u = V (c) into the full dual dissipation potential * X , we use the form (4.5b) and the relations ( With this, we find an approximation of the reduced dual dissipation potential * Y , namely * Thus, the gradient system (Y, E Y , * Y ) obtained by the abstract reduction procedure is exactly given by (C, E, K), which is the gradient system for the RRE (2.2) studied in Theorem 2.2.

Coupling a RRE to a Fokker-Planck Equation
In many applications one is interested in the microscopic description of some variables c j , while other variables c i can be described more macroscopically. We first start from the simplified FPE (5.5) as the gradient system (P (C . It can be shown that A(v) ≥ 1 and e V ( a, a * ) ≥ a * λ B ( a/a * ) for all V , and for V → ∞ we obtain e V ( a, a * ) → a * λ B ( a/a * ). To simplify the model we are therefore allowed to replace the last term in E V Y by the relative entropy E m ( c m ) = I j=J +1 c * j λ B ( c j /c * j ) for the RRE. Neglecting the irrelevant constant term Z (V )/V , we obtain the hybrid energy again as a relative entropy, namely For the Onsager operator we also use a cruder reduction than the minimization advocated in Sect. 6.1. We simply postulate the Onsager operator K V via the dual dissipation potential * V ,FP-RR ( s , c m ; ξ, ζ ) = where ξ ∈ C 1 (C s ) and ζ ∈ R I −J . Indeed, in the sense of the general reduction method explained in Sect. 6.1 we see that K FP-RR This reveals that the system is a classical mean-field model, which is linear in the density ρ for the component c s while it is nonlinearly coupled to the ODE for the component c m .

Coupling a RRE to a CME
In analogy to the coupling of an RRE for some macroscopic c m to a Fokker-Planck equation we can directly couple the CME to an RRE, which leads to hybrid system defined on Instead of given the general derivation as in Sect. 6.3, we just give an explicit example.

Combining CME and Fokker-Planck Descriptions
We consider the simplest nontrivial model, namely the scalar RREċ = a − bc with a, b > 0, which is induced by the reaction ∅ b − − a X . This corresponds to α = 0, β = 1, k fw = a, and k bw = b. We have the following three derived gradient systems: (1) The RREċ = a − bc is generated by the gradient system (R + , K, E) with steady state c * = a/b, K(c) = (a, bc), and E(c) = a b λ B (bc/a). (2) The associated chemical master equationu = B V u is generated by the gradient system (P(N 0 ), E V , K V ) and readṡ u n = V au n−1 − V a+bn u n + b(n+1)u n+1 for n ∈ N 0 (with u −1 = 0) (6.7) and has the steady state w V = (e −V a/b (V a/b) n /n!) n∈N 0 . The entropy and Onsager operator are where Z V ,N is uniquely determined by asking N dw V ,N = 1. The entropy functional is defined via the obvious relative entropy per volume, namely where du dw denotes the Radon-Nikodym derivative. The difficult part is the modeling of the Onsager operator K V ,N (u) as it includes the crucial transfer between the discrete and the continuous parts of the hybrid model. We define K in terms of its associated quadratic form acting on smooth functions ξ : N → R, where we write ξ n for ξ(n) and W (c) for W V (c): While the first and the third terms on the right-hand side give the purely discrete and the continuous parts of the state space, respectively, we see that the second term is the new term that couples the discrete and the continuous parts. The parameter a is still to be chosen, the natural parameter being a. The evolution equation for u is again a linear equation of the formu = B V ,N u, i.e., it corresponds to a continuous-time Markov process. It consists of a discrete part, as in (6.7) but only for n = 0, . . . , N − 2, and a continuous part, as in (6.8) but only for c > N /V . The new structure is the coupling between the two subsystems which gives rise to the following conditions: By our definition of w V ,N we have w N −1 W (N /V ) ≈ N b/(aV 2 ) and see that for a = a these conditions take the approximate forṁ where the second relation clearly shows the corresponding Robin boundary condition connecting the parabolic Fokker-Planck equation to the discrete system on {0, . . . , N −1}. Note that u n and U are scaled such that V u N −1 is comparable to U (N /V ).