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 Γ -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 to , which may represent different molecules, e.g., , , and . We assume that these I species undergo R different reactions of mass-action type:where the vectors contain the stoichiometric coefficients, and are the forward and backward reaction rates, see Sect. 2. The reaction would lead to the vectors and .

Denoting by the vector of nonnegative densities, the simplest model is the macroscopic (RRE), which is a system of ODEs on the state space :Here the monomials indicate that the probability for the right number of particles for the rth 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 for each species and consider the random vector . A forward or backward reaction of type r is modeled as an instantaneous event where the particle numbers jump from to or vice versa. The corresponding jump rates in a volume of size are given by and respectively; see (3.1) for the definition of .

Here we study the vector of probabilitiesthat describes the probability distribution of the random variable . The time evolution of is given by the (CME), i.e., the Kolmogorov forward equation associated with the continuous time Markov chain above. This is a countable linear system of ODEs:where is an (unbounded) linear operator on , 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 such thatDefining the Boltzmann entropy E and the Onsager operator viawhere is the logarithmic mean, we see that (RRE) is generated by the gradient system , namely . In Sect. 2.5 we also discuss further gradient structures, e.g. those used in [43, 44, 46, 47].

If (RRE) satisfies the detailed-balance condition, then (CME) does so with an equilibrium distribution that is explicitly given as a product of one-dimensional Poisson distributions with mean , namely (cf. Theorem 3.1),Consequently, we are also able to interpret (CME) as a gradient flow induced by a gradient system , see (3.7). Here is again the Boltzmann entropy with respect to , but now divided by the volume V:

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-33] by connecting the particle numbers to the concentrations and showing thatwhere is the solution of (RRE) with . 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 , which leads to a rate functional that generates a gradient structure ; 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 for the gradient system , 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 (Sect. 4) and obtain the coarse grained gradient system withIn particular, the coarse grained gradient flow equation is the associated with (RRE); here we used that and . 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 we obtain . While Kurtz works directly on the Markovian random variables, we work at the level of their distributions:

Our convergence result for the gradient systems to the limiting gradient system 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 it is also important to study the case of intermediate V, where already behaves continuously but still shows some fluctuations of standard deviation , see [61] for a numerical approach to treat the hierarchy via a suitable hybrid method. In [34] it is shown that the random vector obtained by solving the stochastic differential equation(see [34, Eq.  (1.7)]) yields an improved approximation because , while . 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 (CLE).

The associated Kolmogorov forward equation takes the form Here the diffusion matrix can be written in the explicit formthat is different from , 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 withwhich leads to the linear In Sect. 5 we show that by systematically deriving higher-order corrections to and we can recover the asymptotically correct diffusion matrix 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 by a simpler one via an embedding mapping . 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: . 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 . The resulting gradient flow equation turns out to be a mean-field equation, in which the density of the component satisfies a linear equation which is nonlinearly coupled to an ODE for the component .

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 . Finally we analyze a mixed CME / Fokker–Planck model with state space , in which the underlying space 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–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.The space of all signed Borel measures of bounded variation on is denoted by .

Reaction-rate equation The RRE is denoted by :

state and state space , steady state , dual variable

energy functional , Onsager operator

conserved quantities , stoichiometric subsets .

Chemical master equation The CME is denoted by :

state and state space , steady state , dual variable

energy functional , Onsager operator

invariant subsets .

Liouville equation The LE is denoted by :

state and state space , steady state , dual variable

energy functional , Onsager operator

Fokker–Planck equation The FPE is denoted by :

state and state space , steady state , dual variable

energy functional , Onsager operator .

Hybrid systems are denoted by “mathfrak” letters:

for coupling FPE and RRE

for coupling CME and RRE

for merging discrete and continuous modeling for one species.

Reaction Rate Equations

We denote by the concentrations of I different chemical species reacting according to the mass action law, i.e., the reactionsfor , where R is the number of possible reactions, are the vectors of the stoichiometric coefficients, and are the forward and backward reaction-rates. In general these rates may depend on , 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 systemwhere , see [17, 19, 26].

Stoichiometry, Conservation, and Decomposition of the State Space

The stoichiometric subspace and its orthogonal complement are defined viaFor each the function defines a first integral, which easily follows from . These conservation laws often go under the name conservation of atomic species, see [17]. Suppose now that is a non-trivial subspace of . We shall argue that the RRE induces a decomposition of the state space into affine invariant subsets. (If , the only invariant set is itself.)

Choosing a basis of we define the matrix , which has the rows . By construction we have , and we conclude that the solutions of (2.2) conserve as follows:By construction every affine conserved quantity is of the form for some and . This allows us to decompose the full state space into the invariant, affine subsets for . Using the notationwe define, for all , the setsThen, implies , and we have . Let us note that this decomposition does not depend on the choice of the orthonormal basis which determines the matrix , although the set does depend on . Note also that we can always write for some arbitrary satisfying .

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 such that all reactions are simultaneously in equilibrium, i.e.,This condition implies that , 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 be the matrix which has the row vectors , . We call the Wegscheider matrix because of the pioneering work in [60]. We then havewhich explains the abbreviation . Since is strictly positive, we can take the logarithm of the polynomial conditions (2.6) and find the equivalent linear systemBy Fredholm’s alternative, (2.7) is solvable if and only ifThese conditions on the reaction coefficients and are called Wegscheider conditions (see, e.g., [25, 57, 59, 60]). By choosing a basis of and exponentiation they can be rewritten as polynomial conditions without referring to the equilibrium state .

