Convergence to equilibrium of renormalised solutions to nonlinear chemical reaction-diffusion systems

The convergence to equilibrium for renormalised solutions to nonlinear reaction-diffusion systems is studied. The considered reaction-diffusion systems arise from chemical reaction networks with mass action kinetics and satisfy the complex balanced condition. By applying the so-called entropy method, we show that if the system does not have boundary equilibria, then any renormalised solution converges exponentially to the complex balanced equilibrium with a rate, which can be computed explicitly up to a finite dimensional inequality. This inequality is proven via a contradiction argument and thus not explicitly. An explicit method of proof, however, is provided for a specific application modelling a reversible enzyme reaction by exploiting the specific structure of the conservation laws. Our approach is also useful to study the trend to equilibrium for systems possessing boundary equilibria. More precisely, to show the convergence to equilibrium for systems with boundary equilibria, we establish a sufficient condition in terms of a modified finite dimensional inequality along trajectories of the system. By assuming this condition, which roughly means that the system produces too much entropy to stay close to a boundary equilibrium for infinite time, the entropy method shows exponential convergence to equilibrium for renormalised solutions to complex balanced systems with boundary equilibria.


Introduction and Main results
The large time behaviour of reaction-diffusion systems is a long standing and yet highly active topic in the analysis of partial differential equations. Classical methods include dynamical systems, invariant regions or linearisation methods. Recently, the so-called entropy method, a fully nonlinear approach, proved to be very useful in studying convergence to equilibrium for many PDE systems, in particular reaction-diffusion systems which feature a suitable dissipative structure.
The aim of the present paper is to prove convergence to equilibrium for nonlinear reaction-diffusion systems arising from chemical reaction networks. A chemical reaction network is a quadruple {S, C, R, K} in which S = {S 1 , . . . , S N } denotes the set of chemical substances, C = {y 1 , . . . , y |C| } with y i ∈ ({0} ∪ [1, ∞)) N , i = 1, . . . , |C| is the set of chemical complexes, which are either reactants and/or products of a chemical reaction and y = (y j ) N j=1 denotes a vector of stoichiometric coefficients for the substances S 1 , . . . , S N , where y j = 0 iff the substance S j is part of the complex y. Correspondingly, R is the set of all considered chemical reactions y r → y ′ r with y r , y ′ r ∈ C and for r = 1, . . . , |R|. Moreover, K = {k r : r = 1, . . . , |R|} is set of the associated reaction rate constants with k r > 0 being the rate of the reaction y r → y ′ r for all r = 1, . . . , |R|. The reaction network {S, C, R, K} is assumed to satisfy the following natural conditions: (1) for each S i ∈ S, there exists at least one complex y ∈ C for which the corresponding stoichiometric coefficient y i is nontrivial, i.e. y i ≥ 1, (2) there exist no trivial reaction y → y ∈ R for any complex y ∈ C, (3) for any y ∈ C, there must exist a y ′ ∈ C such that y → y ′ ∈ R or y ′ → y ∈ R, i.e. every complex must be either reactant or product of at least one reaction.
In the following, we shall use the convention that the primed complexes y ′ r ∈ C (respectively y ′ ∈ C) denote the product of the r-th reaction while the unprimed complexes y r ∈ C (respectively y ∈ C) denote the reactant; except when specified otherwise.
As an example for a reaction network, we find for the single reversible reaction To specify a reaction-diffusion system modelling a chemical reaction network {S, C, R, K}, we assume the reactions to take place in a bounded vessel (or reactor) Ω ⊂ R n , where Ω is a bounded domain with Lipschitz boundary. We also assume (w.l.o.g. after a suitable rescaling of the space variable) that Ω has normalised volume, i.e. |Ω| = 1.
We denote by c(x, t) = (c 1 (x, t), . . . , c N (x, t)) the vector of concentrations where c i (x, t) is the concentration of S i at time t > 0 and position x ∈ Ω. Each substance S i is assumed to diffuse in Ω with a strictly positive diffusion coefficient d i > 0. The corresponding reaction-diffusion system then reads as where the diffusion matrix D = diag(d 1 , . . . , d N ) is positive definite since d i > 0 for all i = 1, . . . , N , and R(c) represents all the reactions in R. We shall apply the law of mass action to get an explicit form of R(c) which reads as where k r > 0 denotes the reaction rate constant of the r-th reaction. Finally, system (1) is subject to nonnegative initial data c 0 (x) ≥ 0 (by which we mean c 0 (x) := (c 1,0 (x), .., c N,0 (x)) and c i,0 (x) ≥ 0 for i = 1, .., N and x ∈ Ω), and homogeneous Neumann boundary conditions c(0, x) = c 0 (x) for x ∈ Ω, and ∇c · ν = 0 for (x, t) ∈ ∂Ω × R + , where ν := ν(x) is the outward normal unit vector at point x ∈ ∂Ω.
Therefore it follows from (1) that ∂ t (Q c) − QD∆c = Q R(c) = 0 and hence due to the homogeneous Neumann boundary condition, the co-dimension of the Wegscheider's matrix W leads to m (linearly independent) mass conservation laws of the following form where c = (c 1 , . . . , c N ) and c i = Ω c i dx (after recalling that |Ω| = 1), and M is called an initial mass vector, which depends on the choice of Q. By changing the signs of some rows of Q if necessary, we can always consider (w.l.o.g.) a matrix Q such that the initial mass vector M is non-negative, i.e. M ∈ R m + . To state the main results of this paper, we need the following definitions concerning equilibria of chemical reaction networks.
• c ∞ is called a detailed balanced equilibrium if for each forward reaction y in R, there exists in R also the corresponding backward reaction y ′ k b − → y (with k b > 0) and • c ∞ is called a complex balanced equilibrium if the total outflow and inflow at the equilibrium c ∞ are equal for every complex y ∈ C, i.e. for any complex y ∈ C, we have where {s : y ′ s = y} denotes the set of all reactions y s ks>0 − −− → y ′ s with fixed product complex y ′ s = y ∈ C. • c ∞ is called a boundary detailed/complex balanced equilibrium (or shortly a boundary equilibrium) if c ∞ is a detailed/complex balanced equilibrium and c ∞ ∈ ∂R N + . • It follows directly from the above definitions that c ∞ is a detailed balanced equilibrium =⇒ c ∞ is a complex balanced equilibrium =⇒ c ∞ is an equilibrium, but the reverse is in general not true.
A chemical reaction network is called complex balanced if for each strictly positive mass vector M ∈ R m >0 it possesses a strictly positive (i.e. not a boundary) complex balanced equilibrium.
The concept of detailed balance goes back as far as Boltzmann for modelling collisions in kinetic gas theory and for proving the H-theorem for Boltzmann's equation [4]. It was then applied to chemical kinetics by Wegscheider [58]. The complex balanced condition was also considered by Boltzmann [5] under the name semi-detailed balanced condition or cyclic balanced condition, and was systematically used by Horn, Jackson and Feinberg in the seventies for chemical reaction network theory, see e.g. [21,39,41].
It is well-known that if a chemical reaction network (such as modelled by system (1)-(3)) has one complex balanced equilibrium, then all other possible equilibria (independently of the initial mass vector) are necessarily also complex balanced, see e.g. [39,41]. Moreover, for every fixed positive initial mass vector M ∈ R m >0 , there exists a unique complex balanced equilibrium c ∞ ∈ R N >0 satisfying the mass conservation laws determined by the initial mass vector M ∈ R m >0 . Note that (possibly infinitely) many boundary equilibria may exist as well.
Throughout this paper, we will refer to this strictly positive equilibrium as the complex balanced equilibrium while all other equilibria are simply called boundary equilibria. Moreover, we will consider positive initial mass vectors M ∈ R m >0 in order to ensure that any considered complex balanced network features a positive complex balanced equilibrium c ∞ ∈ R N >0 . Note that all our results hold equally true for non-negative initial mass vectors M ∈ R m + as long as there exists a unique positive complex balanced equilibrium c ∞ ∈ R N >0 , which will typically (but not always) be the case.
This paper aims to prove exponential convergence to equilibrium of solutions to the nonlinear reactiondiffusion system (1)-(3) under the assumption the considered chemical reaction network is complex balanced.
The method of proof is the so-called entropy method. The main idea of the entropy method is to qualitatively exploit the decay of a suitable entropy (e.g. convex Lyapunov) functional E[f ] along a trajectory f of an evolution process: where D[f ] is called entropy production functional or also entropy dissipation functional in cases when E[f ] is physically an energy functional. The latter is the case for nonlinear complex balanced reactiondiffusion systems of the form (1)-(3), where the following logarithmic relative free energy functional constitutes a suitable entropy functional. Note that for nonlinear reaction-diffusion systems, the above logarithmic relative entropy is the only generally existing Lyapunov functional, while for linear complex balanced systems, other generalised relative entropy functional do also exist, see e.g. [25,47]. The following explicit form of the entropy dissipation functional D(c) associated to (7) along the flow of system (1)-(3) was derived in [15]: The entropy method applies to general evolution processes, which are well behaved in the sense that This condition holds true for the system (1)- (3), where D(c) = 0 is satisfied by all constant states which balance the reactions of the complex balanced network. Thus, provided no boundary equilibria exist, taking into account all conservation laws uniquely identifies the complex balanced equilibrium c ∞ . Given such a well-behaved evolution process, the entropy method aims to quantify the decay of the entropy functional E[f ] in terms of the relative entropy towards the equilibrium state. More precisely, the goal is an entropy-entropy production estimate (which is a functional inequality independent of the flow of the evolution process) of the form where Φ(x) ≥ 0 and Φ(x) = 0 ⇔ x = 0. More specifically for system (1)-(3), the first key result of this paper is to prove the following entropy-entropy dissipation estimate (or rather free energy-free energy dissipation estimate) for some constant λ > 0. Assuming that such a functional inequality is proven and that a suitable concept of solutions to system (1)-(3) satisfyies a weak entropy-entropy dissipation law (i.e. an integrated version of the formal relation (8)) of the form for almost all 0 ≤ s < t, then a Gronwall argument implies exponential convergence to equilibrium first in relative entropy, i.e.
and consequently in L 1 -norm, thanks to a Csiszár-Kullback-Pinsker type inequality, see Lemma 2.2 below.
In the first main results of this paper, we prove for general, complex balanced reaction-diffusion systems (1)-(3) without boundary equilibria, that any so-called renormalised solution (which is the only existing solution concept for such a general class of nonlinear reaction-diffusion systems, see Theorem 2.1) converges exponentially to the complex balanced equilibrium with a rate which can be explicitly estimated in terms of the systems' parameters and a constant obtained from a finite dimensional inequality with mass conservation constraints. More precisely, our first main theorem reads as Theorem 1.1 (Convergence to equilibrium for general complex balanced reaction-diffusion systems without boundary equilibria).
Let Ω be a bounded domain in R n with Lipschitz boundary ∂Ω. Assume that the diffusion matrix D is positive definite, i.e. d i > 0 for all i = 1, . . . , N . Moreover, we assume that system (1)-(2) is complex balanced. Consequently, for each positive initial mass vector M ∈ R m >0 there exists a unique positive complex balanced equilibrium c ∞ ∈ R N >0 . Assume in addition that system (1)-(2) does not have boundary equilibria.
Then, for all states c ∈ R N >0 satisfying E(c|c ∞ ) < +∞ and Q c = M , there exists a constant H 1 > 0 depending only on Q, the stoichiometric coefficients y ∈ C, M and E(c|c ∞ ) such that Here Further, inequality (11) implies that for all measurable vector functions c : Ω → R N + satisfying E(c|c ∞ ) ≤ K and Q c = M , the entropy-entropy dissipation inequality (9), i.e.
..,N {d i } and C LSI is the constant in the Logarithmic Sobolev inequality (see Lemma 2.4), and K 1 and K 2 are constants (see (18) and (29)) depending explicitly on the domain Ω, the diffusion matrix D, the stoichiometric coefficients y ∈ C, the reaction rate constants k r , the initial mass M , the complex balanced equilibrium c ∞ and the constant K.
Finally, as a consequence of functional inequality (9), any renormalised solution c(x, t) to (1)-(3) (see Theorem 2.1) associated with initial data c 0 satisfying Q c 0 = M and E(c 0 |c ∞ ) < +∞, converges exponentially to c ∞ in L 1 -norm with the rate λ/2, that is where C CKP is the constant in a Csiszár-Kullback-Pinsker type inequality (see Lemma 2.2).
Remark 1. Note that while an explicit bound for H 1 in (11) can certainly be obtained near the equilibrium c ∞ via Taylor expansion, such bounds far from equilibrium are highly nontrivial and an open problem due to the non-convexity of the involved nonlinear terms. Moreover, an additionally difficulty stems from the lack of a constructive approach to characterise and exploit the matrix Q.
Theorem 1.1 comprises in our opinion the most general equilibration result for complex balanced reaction-diffusion systems, which is currently feasible. It generalises previous results on the exponential convergence to equilibrium for reaction-diffusion systems, partially in terms of considering complex balanced instead of detailed balanced systems, partially in terms of applying to renormalised solutions rather than weak-or classical solutions, and partially that the obtained convergence rate λ is explicitly stated in terms of the key constant H 1 .
At this point, we review some previous results concerning the large time behaviour of reaction-diffusion systems arising from chemical reaction networks: • The first results on the entropy methods for nonlinear reaction-diffusion systems trace back to works of Gröger, Glitzky and Hünlich [30,31,32,33,34], where the authors consider electrochemical drift-diffusion-recombination models. However, the proof of associated entropy-entropy dissipation estimate was based on a contradiction argument in combination with a compactness method, thus provided only convergence to equilibrium in space dimension two and without explicit control of the rate of convergence.
• For detailed balanced systems without boundary equilibria, a first general approach to prove exponential convergence to equilibrium for (1)-(3) was presented in [44]. The inspired key idea of [44] was to prove an entropy-entropy dissipation estimate via a suitable convexification argument (of the non-convex sum of reaction terms in (8)). The disadvantage, however, is that except in special cases (e.g. 2S 1 ⇌ S 2 ) the convexification argument seems not to allow for explicit estimates on the rate of convergence. The results of [44] were extended in [15] to complex balanced systems thanks to the derivation of the entropy dissipation (8).
• In a recent work [26], we proposed a constructive approach to show exponential convergence to equilibrium for general detailed balanced reaction-diffusion systems, which allows to obtain explicit bounds on the rates of convergence in contrast to the convexification argument of [44]. The applicability of the constructive approach was demonstrated for two typical example systems: i) a reversible reaction of arbitrary many chemical substances α 1 S 1 + · · · + α I S I ⇌ β 1 B 1 + · · · + β J B J ( * ) and ii) a reversible enzyme reaction S 1 + S 2 ⇌ S 3 ⇌ S 4 + S 5 . This approach is also applicable to complex balanced systems as demonstrated in [15] for a cyclic reaction α 1 S 1 → α 2 S 2 → · · · → α N S N → α 1 S 1 . Also in [26], we provided an If-Theorem that for any detailed balanced systems, under the assumption of a finite dimensional inequality (like (11)) and a technical non-degeneracy assumption on the entropy dissipation, then the solutions converge exponentially to the positive equilibrium with explicit rates. In this paper, we are able to remove these technical assumptions as well as generalise the result to complex balanced systems. It is also worth mentioning that the reversible reaction with arbitrary chemical substances ( * ) was also recently treated in the paper [51]. Altogether, these previous results prove either exponential convergence for general systems at the price of a lack of explicitness of convergence rates, or they showed explicit rates of convergence for some special classes of reaction-diffusion systems.
The results of Theorem 1.1 improve the previous results in several directions: i. We prove the functional inequality (9) explicitly up to the finite dimensional inequality (11). More precisely, Theorem 1.1 states that the constant λ in (9) scales with the minimum of λ 1 (derived from the diffusion coefficients and the Logarithmic Sobolev Inequality) and the constant H 1 from (11) times the structural constant K 2 /K 1 with K 1 and K 2 given in (18) and (29). We note that the idea of proving (9) by using a finite dimensional version was already considered in [44]. However, the approach therein lacks explicitness due to the use of the convexification argument. ii. We provide a general result of exponential convergence to equilibrium for complex balanced systems without boundary equilibria. In particular, the rate of convergence is explicitly controlled in terms the constant H 1 of the finite dimensional inequality (11) (and other explicit parameters). It is emphasised that although the constant H 1 is not explicit in general, we believe it is possible to explicitly estimate H 1 in any concrete system once the mass conservation laws are explicitly known (see Section 2.2 for such a system arising from reversible enzyme reactions). iii. Another important advantage of Theorem 1.1 and our method of proof is its role in a potential strategy to consider systems with boundary equilibria. This leads to the second main result of this paper, which is discussed in the following paragraphs.
It is important to point out that the entropy-entropy dissipation inequality (9), and consequently the finite dimensional inequality (11), cannot hold for general systems with boundary equilibria: If a solution trajectory of such a system should approach a boundary equilibrium, then the entropy dissipation D(c) tends to zero while the relative entropy to the complex balanced equilibrium E(c|c ∞ ) remains positive, see e.g. [15,26] for the details. Consequently, an entropy-entropy dissipation estimate of the form (9) cannot hold.
This structural difficulty is already encountered in complex balanced reaction networks in the ODE setting, i.e. by considering the solution u(t), which satisfies the ODE system where R(u) is defined as (2) with u in place of c. There is an extensive literature concerning the large time asymptotics of complex balanced systems of the form (13). Indeed, it is proven that the unique strictly positive complex balanced equilibrium of an ODE reaction network is locally stable (cf. [41]). Moreover, it is conjectured that the positive complex balanced equilibrium is in fact globally stable, i.e. it is the unique global attractor for the dynamical system given by the ODE network (with the exception of initial data starting on ∂R N + ). This statement is usually called the Global Attractor Conjecture (GAC) and has remained one of the most important open problems in the theory of chemical reaction networks, see e.g. [1,8,36,45] and the references therein. A recently proposed proof of this conjecture in the ODE setting is currently under verification [9].
For reaction-diffusion systems of the form (1)-(3), it was pointed out in [15,Remark 3.6] that if the boundary equilibria are unstable in the sense that solution trajectories cannot stay too close to those equilibria (in L 1 -norm distance) for too long, then the convergence to the complex balanced equilibrium follows via a contradiction argument. However, proving such an instability for boundary equilibria is usually a subtle issue, in particular in the PDE setting (1)- (3).
In this paper, by using elements of the proof of Theorem 1.1, we establish a weaker condition entailing instability of boundary equilibria and convergence to the complex balanced equilibrium. More precisely, our condition is based on a quantitative estimate that solution trajectories do not converge to a boundary equilibrium "too fast" (if it should converge at all), see Theorem 1.2. To explain this approach further, we remark at first that our proof of deriving the entropy-entropy dissipation inequality (9) from the finite dimensional inequality (11) is independent of the presence of boundary equilibria. Thus, instead of trying (or rather failing) to prove (11) as a pure functional inequality, we look for a generalisation with a time-dependent coefficient H 1 (t) along the trajectories of solutions, where H 1 (t) may tend to zero in case a solution trajectory would converge to a boundary equilibrium. Therefore, we look for a modified entropy-entropy dissipation inequality along solutions c(x, t) of (1)-(3) of the following form (which is no longer a pure functional inequality like (9)) with λ(t) = 1 2 min{λ 1 , K 2 H 1 (t)/K 1 } where K 1 and K 2 are given in (18) and (29). Intuitively, the time dependent function λ(t) (which may decay to 0 as t → ∞) gives a lower bound for the entropy dissipation D(c(t)) or equivalently for the convergence of a trajectory towards a boundary equilibrium (where E(c(t)|c ∞ ) remains bounded below). Therefore, if λ(t) satisfies +∞ 0 λ(s)ds = +∞ or equivalently the function H 1 (t) satisfies +∞ 0 H 1 (s)ds = +∞, then it follows from Gronwall's inequality and the weak entropy-entropy dissipation law (10) that By using this (so far non-exponential) convergence to the complex balanced equilibrium c ∞ , we obtain the L 1 -instability of the boundary equilibria. In return, this instability allows to show an entropy-entropy dissipation estimate of the form (9) on a reduced domain of states, which is strictly bounded away from the boundary equilibria. Thus, we recover exponential convergence to the complex balanced equilibrium c ∞ after a sufficiently large time. Our second main result reads as follows: Theorem 1.2 (Conditional convergence to equilibrium for complex balanced reaction-diffusion systems with boundary equilibria).
Let Ω be a bounded domain in R n with Lipschitz boundary ∂Ω. Assume that the diffusion matrix D is positive definite, i.e. d i > 0 for all i = 1, . . . , N . Moreover, we assume that system (1)-(3) is complex balanced. Consequently, for each positive initial mass vector M ∈ R m >0 there exists a unique positive complex balanced equilibrium c ∞ ∈ R N >0 . Note that the system may possess (possibly infinitely) many boundary equilibria.
Let c(x, t) be a renormalised solution to (1)-(3) with initial data satisfying Q c 0 = M and E(c 0 |c ∞ ) < +∞. Note that any such renormalised solution satisfies the mass conservation laws Q c(t) = M , [26,28]. Assume that there exists a function H 1 : [0, +∞) → [0, +∞) with the property +∞ 0 Then, the renormalised solution c(x, t) converges exponentially to the positive complex balanced equilibrium c ∞ in the L 1 -norm with a rate, which can be explicitly computed in terms of the function H 1 , the domain Ω, the diffusion matrix D, the stoichiometric coefficients y ∈ C, the initial mass M , the complex balanced equilibrium c ∞ and the reaction rate constants k r .
The main progress of Theorem 1.2 is that the question of convergence to equilibrium for complex balanced reaction-diffusion systems with boundary equilibria is reduced to proving the finite dimensional inequality (16). Moreover, if the function H 1 (t) is explicitly computable (i.e. for some specific systems), then the rate of equilibration of the renormalised solution c(x, t) to (1)-(3) can also be computed explicitly. However, proving (16) for general systems with boundary equilibria remains a difficult problem since it requires suitable estimates on renormalised solutions, more precisely, on the behaviour of the L 1 -norm of renormalised solutions near the boundary ∂R N >0 , which is already a hard problem for ODE systems with boundary equilibria. Nevertheless, we will show in Subsection 3.2 how to apply Theorem 1.2 to specific systems.