Let denote the number of Wegscheider conditions. Then the following assertions hold: Since by standard linear algebra, the number of Wegscheider conditions can be expressed asHence, if the number R of reactions is smaller than the number I of species, the Wegscheider conditions can usually be satisfied easily.

If the stoichiometric vectors , , are linearly independent, then , hence there is no Wegscheider condition.

If , , are linearly dependent, then and non-trivial Wegscheider conditions appear.

Remark 2.1

(Wellposedness of RRE) We conclude this subsection with a statement concerning the well-posedness of the RRE given as in Theorem 2.2 below. For all there exists a unique global solution . Local existence for solutions starting in the interior of is trivial, as 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 for some .

Moreover, solutions cannot leave this region via the boundary , since the vector field is either tangential to or points inwards. Indeed, if for some j, thenbecause each term in the sum is nonpositive: If or , then the term is 0. Thus, we are left with the cases for some positive n. In the first case implies 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 . Here, the state space contains all possible concentration vectors . The driving functional is the relative entropy and the Onsager matrix is chosen suitably (recall that ):where the logarithmic-mean function is given viaThe following result shows that a RRE (2.2) satisfying the detailed-balance condition (2.6) is indeed generated by the gradient system . This was first established in [62, Sect. VII] to derive entropy bounds for hyperbolic conservation laws in reactive flows and was rederived in [40] in the context of reaction diffusion systems including electric charge-interactions. It is interesting to note that for continuous time Markov chains (CTMC), which form a special subclass of RRE with linear reactions, there are several distinct gradient structures, see [37, Prop. 4.2] and [42, Thm. 3.1] and Sect. 2.4. However, in the case of nonlinear reactions according to the mass-action law, only the gradient structure with the Boltzmann entropy remains. The key fact is the logarithm identity .

Theorem 2.2

(Gradient structure for RRE) If the RRE (2.2) satisfies the detailed-balance condition (2.6) for a positive steady state , then it has the gradient structure defined in (2.9), namely .

Proof

Multiplying by we obtainwhich is the denominator of . Hence, using giveswhere we used the detailed-balance condition (2.6) in . Thus, the assertion is established.

Summarizing the above derivations, we have rewritten the RRE in thermodynamic formwhich 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 and the thermodynamic driving force , see e.g. [10, Ch. X, § 4]. The symmetry is related to microscopic reversibility, i.e., detailed balance, see also [44, 45]. Subsequently, we refer to as the Onsager operator or matrix.

Here we clearly see the advantage of using the Onsager operator 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 is not invertible except if .

Continuous Time Markov Chains as a Gradient System

The forward equation for a reversible CTMC on a discrete space is a special case of the RRE considered above. In this case all reactions are of the formand the reaction rates (resp. ) are interpreted as the transition rates from i to j (resp. from j to i). The reaction-rate equation is the linear system of ODEsand the detailed-balance condition for the equilibrium state takes the formUsing this condition, the RRE can be written coordinate-wise asHere we used the notational convention that for . The relative entropy E is as above and the Onsager matrix takes the formwhere denotes the i-th unit vector. We then have the gradient structure , namelyThis 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 we do not need the summation rule for logarithms. Hence, following [37, 42] there are more general gradient structures. Choosing a strictly convex function that is smooth on we set where for and . The gradient flow structure corresponds to the case where .

Proposition 2.3

(Gradient structure for CTMC) If the CTMC (2.13) satisfies the detailed-balance condition (2.14) for a positive steady state , then it has the gradient structures , namely .

Remark 2.4

The construction in Proposition 2.3 does not extend to general RRE. There one would need to replace the quantity in (2.9) by , 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 and ), by virtue of the convexity of .

In the following we will mainly concentrate on the gradient structure with the logarithmic entropy, as it is the only one that connects with the RRE.

Generalized Gradient Structures

For Markov chains and RRE there are several families of generalized gradient structures where the quadratic function is replaced by a general dual dissipation potential that is continuous and convex and satisfies .

In the case of RRE, the monomial terms can only be generated by the logarithmic summation rule . Hence, we stick to the relative entropy E defined in (2.9), i.e., . However, we may replace the linear Onsager principle by the more general nonlinear form .

To define we choose an arbitrary family of smooth dissipation functionals , i.e., and and define the dissipation potentialUsing (2.11) we easily obtain , i.e., is generated by the generalized gradient system .

The case leads to the quadratic dissipation potential in (2.9), i.e., the functions are given in terms of the logarithmic mean. In [3] the choices is used. Based on a derivation via the large deviation principle (see [44-46]) a special role is played by the choice of a “cosh-type” function :Here is normalized such that . Hence, the dual dissipation potential takes the formIt 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].