Remark 2 (Towards a Global Attractor Conjecture for Reaction-Diffusion Systems).
It is worthwhile to remark on the key assumption +∞ 0 H 1 (s)ds = +∞. Note first that if H 1 (t) is defined as the fraction of the left-hand-side sum and the right-hand-side sum of (16), then H 1 (t) = 0 if and only if c(t) is a boundary equilibrium and otherwise bounded (if c(t) = c ∞ , then this H 1 (t) extends continuously to a positive constant). Hence, it is equivalent to consider in Theorem 1.2 the assumption +∞ t1 H 1 (s)ds = +∞ for a time t 1 arbitrarily large. Secondly, note that since (16) constitutes a finite dimensional inequality for the spatial averages c(t), one could conjecture to prove (16) by assuming the Global Attractor Conjecture for the corresponding ODE system (13), see [9] for a proof under review of the GAC for complex balanced ODE systems.
Indeed at a time t 1 > 0, consider the spatial averages c(t 1 ) as initial data u(t 1 ) = c(t 1 ) of (13). Then, the ODE Global Attractor Conjecture for (13) should imply (via a contradiction argument) the existence of H ODE for u to be the solution of (13), since the ODE system (13) shares the same complex balanced equilibria as the PDE system. Moreover, formal estimates seem to suggest that it is possible to establish bounds on provided that there is a good comparison between of the evolution of the ODE system u(t) and the evolution of the PDE system c(t) via its spatial averages c(t). Next at time t 2 , one restarts the ODE evolution (13) with a second set of initial data u(t 2 ) = c(t 2 ) and use that also this ODE system satisfies the GAC and yields another function H ODE 2 (t) on a time interval (t 2 , t 3 ) and so forth. Assuming that the evolution of c(t) converges sufficiently fast to these family of related ODE solutions, is seems possible to prove a statement like ODE GAC implies GAC for the PDE systems.
However, the problem of deriving good convergence estimates on the difference between the ODE system (13) and the evolution of the spatial averages c(t) seems to be (at least) as hard as understanding directly the evolution of c(t). First, the non-convexity of R(u) prevents any direct comparison between d dt u = R(u) = R(c) = d dt c(t). Moreover, the evolution of the difference c(t) − u(t) is not non-negative and doesn't seem to feature an entropy functional. Hence, it seems that in order to derive estimates on the difference c(t) − u(t), one is brought back to understanding the equilibration of c(t), which is the problem to solve at first.
Notations: Throughout this paper, we will use the following set of convenient notations: • Capital letters for square roots of corresponding normal letter, that is • For two vectors y = (y 1 , . . . , y N ) and z = (z 1 , . . . , z N ) in R N with z i = 0 for all i = 1, . . . , N , we write • For a function f : R → R and a vector y ∈ R N , we denote by For example, Organisation of the paper: In section 2, we present the proof of Theorem 1.1 and show its application to a reversible enzyme reaction. The proof of Theorem 1.2 and applications to networks with boundary equilibria will be presented in Section 3.  3), it was very recently proven in [28] that all renormalised solutions according to the following definition satisfy the weak entropy-entropy dissipation law (10) and the conservation laws (5). [27,28]).
Let Ω be a bounded domain in R n with Lipschitz boundary ∂Ω. Assume that the diffusion matrix D is positive definite, i.e. d i > 0 for all i = 1, . . . , N . Moreover, we assume that system (1) is complex balanced and thus possesses the entropy dissipation structure (7) and (8).
Moreover, any renormalised solution c(x, t) to (1)-(3) satisfies the weak entropy-entropy dissipation law for almost all t ≥ s ≥ 0, and the mass conservation laws, i.e. Q c(t) = Q c 0 for a.a. t > 0.