The Chemical Master Equation

Modeling Discrete Particle Numbers via CME

The chemical master equation (CME) is a CTMC that is defined on the set where is the vector of particle numbers, see [39] for an introduction. This means that denotes the number of particles of species in a sufficiently big volume, whose size is denoted by . 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 probabilitiesAll 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 and , respectively. Moreover, the jump intensitiesalso depend on the volume V, as denotes the absolute particle number, while for the reaction the densities matter. The specific form of (cf. [32, 39]) readsTo avoid clumsy notation we defined for all , but if . We also see that for , 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 probability distributions , namelyThe rth forward reaction from to can only happen (i.e., ) if . Hence any occurring with is multiplied by intensity 0, so in (3.2) we may set for all . The operators are the adjoints of the Markov generators given byfor .

We emphasize that the RRE as well as the CME are uniquely specified if the reaction network (2.1), the reaction rates and , and the volume are given. Hence, there are obviously close relations between both models, in particular for , see [4, 5, 23, 32, 61].

So far, we have not used the detailed-balance condition, i.e., we can even allow for in the above considerations. In all cases, the Kolmogorov forward equation is an infinite-dimensional 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 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 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.

Theorem 3.1

(Detailed balance for CME) Let be given in the form (3.1). Assume that (2.2) has the equilibrium satisfying the detailed-balance condition (2.6). Then the equilibrium given bysatisfies the detailed-balance condition for the CME (3.2), namely

For each reaction we obtain the relationAnalogously we obtain the same result for , where the detailed-balance condition (2.6) is used in the definition of .

Using the detailed-balance coefficients we can rewrite the operator from (3.2) in a symmetrically balanced form aswhere is the unit vector, i.e., .

It is important to realize that in general the steady state for the detailed-balance condition is highly non-unique, because of the discrete versionsof the invariant stoichiometric subspaces . Indeed, choosing arbitrary and defining via for and 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.

Example 3.2

(Equation without detailed balance) For we consider the RREwhich consists of two individual reaction pairs, namely and with the individual steady states and . The joint steady state of (3.5) is , and we have detailed balance if and only if .

Building the CME according to (3.2) based on the two reaction pairs we obtainFor the case and , where the detailed-balance condition holds with , we obtainand it is easy to check that is a steady state.

However, for and the detailed-balance condition fails with . The CME readsAn explicit calculation shows that the Poisson distribution based on , i.e., , is not a steady state. Indeed, inserting into the right-hand side of the last equation we find (for )

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 and produce the RRE that has the equilibrium . However the CME readsAgain, the Poisson distribution with is not the equilibrium:Note that the reversible reaction pair 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 we construct a special Green’s function . General Markov chain theory (e.g., [36, Ch. 2]) implies that there exist a unique minimal solution to the backward equationassociated with the CME with initial condition . This minimal solution is non-negative and satisfies and , 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 spacesas well as the weighted spaceswith the corresponding norms and the usual modification for . Now, we consider the transition semigroup defined bywhich 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 .

Theorem 3.4