Preliminary estimates.
We present in this part some useful preliminary estimates which are needed for the sequel proofs.
The following Csiszár-Kullback-Pinsker type inequality shows that convergence to equilibrium in relative entropy implies convergence to equilibrium in L 1 -norm. For its proof, even in more general settings, we refer the reader to e.g. [3,11,12,15].  Proof. The proof follows from direct computations, hence we omit it here. Lemma 2.3 allows to prove the entropy-entropy dissipation inequality (9) by estimating E(c|c) and E(c|c ∞ ) separately. The first part E(c|c) can be easily controlled by D(c) thanks to the Logarithmic Sobolev inequality as in the following Lemma 2.4. For all measurable c : Ω → R N + with finite relative entropy E(c|c ∞ ) < +∞ holds D(c) ≥ λ 1 E(c|c) for λ 1 = min i=1,...,N {d i } C LSI where C LSI is the best constant in the Logarithmic Sobolev inequality.
Proof. By using the Logarithmic Sobolev inequality Thanks to the Lemmas 2.3 and 2.4, the remaining part of this section is dedicated to control the second part E(c|c ∞ ) of the relative entropy E(c|c ∞ ). Note that such a control has to quantify the system behaviour of the reacting concentrations c as well as the conservation laws Q c(t) = M = Q c 0 . Therefore, the control of E(c|c ∞ ) is much more challenging and technical.
We first show that E(c|c ∞ ) is bounded above by the right hand side of the finite dimensional inequality (11).
for an explicit constant K 1 > 0 depending on K and c ∞ (see (18)).
Proof. First, by using the elementary inequalities Next, we introduce the function with the extensions Φ(0) = lim z→0 Φ(z) = 1 and Φ(1) = lim z→1 Φ(z) = 2, and monotone increasing. By using now the bound c i ≤ K, we can estimate By using Lemmas 2.3, 2.4 and 2.5, where the latter establishes the right hand side of (11), our proof of Theorem 1.1 still requires i) to control the left hand side of the finite dimensional inequality (11) in terms of D(c), and ii) to prove (11). These will be done in Lemmas 2.6, 2.7 and Lemma 2.8 respectively.
As the first step, we observe that the entropy dissipation D(c) is a combination of the diffusion and reaction processes of the system (1), see (8). The reaction term seems hard to control due to the nonconvex nonlinearities (with arbitrary high order polynomials) and the very low regularity of renormalised solutions. In fact, we will not show that the entropy dissipation D(c) is bounded for renormalised solutions, but only that it constitutes an upper bound even while potentially unbounded. We will prove in the following lemma that, with the help of the diffusion terms, the reaction part is bounded below by "reactions of averaged concentrations". Herein, we recall the convention of square roots C i = √ c i and C i,∞ = √ c i,∞ and denote by C = (C 1 , . . . , C N ) and C = (C 1 , . . . , C N ) with C i = Ω C i dx.
Proof. By using the identity ∇ √ To prove (19), we use similar arguments to [22,15,26]. Fix a constant L > 0. The proof uses a domain decomposition corresponding to the deviation of C i around the averages C i , i.e. by denoting where S = {x ∈ Ω : |δ i (x)| ≤ L for all i = 1, . . . , N } and S c = Ω\S. We will see that on S the reaction part is crucial while the diffusion part is sufficient on S c .
Remark 3. The constant L in Lemma 2.6 can be chosen arbitrary. One certainly can choose L in order to optimise (i.e. maximise) the constant K 3 in (25). This may help to improve the rate of convergence. However, due to the multistage proof of the entropy-entropy dissipation inequality, the estimated rates are not optimal.
Now we are able to control the left hand side of (11) by D(c).
Lemma 2.7. For any measurable c : Ω → R N + satisfying E(c|c ∞ ) ≤ K and Q c = M , there exists an explicit constant K 2 > 0 (see (29)) such that Remark 4. Lemma 2.7 is a crucial step in proving Theorem 1.1. As mentioned in the introduction, we are here able to remove a technical assumption on cases when the L 1 -norm of the concentrations approaches the boundary ∂R N + , which was needed in [26,Theorem 1.4]. The key observation is the remainder estimate (27). Note that this idea was also used in [38, Lemma A.5] for energy-reactiondiffusion systems.