Assume that the detailed balance condition (2.6) holds. Then, the semigroup extends to a -semigroup of contractions on for all . Moreover, the semigroup is selfadjoint on and Markovian, i.e., for all .

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.

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 . Each of the irreducible components of the state space is then one-dimensional (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 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 for all , and at least one component is strictly positive. Then each of the infinite irreducible components of is of the formfor some , and the restricted Markov process is a birth-death process with birth rate and death rate given byReuter’s criterion [53, Thm. 11] gives a characterization of non-explosion for birth-death chains; it asserts that the chain is non-explosive if and only ifIn our setting we have , so that , and thereforeSince the summands tend to as , we infer that the latter sum is infinite; hence the Markov chain is non-explosive, or equivalently (see [36, Thm. 2.33]).

The case of multiple reactions follows by induction on the number of reactions R. Indeed, for , let denote the semigroup corresponding to the reactions . Then the Trotter product formula for contraction semigroups on (see e.g., [9]) asserts thatstrongly in . 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 . We also observe that the class of finitely supported functions is a core for each of the generators. The Markovianity of thus follows from the Trotter formula and the Markovianity of and .

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 to . Denoting by the generator of the -semigroup on , we define the operator byThis definition of is consistent with the explicit formula for given above. Since generates a -semigroup of contractions on , it follows that generates a -semigroup of contractions on . Furthermore, since preserves positivity and , it follows that is invariant under the semigroup generated by .

As an immediate consequence we obtain global well-posedness for the CME in .

Theorem 3.6

(Global well-posedness of the CME) Let the detailed-balance condition (2.6) hold. Then, for all there exists a unique mild solution to the CME (3.2) satisfying .

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 that is smooth on , let us writewhere is defined after (2.16), is given in Theorem 3.1, and denotes the -th unit vector in . The following result is then a special case of Proposition 2.3.

Proposition 3.7

(Quadratic gradient structures for CME) If the RRE (2.2) satisfies the detailed-balance condition (2.6) for a positive steady state , then the associated CME has the gradient structure defined in (3.6), namely

In the following we will mainly be concerned with the case that is the logarithmic entropy, where is the Boltzmann function . In that case we obtainwhere the logarithmic mean is defined in (2.10). The above definitions do not only restrict to the entropy function , but also introduce a normalization with respect to the volume V. Hence, can be seen as an entropy per unit volume. The corresponding scaling of was chosen such that the evolution equation is the same as .

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 in (2.18) and note the special scaling via the volume V in (3.8) below, which is needed because 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 , then the associated CME has the gradient structure with from (3.7) and

The desired formula follows easily by recalling from (3.4) and by using and .

Liouville and Fokker–Planck Equations

For general evolutionary equations one can define a measure-valued flow in the phase space that is given by transporting the measures according to the semiflow of the original equation. The evolution equation describing this measure-valued flow is the Liouville equation. For our RRE in we assume that we have a global semiflow and consider probability measures that are obtained by transporting with , namelyIn particular, if , then .

It is now easy to see that satisfies the Liouville equationin the sense of distributions. We will regard (4.1) as an evolution equation in the space . 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 with .

The goal of this section is to give a rigorous connection between the CME for and the Liouville equation in terms of the associated gradient structures.

The Liouville Equation as a Gradient System

We show that the gradient structure 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 that acts on functions viawhere and are taken with respect to . We also consider the affine potential energy functional defined byIn the next result we identify the formal gradient structure for the Liouville equation.

Proposition 4.1

(Gradient structure for the Liouville equation) If the RRE (2.2) satisfies the detailed-balance condition (2.6) for a positive steady state , then the associated Liouville equation has the gradient structure , namely

Let and let be a signed measure of finite total variation such that and for |h| sufficiently small. Then we havehence for all . Therefore, The gradient flow equation 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 .

More precisely, we will show that after a suitable V-dependent embedding of into the proper scalings of the functionals and converge in the sense of -convergence to the corresponding structures for the Liouville equation given by the gradient system , 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 of solutions of the CME to the solution of the Liouville equation , thereby recovering Kurtz’ result (1.3), see Sect. 4.6.

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 solves the CME if and only if for all the following energy-dissipation estimate holds:where we use the quadratic dissipation potential and its Legendre dual defined via with from (3.7), namely where the second form uses Theorem 3.1 and is especially useful to perform the limit , see the proof of Proposition .

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 may not have densities; indeed, because we want to recover the Kurtz result (1.3) we are interested in the “deterministic case” . So our analysis has to be more careful in dealing with general limit measures. For this, we use the dualization approach introduced in [35] where is estimated from below by for suitably chosen recovery functions .

In order to compare probability measures on different spaces and , we consider a suitable embedding . Here is simply obtained by assigning the mass of at uniformly to the cubeMore explicitly, is given bywhere denotes the indicator function with for and 0 otherwise. The corresponding dual operation acting on functions is given byThe final convergence result will be formulated in Theorem 4.7, which will be a direct consequence of the following three estimateswhere the dual dissipation potential is defined viaWe will see in Sect. 4.6 that the limsup estimate for the dual potential in Sect. 4.5 provides a weak form of a liminf estimate for the primal potential .

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 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 and are in fact linear in .

-Limit of the Relative Entropies

We also define and and consider the functionalswhere is defined viaThese definitions are chosen such that for all .

Finally we define a natural inverse of , namelysuch that is a projection from onto .

To understand the limit of for we will use the representationIn Lemma 4.2 below we will show that converges pointwise to E as defined in (2.9). To quantify the latter convergence, we use the classical lower and upper bounds of [49] for Stirling’s formula:Using this estimate and recalling E from (2.9) we obtain the following estimate.

Lemma 4.2

(Pointwise bound for ) For all there exist and such that for all the following bounds hold:

We decompose the error viawith defined by . For the second term we use the convexity of and the estimate . Hence, we haveSumming this inequality over we obtain the upper boundFor the opposite direction we use (a) that decreases on [0, 1] and the convexity of which implies (b) for . This yieldsif , where and 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 identitywith from (4.10). Because of we obtain, for all , the lower boundFor the upper bound we use and for . Hence for we obtain, using again the estimate ,Summation over yields, for all and , the upper boundTogether with the lower estimate we control the first error term in (4.12) viawhere .

Adding the estimates for first and the second error term (4.12) we obtain the desired estimate (4.11) with the choices and .

For consistency of notation we remark that can be rewritten asprovided that this integral exists. The limit functional is given bywhere we use that E is a continuous and non-negative function, so that can be defined everywhere but attains the value if does not decay suitably at infinity. We will use the following semi-continuity result.

Lemma 4.3

(Lower semi-continuity of ) For sequences with , we have .

For cut-off functions with we have by weak* convergence and continuity of E. Using yields . Choosing a non-decreasing sequence with for all we have by Beppo Levi’s monotone convergence, and the assertion follows.

The following result gives the -convergence of to with respect to the sequential weak* convergence as well as the equi-coercivity.

Theorem 4.4

(-convergence of to ) Let and be defined on as above. Then we have the following properties:

Compactness / equi-coercivity:

Weak* liminf estimate:

Limsup estimate / recovery sequence: where we may take .

Obviously it is sufficient to show the lower bound (a) and the liminf estimate (b) for the smaller functional , and for with (resp. with ). We use the elementary convexity estimateWe choose , , and . Note that and is bounded from above, for any fixed V. Hence, is bounded from below, and we can integrate the above estimate to obtain the lower boundSince there exists a constant such that and since satisfies the lower bound in (4.11), we obtain the lower boundThis immediately implies (4.18) in part (a) with . Moreover, if then we have the lower bound and the liminf estimate (4.19) follows from Lemma 4.3.

To show part (c) we use the indicated recovery sequence and the upper bounds for from (4.11). For a given we define . For an arbitrary continuous and bounded test function we define the piecewise constant approximation via averaging over . We obtainwhere the convergence follows via Lebesgue’s dominated convergence from the pointwise convergence and the uniform boundedness of . Thus, we conclude .

To show convergence of it suffices to prove the upper bound . For this we use the bound and the fact that and are constant on the same cubes to obtainwhere now only the measure is left. The first term tends to 0 for , and the second can be estimated from above using the upper estimate in (4.11), which yieldsThis implies the desired upper bound for , and the proof is complete.

A Liminf Estimate for the Dual Dissipation Functional

Here we provide the liminf estimate for the dual dissipation potential based on the lower boundWe observe that the latter term is linear in while the former term is convex in . Indeed, introducing the convex function for and noting the relation we haveTo establish the linear lower bound we use the elementary, affine lower boundThis estimate follows easily by convexity, , and 1-homogeneity giving . Note that equality holds in (4.23) if . Moreover, we have , so a careful choice of depending on will be necessary to obtain a good lower bound with a positive leading term.

Proposition 4.5

We have the liminf estimate

The special forms of , , and in (4.21) give the formulaSince and are defined as sums over 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 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 for we findWe now choose for and otherwise and find, for all with or , the relationThe idea is now that as we have the convergenceswhich yields and as desired. To be more precise we define, for all , the functionswhich converge monotonely to G(a, b) for . A lengthy calculation using the explicit structure of shows that for all there exists such that for all and all . Even more, if we define the functions , then, for all there exists such thatHence, using the definition of we find the lower boundSince is lower semi-continuous and bounded, this implies the liminf estimateBecause was arbitrary we can use the monotone convergence to conclude the desired result for each of the R reactionsSummation over yields the full result for .

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 which in our case is defined via the Legendre transformHowever, as our theory relies on the dualization it will be sufficient to have the following limsup estimate for , which crucially relies on the concavity of the map .

Proposition 4.6

Consider any pair and set with defined in (4.7). Then, for every family we have the limsup estimate

As in the proof of Proposition 4.5 we can exploit that is a sum of non-negative terms over . Hence, it is sufficient to show the desired limsup estimate for each reaction individually. For notational simplicity we drop the reaction index r.

Defining , relation (4.5b) leads us to the integral representationwhereand the functions , , and are givenUsing there exists such that , and we have uniform convergencewhere , , and . Using , the uniform boundedness of and on , and that is a probability measure, we see that in the limsup of we can replace by without changing the limsup in the left-hand side of (4.25).

Next we consider the functionals given byNote that . Moreover, f is convex and positively homogeneous of degree 1. Thus, F is weak* lower semi-continuous on , cf. [20, Thm. 6.57]. Now using the convergenceswe obtain the liminf estimate .

Thus, in the view of the identityand observing that the first term on the right-hand side is weak continuous, the limsup for givesThis is the desired result for one reaction, and the full result follows by summation over and the definition of , namely .

Convergence of Solutions

Here we provide the general convergence result as for the appropriately embedded solutions of the CME to the solutions of the Liouville equation, which is a simple transport along the solutions of the RRE . Our approach follows the strategy of evolutionary -convergence as initiated in [54, 58] with the new idea of dualization as introduced in [35].

Theorem 4.7

(Evolutionary -convergence of CME to Liouville) For all consider a solution of the CME (3.2). Assume that the initial conditions are well-prepared in the sense thatThen, for all , we have the convergencewhere is the unique solution of the Liouville equation (4.3) starting at , i.e., for all with we haveMoreover, for all with we have the energy identity

For the proof we use the energy-dissipation principle for and pass to the limit in each of the terms. If is a solution of the CME, then for all we haveFollowing 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 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 . A major difficulty is even to obtain a suitable equi-continuity for the solutions to be able to extract a subsequence converging at all times. In particular, it is unclear how to pass to the limit in by a direct argument.

Hence, following [35], we estimate the primal dissipation potential from below using the definition in terms of the Legendre transform of . Using additionally an integration by parts we havewhere . With this argument we can replace the energy-dissipation principle (4.28) by the estimatewhich holds for all differentiable . In this equation we are then able to pass to the limit , when choosing 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 and embed it into via 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 . We define the maximal order p of all reactions viaFor and for we setwhere .

Using the definition of the Markov generators in terms of the coefficients , see (3.3), it is easy to derive the uniform estimate independently of the initial conditions and (one simply needs ). Hence, we obtain the uniform Lipschitz boundMoreover, as by well-preparedness, the equi-coercivity of established in (4.18) yields the uniform boundStep 2: Extraction of a subsequence The subset of defined by the boundedness of the above first moment is a compact subset of the metric space . Indeed, using Prokhorov’s theorem one finds that this set is weak sequentially compact. Since is dominated by the bounded Lipschitz metric (which metrizes weak convergence), the compactness of follows.

Hence, we can apply the abstract Arzelà-Ascoli theorem in to extract a subsequence and a limit function such that At first, in place of (4.31a) one obtains . 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 . Similarly, (4.31d) follows from (4.31b). In particular, combining (4.31d) and the assumption we conclude . Finally, (4.31c) follows via (4.31a) from Theorem 4.4:Step 3: Limit passage in (4.29) Combining (4.31a) for and Theorem 4.4 (cf. (4.19)), the first term satisfies the liminf estimate . For the last term we use the assumption .

For the third term we employ Proposition 4.5 for each based on (4.31a). Using Fatou’s lemma we conclude the liminf estimateThus, it remains to pass to the limit in . For this we choose an arbitrary and define , cf. (4.7). With this choice we can apply Proposition for all based on (4.31a). Now, Fatou’s lemma yieldsIn summary, we conclude that the limit function satisfiesfor all .

Step 4: Energy balance By inserting in (4.32) we obtain the upper boundWe want to show energy balance, i.e., equality when the factor 2 is omitted. For this purpose, we observe that the measures decay at infinity such that (4.31c) holds. Hence, we may also use as testfunctions in (4.32). Writing shortly we find and obtainMaximizing with respect to leads to which implies , or more explicitly , which is the desired energy balance (4.27) for and .

Moreover, we can repeat the calculation on [0, s] with instead of [0, T]. The full result (4.27) follows by subtracting the identity on [0, r] from that on [0, s].

Step 5: Weak form of gradient flow equation With Step 4 we rewrite (4.32) asand know that the left-hand side is maximized by . Inserting the test functions with small and with we arrive, after some cancellations and after dividing by , atTaking the limit 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 . 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 and 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 as a Riemann sum and replace it by a corresponding integral. The main point of the improvement is that we keep the entropy term in the definition of , which is in contrast to the limit obtained in Theorem 4.4. Working with absolutely continuous probability measures with , we can define the V-dependent entropy bywhere the equilibrium density has to be chosen suitably. A first simple approximation is with as above and . However, a better and more refined is obtained using the next order of expansion in Stirling’s formula (4.10) as well. For this we use the approximation , i.e., for . Hence, taking the limits such that , we obtainwith the V-dependent correction for E.

We now take a probability measure and a discrete approximation , where is the natural projection defined in (4.8). Then the Riemann-sum approximation results inThe probability density is then defined by normalizing . We thus set andThis yields the expansionwith . In summary, for defined via (5.1) and (5.2) we haveand .

Simple Fokker–Planck Approximation

Here we keep the V-independent Onsager operator of the Liouville equation and obtain the V-dependent continuous gradient system . The associated gradient-flow equation is the FPEwhere we used and set .

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 , which is close to the discrete steady state using the embedding as above. In contrast, the only steady states of the Liouville equation (4.1) are concentrated on the equilibria of . Of course, the FPE still respects the invariant sets , because the mobility of the Onsager operator is the same as for the Liouville equation. In particular, is the unique equilibrium density if and only if has full rank, i.e., for all .

The simpler choice for the equilibrium yields the relative entropyThe flow equation induced by the gradient system is the simplified FPEwhich is the same as (5.4) but with . The simplified equation will be used below as well, since 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 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 formwhere is given in (1.6) and differs from as the logarithmic mean between and is replaced by the arithmetic mean . Obviously, (5.6) does not have a gradient structure with respect to , because there is no function such that .

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 defined in Sect. 5.1. Up to an irrelevant V-dependent constant, this functional approximates from (3.7) up to order .

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), namelyFor a smooth function we use the second-order accurate midpoint approximation with to obtain the expansionwhere we used symmetric difference quotients to obtain second order accuracy. Moreover, for a smooth and sufficiently fast decaying we define the associated discrete via , which yields ,and similarly for . Hence, for the arguments of we find the expansionFor all smooth functions with compact support in int, the trapezoidal rule for Riemann integrals givesHence, for smooth and we find the expansionwhere the correction term takes the explicit formHere we used the relation , givingNow we are in the position to calculate the first-order correction to the Liouville equation from the approximate entropy (cf. (5.3)) and the dual dissipation potential , namely , which yieldswhere the coefficients are given byIt is interesting to see the cancellation in the term , where did not have a sign, but after multiplication with it becomes positive and increases the logarithmic mean to the arithmetic mean . Moreover, the coefficient consists of two terms, the first of which corresponds (up to order ) to the correction in (5.4) arising from the improvement of , while the second term arises from improving the dissipation potential , namely via .