Proof. By denoting
with Note that Q(C i ) ≥ 0 and that Similarly to (21), we use Taylor expansion and ansatz (26) in which the constant H i is estimated as Hence, it follows from Lemma 2.6 and (28) and the Poincaré inequality that for any θ ∈ (0, 1) by choosing θ = min 1 2 ; C P (N |R| max i {H i }) −1 . The last step now is to prove the finite dimensional inequality (11). Let us recall that until this point, we have not used the fact that the system under consideration possesses no boundary equilibria. This fact turns out to be very useful when dealing with systems having boundary equilibria (see Section 3).
Lemma 2.8. Assuming that the chemical reaction network (S, C, R, K) is complex balanced and does not have any boundary equilibria. Then, for any c ∈ R N >0 satisfying E(c|c ∞ ) < +∞ and Q c = M , the inequality (11) holds for some constant H 1 > 0.
Remark 5. We remark here that while all the constants in previous lemmas can be explicitly estimated, the constant H 1 in (11) (as established in the this lemma) is in general not explicit since the proof utilises a contradiction argument. However, we believe that for any concrete system, where the conservation laws are explicitly known, H 1 can be computed explicitly via only elementary calculations (see Section 2.2). Estimating H 1 for general systems is a subtle issue since the structure of conservation laws, which is crucial for an explicit estimate, is unclear in general and remains thus an open problem.
Proof. Observe that the right hand side of (11) equals zero if and only if c = c ∞ . Therefore we first prove that the left hand side of (11) can only be zero when c ≡ c ∞ . Indeed, assuming that the left hand side of (11) is zero, then we have Thus, for any y ∈ C we have  (30)).
That means that c is a complex balanced equilibrium. Since the chemical reaction network has no other complex balanced equilibrium than c ∞ , we obtain the desired claim that c ≡ c ∞ . It is obvious that Ξ ≥ 0. Now assume by contradiction that Ξ = 0. By linearising both the nominator and denominator around c ∞ , and by setting σ = c − c ∞ and η = σ c∞ = σ1 c1,∞ , . . . , σN cN,∞ , we obtain Note that η is the same vector for all r = 1, . . . , |R| in the numerator. Note moreover that both numerator and denominator are of homogeneity two. We can thus rescale and normalise η w.l.o.g. and only consider η on the unit ball, that is |η| = 1. Moreover, Ξ = 0 if and only if the nominator is zero: which is only possible when η ∈ ker(W ), where we recall that W is the Wegscheider matrix W = [(y ′ r − y r ) r=1,...,|R| ] ⊤ ∈ R |R|×N . Recall that m = codim(W ) = dim(ker(W )) is the number of conservation laws. If m = 0 and the system (1)-(2) does not have a conservation law and equivalently ker(W ) = {0}, then it follows that η = 0, which is a contradiction to |η| = 1. If m > 0, then by using η ∈ ker(W ) and the fact that the rows of Q form a basis of ker(W ), it follows that η = Q ⊤ γ with some γ ∈ R m . Since which implies γ = 0 since Q has full rank. Thus η = 0 which again contradicts with |η| = 1.
In conclusion, we have proved that Ξ > 0, which implies the existence of a constant H 1 > 0 and hence completes the proof.
We can now begin the Proof of Theorem 1.1. From Lemmas 2.5, 2.7 and 2.8 we get which in combination with Lemma 2.4 leads to the desired estimate (9). Next, thanks to Theorem 2.1, any renormalised solution satisfies the conservation laws and the weak entropy-entropy dissipation law (10). Hence we can apply a variant version of Gronwall's inequality (see e.g. [18] or [59]) to get the exponential decay for almost all t > 0. This convergence in a combination with the Csiszár-Kullback-Pinsker type inequality in Lemma 2.2 leads to the claimed convergence to equilibrium (12).

2.2.
Applications to reversible enzyme reactions. Theorem 1.1 shows that any renormalised solution of complex balanced reaction-diffusion systems without boundary equilibria converges exponentially to equilibrium with a constant and a rate, which can be explicitly estimated up to the finite dimensional inequality (11). Proving (11) with an explicit constant H 1 seems to be a difficult task in full generality due to the non-convex nonlinear reaction terms and the non-explicit structure of conservation laws, i.e. due to the fact that we have no explicit structure of the constraints imposed by the matrix Q .
In this section, however, we will show that for a specific system, where the conservation laws are explicitly known, we can prove inequality (11) with an explicit constant H 1 by using elementary estimates. Hence we obtain convergence to equilibrium for (1) with explicit bounds for the convergence rates and constants in a highly relevant model of enzyme reactions.
For notational convenience, we use a change of variables and rewrite the finite dimensional inequality (11) in a form, which is easier to handle in the specific case at hand. By denoting for µ i ∈ [−1, +∞) and µ = (µ 1 , . . . , µ N ), inequality (11) rewrites as follow: where µ satisfies the following constraint inherited from the mass conservation laws and where we recall the convention c ∞ (µ 2 + 2µ) = (c i,∞ (µ 2 i + 2µ i )) i=1,...,N . We apply our approach to a reversible variant of the famous Michaelis-Menten enzyme reaction For the sake of clarity, we shall assume k 1 = k 2 = k 3 = k 4 = 1, but we emphasise that the subsequent analysis can be equally carried out for general k i > 0, i = 1, . . . , 4 without additional technical difficulties. The corresponding mass action reaction-diffusion system reads as x ∈ Ω, t > 0, (34) with homogeneous Neumann boundary conditions ∇c i · ν = 0 on ∂Ω and non-negative initial data c i (x, 0) = c i,0 (x) ≥ 0, i = 1, . . . , 4, in which Ω is a bounded domain with sufficiently smooth boundary (e.g. ∂Ω ∈ C 2+ǫ with ǫ > 0) and normalised volume |Ω| = 1. The large time behaviour of various reaction-diffusion models of reversible enzyme kinetics has also been recently studied in e.g. [16,Section 8] or [26]. It is easy to check that there are two linear independent mass conservation laws for (34) and that the matrix Q can be chosen as It is straightforward to check that this equilibrium is a complex balanced equilibrium (and even a detailed balanced equilibrium) and that system (34) possesses no boundary equilibria. The existence of global renormalised solution to (34) follows immediately from Theorem 2.1. Moreover, since the nonlinearities in (34) are quadratic, it is well-known (see e.g. [48,49]) that (34) has a global weak solution. Moreover, thanks to the special structure of (34), we show in the following that these weak solutions are in fact strong solutions and grow at most polynomially in time.
Remark 6. Note that the L ∞ -bounds of Proposition 2.9 are sufficient to apply standard parabolic bootstrap arguments and show that c(x, t) is indeed a classical solution (or even smooth if ∂Ω is smooth) and thus unique.
Proof of Proposition 2.9. The proof relies on duality estimates and comparison principle arguments for scalar parabolic equations, which exploit the special structure of (34). In this proof we always denote by C T a general constant depending polynomially on T > 0. First, it follows from (34) that By a classical duality estimate (see e.g. [50]) and by denoting L 2 (Q T ) = L 2 (0, T ; L 2 (Ω)), we have Moreover, (34) is quasi-positive in the sense of e.g. [49] and thus preserves non-negativity of weak solutions c 1 , . . . , c 4 from non-negative initial data. This, implies Next, we by recalling [6, Lemma 3.3], there exists a constant C T , which depend polynomially on T and quantifies the smoothing effect of the heat operator in the following sense: Given f ∈ L p (Q T ) and let v be the solution to v t − d∆v = f subject to homogeneous Neumann boundary condition. Then, • if p < (N + 2)/2 then v L s−ǫ (QT ) ≤ C T for any ǫ > 0 with s = (N +2)p N +2−2p , • if p ≥ (N + 2)/2 then v L r (QT ) ≤ C T for all 1 ≤ r < +∞.
Therefore it follows from (37) and c 3 ∈ L 2 (Q T ) in particular that and On the other hand, by another duality estimate (see e.g. [52,Lemma 33.3]), it follows from that the regularity and the polynomial dependence of C T on T are transferred from c 1 to c 3 , which implies that By repeating this procedure, we obtain after finitely many steps that c 3 L q (QT ) ≤ C T with q ≥ N +2 2 . Then, (37) implies for all r ∈ [1, ∞), which yields in return c 3 L r (QT ) ≤ C T for all r ∈ [1, ∞). Hence, after one application of the classical smoothing effect of heat operator, the proof of the Proposition is completed. i=1 Ω c i,0 log c i,0 dx < +∞, converges in L 1 exponentially to the unique positive equilibrium c ∞ as defined in (36): where the constant C and the rate λ can be explicitly estimated in terms of Ω, the equilibrium c ∞ and initial masses M 13 and M 234 . Moreover, if the initial data c 0 belongs to L ∞ (Ω) 4 , then (34) has a unique global classical solution, which converges exponentially to c ∞ in any L p -norm for 1 ≤ p < ∞, i.e.
for all t ≥ 0, (38) with explicit constant C and rate λ ′ .
Proof. Since the system satisfies the complex balanced condition and possesses no boundary equilibria, Theorem 1.1 implies immediately that any renormalised solution converges exponentially to the equilibrium defined in (36). It remains to bound of convergence rate explicitly. Thanks to Theorem 1.1 and (32) that means to compute explicitly a constant H enzyme 1 > 0 in the finite dimensional inequality satisfying the following constraints, which are equivalent to the mass conservation laws (35): Note that (39) is the speficic form of inequality (32) in case of the reversible enzyme reaction (E). Let G denote the left hand side of (39). First, the elementary inequality From (40a) and by observing that (µ 1 + 2), (µ 3 + 2) ≥ 1, it follows directly that µ 1 and µ 3 must have different signs, which leads to the following two cases: i) Consider µ 1 ≥ 0 and µ 3 ≤ 0: First, we have From (40b) and µ 3 ≤ 0, we infer that at least either µ 2 ≥ 0 or µ 4 ≥ 0, which leads to two subcases: ia) Suppose µ 2 ≥ 0. Then, since µ 1 µ 3 ≤ 0. Similarly, It therefore follows from (41), (42) and (43) that where we have used Young's inequality and the factor 1 18 is not sharp but chosen in order to obtain the same lower bound (44) as in the second case below. ib) Suppose µ 4 ≥ 0. In this case the term [(1 + µ 3 ) − (1 + µ 1 )(1 + µ 4 )] 2 can be estimated analog to case ia). ii) Consider µ 1 ≤ 0 and µ 3 ≥ 0: Since µ 3 ≥ 0, we obtain from (40b) that at least either µ 2 ≤ 0 or µ 4 ≤ 0. We can then use the same arguments to (42) and (43) to imply that if µ 2 ≤ 0 then and if µ 4 ≤ 0 then However, because µ 1 ≤ 0, the inequality (41) is not valid anymore. In order to bypass it, we need to consider two subcases concerning the closeness of µ 1 to −1.
3. Proof of Theorem 1.2 and applications to systems with boundary equilibria 3.1. Proof of Theorem 1.2.
Proof of Theorem 1.2. We already mentioned in the proof of Lemma 2.8, that the validity of the Lemmas 2.3, 2.4, 2.5 and 2.7 is independent of the presence or absence of boundary equilibria. We recall here the key estimates of the Lemmas for the sake of readability: The additivity of the relative entropy allows to control the term E(c(t)|c(t)) via the Logarithmic Sobolev inequality in terms of the entropy dissipation, i.e. E(c(t)|c ∞ ) = E(c(t)|c(t)) + E(c(t)|c ∞ ) and λ 1 E(c(t)|c(t)) ≤ D(c(t)).