Putting these derivations together, summing over different reactions, and dropping all terms of order , we find the following approximative Fokker–Planck equation:where , , and 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 that is non-negative and coincides with from above to lowest order. To keep the notation light, we again explain the construction for the case of one reaction only and set . Our simplest choice iswhere the higher-order corrections and need to be chosen such that is still coercive. Choosing with , we may requirefor all , so that . The bounds for and hold for the following choices (or any bigger ones)Of course, we fix the energy functional to be the improved entropy functional from (5.3), and the gradient system has the associated gradient-flow equation . With we findBecause is of order 1/V, we see that this equation involves terms up to order , namely through and through .

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 is monotone, which amounts to asking that for all . This can be always be achieved by making very big while keeping constant, since only enters once via .

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 :Obviously, we have the explicit solution .

The associated CME for is given bywhere . Using the linearity in (5.9), which leads to the linearity in n of the coefficients in (5.10), we obtain explicit closed form relations of the evolution of the rescaled expectation and variance , namelyMoreover, it can be easily checked that for any solution 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 and .

The Liouville equation and the simple Fokker–Planck equation readThe Fokker–Planck equation for the chemical Langevin equation (cf. (5.6)) takes the formTo compare the solutions of (FP) and (FP) 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 the ansatz with and leads to the necessary conditionssee [55] for rigorous results of this type. Applying these formulas to (FP) we obtainhence the ODEs for a and A/V coincide with those for and in (5.11).

A similar argument indicates that solutions to (FP) are well approximated by Gaußians with mean and variance satisfyingOn the one hand, this clearly indicates that (FP) provides a better approximation to the CME for . By formally passing to the limit in (5.14), we see that the ODE for is asymptotically correct. This is not the case for the ODE for , since the arithmetic mean in (5.13) is replaced by the logarithmic mean in (5.14). However, the error of compared to is less than 10 % for and it converges to 0 for , i.e., in the limit . Equations (5.13) are consistent with Kurtz’ central limit theorem, which asserts that the normalized process has fluctuations around of order , and the rescaled process converges to a Gaußian process with covariance matrix satisfying , see, e.g., [34, Eq.  (1.9)].

On the other hand, (FP) makes a better prediction for the equilibrium distribution that is attained for . For (FP) we have the unique steady stateThus, grows only like c, such that decays exponentially only. In contrast, the equilibrium 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 defined in (3.8) instead of the quadratic dual potential defined in (4.5) we obtain the counterparts to Propositions 4.5 and but now withWithout any need to justify the approximation procedure in the sense of Sect. 4.2 we easily obtain the following result.

Proposition 5.1

(Cosh-type gradient structure for the Liouville equation) The Liouville equation (4.1) has the gradient structure with from (4.2) and from above.