The second term E(c(t)|c ∞ ) satisfies the upper bound
while the entropy dissipation obeys the lower bound These two estimates are connected by assumption (16) and we obtain all together Note that +∞ 0 λ(s)ds = +∞ since +∞ 0 H 1 (s)ds = +∞ and λ 1 > 0. Moreover, it follows from the weak entropy-entropy dissipation law (10) and Gronwall's inequality that Thus the trajectory c(t) converges to c ∞ in relative entropy and, consequently, in L 1 -norm due to the Csiszár-Kullback-Pinsker type inequality in Lemma 2.2. Therefore, after some finite time T > 0, the solution trajectory will always stays outside of any small enough neighbourhood of all boundary equilibria. It then follows from [15,Remark 3.6] that the solution converges exponentially to the positive complex balanced equilibrium.

3.2.
Application to a specific system possessing boundary equilibria.
In order to show convergence to equilibrium for renormalised solutions c(x, t) of complex balanced reaction-diffusion systems with boundary equilibria, we have to verify (16) as stated in Theorem 1.2.
Similarly to Subsection 2.2, it will be convenient to change variables in the finite dimensional inequality (16). By setting c i (t) = c i,∞ (1 + µ i (t)) 2 , for i = 1, . . . , N, inequality (16) becomes where µ(t) = (µ 1 (t), . . . , µ N (t)) and the function H 1 (t) is required to satisfy +∞ 0 H 1 (t)dt = +∞. Proving (54) for general complex balanced systems would yield a proof of the Global Attractor Conjecture (GAC), which is a very interesting, yet challenging open problem. Our aim in this section is to study a typical class of complex balanced systems with boundary equilibria, in which proving (54) for renormalised solutions is a possible approach to answer the GAC in the associated PDE setting. More precisely, we consider here the reaction-diffusion systems modelling the following reaction network with arbitrary α ≥ 1 and k 1 , k 2 , k 3 > 0. The special case α = 1 was investigated in [15]. Here, we study the entire range α ≥ 1 in order to show the robustness of our arguments. The above network (C) is considered in a bounded domain Ω ⊂ R n with smooth boundary ∂Ω (e.g. C 2+ǫ for any ǫ > 0) and normalised volume, i.e. |Ω| = 1. The corresponding mass action reactiondiffusion system reads as subject to homogeneous Neumann boundary conditions ∇c i ·ν = 0 and non-negative initial data c i (x, 0) = c i,0 (x). This system has one conservation of mass, namely (α + 1) c 1 (t) + c 2 (t) + c 3 (t) = (α + 1) c 1,0 + c 2,0 + c 3,0 =: M, for all t > 0. (56) For fixed M > 0, system (55) features the boundary equilibrium c * = (c * 1 , c * 2 , c * 3 ) = (0, 0, M ) and the unique positive complex balanced equilibrium c ∞ = (c 1,∞ , c 2,∞ , c 3,∞ ), where c 2,∞ is the unique positive solution to Due to the presence of the boundary equilibrium, we will have to apply Theorem 1.2 in order to show convergence to equilibrium for renormalised solutions to (55). More precisely, we need to prove the modified finite dimensional inequality (16) or equivalently (54) along solution trajectories of (55). The existence of global renormalised solutions to the complex balanced system (55) follows readily from Theorem 2.1.
However, to prove inequality (54) along renormalised solutions, we need additional information about these solutions, which we are only able to show for specific renormalised solutions constructed via a typical approximation scheme as already used in [27]. The following lemma shows that if such a renormalised solution to (55) should converge to the boundary ∂R 3 >0 , then not faster than with a specific algebraic convergence rate in terms of the parameter α ≥ 1.
Moreover, assume 1/c α 2,0 L ∞ (Ω) < +∞. Then, any renormalised solution, which is constructed via the below approximative scheme (59), satisfies We also remark that the assumed L p , 1 < p ≤ 2 initial data, which we need for technical reasons in order to apply duality estimates, are slightly more restrictive than the usual L log L initial data assumption for renormalised solutions (see Theorem 2.1). Note that for α sufficiently larger than one, the existence of global weak solutions to system (55) is unclear even with L 2 initial data and that renormalised solutions are the only known global solutions in order to study the large time behaviour.
subject to initial data c ε 2,0 ≥ Cε for a constant C and for all ε > 0 imply the existence of a positive time τ > 0 (possibly depending in ε), such that c ε 2 (x, t) ≥ Cε 2 for a.a. x ∈ Ω and 0 ≤ t ≤ τ . Thus, we can test Thus, the weak comparison principle implies again for a.a. x ∈ Ω and t > 0 which implies that testing with − α (c ε 2 ) α+1 is justified for a.a. t ≥ 0 and the lower bound (61) holds indeed globally in time and independently from ε.
Hence c ε 2 (t) ≥ h(t) with h(t) as defined in (58). Finally, the estimate (58) follows from (61) and the fact that c ε 2 (t) → c 2 (t) in L 1 (Ω). Remark 8. If the diffusion coefficients d 1 , d 2 , d 3 are close to each other, for instance, in the sense that δ = max{d 1 , d 2 , d 3 } − min{d 1 , d 2 , d 3 } is sufficiently small, then any renormalised solution to (55) is in fact a strong solution, see [6]. In these cases the arguments in Proposition 3.1 are justified by classical maximum principle arguments as done in [15]. The benefit of Proposition 3.1 is to prove estimate (58) for suitable renormalised solutions without any assumption on the diffusion coefficients. This is due to the fact that (58) involves only the L 1 -norm of c 2 , which is preserved when passing to the limit in approximating renormalised solutions.
We are now ready to prove the main result of this section. Let Ω ⊂ R n be a bounded domain with smooth boundary ∂Ω (e.g. C 2+ǫ with ǫ > 0). Assume for system (55) that α ≥ 1 and k 1 , k 2 , k 3 > 0.
Then, for any fixed positive initial mass M > 0 and non-negative initial data (c 1,0 , c 2,0 , c 3,0 ) ∈ L p (Ω) 3 for some 1 < p ≤ 2 having initial mass M as defined in (56) and satisfying any global renormalised solution (c 1 , c 2 , c 3 ) as constructed in Proposition 3.1 converges exponentially in L 1 to the complex balanced equilibrium (c 1,∞ , c 2,∞ , c 3,∞ ) as defined in (57), i.e. for almost all t > 0 where C and λ are constants depending explicitly on the domain Ω, the constants α, k 1 , k 2 , k 3 , the initial mass M and 1/c α 2,0 L ∞ (Ω) . Remark 9. The exponential convergence to equilibrium in Theorem 3.2 applies to any renormalised solution constructed via the approximation scheme of Proposition 3.1. Note that it is unknown if any renormalised solution according to the definition in Theorem 2.1 can be approximated via (59). Thus, the convergence of any renormalised solution is open for future investigation.

Aknowledgements.
This work is partially supported by International Research Training Group IGDK 1754 and NAWI Graz.