The result follows by using , , andwhere . Using and the definition of gives which is the desired right-hand side of (4.1) when testing with and integrating by parts.

As in the case of quadratic gradient structure for the Liouville equation we may consider the first-order correction to obtain a Fokker–Planck equation. For this we insert the improved energy defined in (5.3) into the dissipation potential (cf. (3.8)) to obtain a quasilinear Fokker–Planck-type equation, namely . Using the abbreviations and we find (note )Using the identities and the FP equation has the expansionwhere 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 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 can be described by states and that there is a reconstruction mapping , i.e., is a subset (or submanifold) of . We now pull back the gradient structure to an approximative gradient structure . The natural approach is to restrict the energy functional and the (primal) dissipation potential as follows:The solutions of the coarse-grained gradient system will provide good approximations of the true solutions of the full GS , if the set approximates a flow-invariant subset of .

In reaction systems, the primal dissipation potential is usually not known explicitly. Hence, it is desirable to have a method for reducing the dual dissipation potential directly to , in the case where is injective but its adjoint mapping has a large kernel. The following exact result will be the motivation for our modeling approximations in the subsequent subsections.

Proposition 6.1

Consider reflexive Banach spaces and and a real-valued dissipation potential (i.e. lower semicontinuous, convex, and ) that is superlinear, i.e. for . Assume that the bounded linear operator has closed range. Then the dissipation potential satisfieswhere we use the convention .

For the proof we use the saddle-point theory in [14, Ch. VI.2].

Fix and assume first that . Since is closed, the Closed Range Theorem yields that is closed as well, and . Consequently, there exists such that , and we obtainThis yields (6.2), since the right-hand side is clearly infinite as well.

Fix now and define the Lagrangian function viaFor notational convenience we setClassical duality theory yields the trivial inequality . Clearly, is convex and lower semicontinuous, whereas is concave and upper semicontinuous, since the boundedness of implies that is closed.

Using , we find with , so that our assumptions guarantee the coercivity of . 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 and introduced for and otherwise. Thus, we concludeThus, taking the minimum over all of is the same as taking it over , namelyOn the other hand, the definition of immediately gives . Hence, we arrive atAs a result, formula (6.2) follows from .

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 and setswhere is the linear version of . Of course, when constructing or one should keep (6.2) in mind to preserve all interesting properties inherited by the coarse-graining 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 is the time scale of the fast reactions, then in the limit a linear constraint is imposed on defining a slow manifold , and is constraint to lie on this linear or nonlinear manifold. The reduced or coarse-grained gradient system is then given by and , which means , see [43, Thm. 5.6] and [47, Prop. 4.2].

A Simple Example: from CME to RRE

We apply the above idea with being and with being . The embedding mapping is given by the Poisson distributionsIn the simple example treated in Sect. 5.4 the image of defines an exactly invariant submanifold, but this is no longer true for nonlinear equations or systems. Nevertheless our construction provides the surprising identitywith the old E defined in (2.9) which is independent of V.

To reduce the dual dissipation potential defined via we use the derivativeThus, the adjoint operator maps to viaIn general, one is not able to solve the minimization problem (6.2) that produces from , so instead we construct a linear mapping that approximates the minimizer for and satisfies . Indeed, we search for in the linear form for and obtainwhere we used the identities and . Thus, we choose the simple operator of the formFor inserting and into the full dual dissipation potential , we use the form (4.5b) and the relations andWith this, we find an approximation of the reduced dual dissipation potential , namelyThus, the gradient system obtained by the abstract reduction procedure is exactly given by , 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 , while other variables can be described more macroscopically. We first start from the simplified FPE (5.5) as the gradient system and partition the components of into stochastic and macroscopic parts, and respectively, viaIn the notation of Sect. 6.1 we let and .

For the mapping we choose the product ansatzwhere the probability densities are given as follows:According to Sect. 6.1 the functional on is then given bywith and .

It can be shown that and for all V, and for we obtain . To simplify the model we are therefore allowed to replace the last term in by the relative entropy for the RRE. Neglecting the irrelevant constant term , we obtain the hybrid energy again as a relative entropy, namelyFor the Onsager operator we also use a cruder reduction than the minimization advocated in Sect. 6.1. We simply postulate the Onsager operator via the dual dissipation potentialwhere and . Indeed, in the sense of the general reduction method explained in Sect. 6.1 we see that is obtained from by inserting and .

Thus, the hybrid model induced by the gradient system is given by the coupled system for and :It is interesting to see that the last terms can be rewritten in terms of the RRE , viz.This reveals that the system is a classical mean-field model, which is linear in the density for the component while it is nonlinearly coupled to the ODE for the component .

Coupling a RRE to a CME

In analogy to the coupling of an RRE for some macroscopic 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.

For we consider the simple reaction with stoichiometric vectors , , and . The associated system of RREs is given byWe have the conservation relation and the detailed-balance steady state . The associated CME on takes the formThe detailed-balance steady state by with . As in (4.5b) the full Onsager operator is defined viaWe partition , i.e., we keep in stochastic description via the distribution , while will be treated macroscopically. Thus, we define the gradient system with relative entropy and Onsager operator defined viafor and . Again, is obtained from by inserting and and performing an approximation for large V. The associated evolution equation is the hybrid systemClearly, this system is consistent with the conservation law , in the sense that satisfies .

Combining CME and Fokker–Planck Descriptions

We consider the simplest nontrivial model, namely the scalar RRE with , which is induced by the reaction . This corresponds to , , , and . We have the following three derived gradient systems: To combine the description via the CME and the Fokker–Planck equation we consider the mixed state space . Hence, counts the number of atoms, while for we use the concentration as a continuous variable to describe the state. A typical choice could be to be sure to capture all small discrete effects.

The RRE is generated by the gradient system with steady state , , and .

The associated chemical master equation is generated by the gradient system and reads and has the steady state . The entropy and Onsager operator are

The associated Fokker–Planck equation (5.4) takes the form This equation has the equilibrium solution (cf. (5.2)) and is generated by the gradient system with

The hybrid gradient system is described by measuresThe idea is now to choose and rather than to model the evolution equation.

We first choose the equilibrium state in the formwhere is uniquely determined by asking . The entropy functional is defined via the obvious relative entropy per volume, namelywhere denotes the Radon–Nikodym derivative.

The difficult part is the modeling of the Onsager operator as it includes the crucial transfer between the discrete and the continuous parts of the hybrid model. We define in terms of its associated quadratic form acting on smooth functions , where we write for and W(c) for :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 is still to be chosen, the natural parameter being a.

The evolution equation for is again a linear equation of the form , i.e., it corresponds to a continuous-time Markov process. It consists of a discrete part, as in (6.7) but only for , and a continuous part, as in (6.8) but only for . The new structure is the coupling between the two subsystems which gives rise to the following conditions:By our definition of we have and see that for these conditions take the approximate formwhere the second relation clearly shows the corresponding Robin boundary condition connecting the parabolic Fokker–Planck equation to the discrete system on . Note that and U are scaled such that is comparable to U(N/V).