Convergence to the complex balanced equilibrium for some chemical reaction-diffusion systems with boundary equilibria

TL;DR

This paper investigates the rate at which certain chemical reaction-diffusion systems with boundary equilibria converge to their complex balanced equilibrium, providing new insights into their stability and convergence behavior.

Contribution

It introduces new convergence results for reaction-diffusion systems with boundary equilibria, including a three-species system and a general two-species network.

Findings

Convergence to equilibrium is established for specific reaction-diffusion systems.

The rate of convergence is characterized for systems with boundary equilibria.

Results apply to systems with quadratic nonlinearities in the positive orthant.

Abstract

In this paper we study the rate of convergence to the complex balanced equilibrium for some chemical reaction-diffusion systems with boundary equilibria. We first analyze a three-species system with boundary equilibria in some stoichiometric classes, and whose right hand side is bounded above by a quadratic nonlinearity in the positive orthant. We prove similar results on the convergence to the positive equilibrium for a fairly general two-species reversible reaction-diffusion network with boundary equilibria.

Equations

$${\mathbf{c}}_t-\mathcal{D}\Delta \mathbf{c}=R(\mathbf{c})\text{mbox}in\Omega \times (0,T),\mathbf{c}(\cdot ,0)={\mathbf{c}}_0\text{mbox}in\Omega ,$$

$$\begin{array}[]{l}\displaystyle a_{t}-d_{a}\Delta a=-ka,\\ \displaystyle b_{t}-d_{b}\Delta b=ka\mbox{ in }\Omega\times(0,T),\\ \displaystyle a(\cdot,0)=a_{0},\ b(\cdot,0)=b_{0}\mbox{ in }\Omega,\end{array}$$

$$\begin{array}[]{l}\displaystyle a_{t}-d_{a}\Delta a=-k_{1}a+k_{2}b,\\ \displaystyle b_{t}-d_{b}\Delta b=k_{1}a-k_{2}b\mbox{ in }\Omega\times(0,T),\\ \displaystyle a(\cdot,0)=a_{0},\ b(\cdot,0)=b_{0}\mbox{ in }\Omega,\end{array}$$

$$\begin{array}[]{l}\displaystyle a_{t}-d_{a}\Delta a=-kab,\\ \displaystyle b_{t}-d_{b}\Delta b=-kab,\\ \displaystyle c_{t}-d_{c}\Delta c=kab\mbox{ in }\Omega\times(0,T)\\ \displaystyle a(\cdot,0)=a_{0},\ b(\cdot,0)=b_{0},\ c(\cdot,0)=c_{0}\mbox{ in }\Omega.\end{array}$$

$$k{\mathbf{c}}^y=k\prod _{i=1}^n{c}_i^{y_i},$$

$$R(\mathbf{c}):=\sum _{y\to y^{\prime }}{k}_{y\to y^{\prime }}{\mathbf{c}}^y(y^{\prime }-y),$$

$$\begin{array}[]{l}\displaystyle a_{t}-d_{a}\Delta a=-k_{1}ab+k_{2}c,\\ \displaystyle b_{t}-d_{b}\Delta b=-k_{1}ab+k_{2}c,\\ \displaystyle c_{t}-d_{c}\Delta c=k_{1}ab-k_{2}c\mbox{ in }\Omega\times(0,T)\\ \displaystyle a(\cdot,0)=a_{0},\ b(\cdot,0)=b_{0},\ c(\cdot,0)=c_{0}\mbox{ in }\Omega.\end{array}$$

$$\sum _{\bar{y}\to y}{k}_{\bar{y}\to y}{c}_0^{\bar{y}}=\sum _{y\to \bar{y}}{k}_{y\to \bar{y}}{c}_0^y$$

$$E(t):=\sum_{i=1}^{n}\int_{\Omega}\bigg{[}c_{i}(x,t)\log\frac{c_{i}(x,t)}{c_{i,\infty}}-c_{i}(x,t)+c_{i,\infty}\bigg{]}dx,$$

$$D(t):=\sum_{i=1}^{n}d_{i}\int_{\Omega}\frac{|\nabla c_{i}(x,t)|^{2}}{c_{i}(x,t)}dx+\sum_{r=1}^{\rho}k_{r}c_{\infty}^{y_{r}}\int_{\Omega}\Phi\bigg{(}\frac{c^{y_{r}}}{c_{\infty}^{y_{r}}};\frac{c^{y^{\prime}_{r}}}{c_{\infty}^{y^{\prime}_{r}}}\bigg{)}dx,$$

$$D(t)\ge \alpha E(t)\text{mbox}forallt\ge 0.$$

$$\{\begin{array}{ll}{a}_t-d_a\Delta a=-k_{1}a{b}^n+k_{2}bc& \\ b_t-d_b\Delta b=-k_{1}a{b}^n+k_{2}bc& \text{mbox}in\Omega \times (0,\infty )\\ c_t-d_c\Delta c=k_{1}a{b}^n-k_{2}bc& \\ \nabla a\cdot \nu =\nabla b\cdot \nu =\nabla c\cdot \nu =0& \text{mbox}on\partial \Omega \times (0,\infty )\\ a(\cdot ,0)=a_0,b(\cdot ,0)=b_0,c(\cdot ,0)=c_0& \text{mbox}in\Omega ,\end{array}$$

$$\{\begin{array}{ll}{a}_t-d_a{a}_{xx}=-a{b}^2+bc& \\ b_t-d_b{b}_{xx}=-a{b}^2+bc& \text{mbox}in(0,1)\times (0,\infty )\\ c_t-d_c{c}_{xx}=a{b}^2-bc& \\ a_x=b_x=c_x=0& \text{mbox}on\{0,1\}\times (0,\infty )\\ a(\cdot ,0)=a_0,b(\cdot ,0)=b_0,c(\cdot ,0)=c_0& \text{mbox}in(0,1).\end{array}$$

$$m_{1}A+n_{1}B\rightleftharpoons m_{2}A+n_{2}B.$$

$$\{\begin{array}{ll}{a}_t-d_a\Delta a=\bar{m}(a^{m_2}{b}^{n_2}-a^{m_1}{b}^{n_1})& \text{mbox}in\Omega \times (0,\infty )\\ b_t-d_b\Delta b=\bar{n}(a^{m_1}{b}^{n_1}-a^{m_2}{b}^{n_2})& \text{mbox}in\Omega \times (0,\infty )\\ \nabla a\cdot \nu =\nabla b\cdot \nu =0& \text{mbox}on\partial \Omega \times (0,\infty )\\ a(\cdot ,0)=a_0,b(\cdot ,0)=b_0& \text{mbox}in\Omega .\end{array}$$

$$\begin{array}{r}E(a,b,c)={\int }_0^{1}a(\ln a-1)dx+{\int }_0^{1}b(\ln b-1)dx+{\int }_0^{1}c(\ln c-1)dx\end{array}$$

$$\begin{array}{rl}& D(a,b,c)=4d_a{\int }_0^1|{\partial }_x\sqrt{a}{|}^{2}dx+4d_b{\int }_0^1|{\partial }_x\sqrt{b}{|}^{2}dx+4d_c{\int }_0^1|{\partial }_x\sqrt{c}{|}^{2}dx\\ & +{\int }_0^1(a{b}^2-bc)\ln (\frac{a{b}^2}{bc})dx.\end{array}$$

$$\begin{array}{rl}& {\int }_0^{1}a(x,t)dx+{\int }_0^{1}c(x,t)dx={\int }_0^1{a}_0(x)dx+{\int }_0^1{c}_0(x)dx=:M_1,\\ & {\int }_0^{1}b(x,t)dx+{\int }_0^{1}c(x,t)dx={\int }_0^1{b}_0(x)dx+{\int }_0^1{c}_0(x)dx=:M_2,\end{array}$$

$$\Vert a{\Vert }_{L^2([0,1]\times [\tau ,\tau +1])},\Vert b{\Vert }_{L^2([0,1]\times [\tau ,\tau +1])},\Vert c{\Vert }_{L^2([0,1]\times [\tau ,\tau +1])}\le C$$

$$\begin{split}\bigg{|}\sqrt{a(x,t)}-\int_{0}^{1}\sqrt{a(y,t)}dy\bigg{|}\leq\int_{0}^{1}\big{|}\partial_{y}\sqrt{a(y,t)}\big{|}dy,\end{split}$$

$$\begin{split}a(x,t)&\leq\bigg{(}\int_{0}^{1}\sqrt{a(y,t)}dy+\int_{0}^{1}\big{|}\partial_{y}\sqrt{a(y,t)}\big{|}dy\bigg{)}^{2}\\ &\leq 2\bigg{(}\int_{0}^{1}\sqrt{a(y,t)}dy\bigg{)}^{2}+2\bigg{(}\int_{0}^{1}\big{|}\partial_{y}\sqrt{a(y,t)}\big{|}dy\bigg{)}^{2}\\ &\leq 2\int_{0}^{1}{a(y,t)}dy+2\int_{0}^{1}\big{|}\partial_{y}\sqrt{a(y,t)}\big{|}^{2}dy.\end{split}$$

$$E(a(t),b(t),c(t))+{\int }_0^{t}D(a(s),b(s),c(s))ds=E(a_0,b_0,c_0),$$

$$\begin{array}{rl}& E(a(t),b(t),c(t))+4d_a{\int }_0^t{\int }_0^1|{\partial }_x\sqrt{a}{|}^{2}dxdt+4d_b{\int }_0^t{\int }_0^1|{\partial }_x\sqrt{b}{|}^{2}dxdt+\\ & 4d_c{\int }_0^t{\int }_0^1|{\partial }_x\sqrt{c}{|}^{2}dxdt\le E(a_0,b_0,c_0).\end{array}$$

$$\begin{array}{rl}& 4d_a{\int }_0^t{\int }_0^1|{\partial }_x\sqrt{a}{|}^{2}dxdt+4d_b{\int }_0^t{\int }_0^1|{\partial }_x\sqrt{b}{|}^{2}dxdt+4d_c{\int }_0^t{\int }_0^1|{\partial }_x\sqrt{c}{|}^{2}dxdt\\ & \le E(a_0,b_0,c_0)+3,\end{array}$$

$$\begin{array}{rl}& \Vert {\partial }_x\sqrt{a}{\Vert }_{L^2([0,1]\times [0,t])}^2+\Vert {\partial }_x\sqrt{b}{\Vert }_{L^2([0,1]\times [0,t])}^2+\Vert {\partial }_x\sqrt{c}{\Vert }_{L^2([0,1]\times [0,t])}^2\\ & \le \frac{E(a_0,b_0,c_0)+3}{4d}:=\text{nolinebreak}[4]C_1\end{array}$$

$$\begin{split}&\int_{\tau}^{\tau+1}\int_{0}^{1}{a^{2}}dxdt\leq\int_{\tau}^{\tau+1}{\|a(t)\|_{L^{\infty}[0,1]}}\bigg{(}\int_{0}^{1}{a(x,t)}dx\bigg{)}dt\\ &\leq 2M_{1}\int_{\tau}^{\tau+1}\bigg{(}\int_{0}^{1}{a(x,t)}dx+\int_{0}^{1}|\partial_{x}\sqrt{a(x,t)}|^{2}dx\bigg{)}dt\\ &\leq 2M_{1}\int_{\tau}^{\tau+1}\bigg{(}M_{1}+\int_{0}^{1}|\partial_{x}\sqrt{a(x,t)}|^{2}dx\bigg{)}dt\leq 2{M_{1}}^{2}+2M_{1}C_{1}=:C\end{split}$$

$$u_{i,t}-{\kappa }_i\Delta u_i=f_i(x,t,u),i=1,...,m,$$

$$\mbox{ There exist }K_{1},\ K_{2}>0,\ 1\leq p<\infty,\ 1\leq r<1+\bigg{[}1-\frac{d}{p(d+2)}\bigg{]}\frac{2p}{d+2}\mbox{ such that }$$

$$\text{mbox}foreach1\le j\le m\text{mbox}thereexist{\alpha }_{j,k},1\le j\le k\text{mbox}with{\alpha }_{j,j}=1\text{mbox}suchthat$$

$$\sum _{k=1}^j{\alpha }_{j,k}{f}_j(x,t,v)\le K_1|v{|}^r+K_2\text{mbox}forallv\text{mbox}inthepositiveorthantof{\mathbb{R}}^m.$$

Full text

1 Introduction

The dynamical behavior of spatially homogeneous mass-action reaction systems has been the focus of much research over the last fifty years. These ODE systems are usually high-dimensional, non-linear, and depend on a large number of parameters, which makes them generally difficult to study. However, a fertile theory started fifty years ago with work of Horn, Jackson and Feinberg [20, 21, 15] has been successful in addressing questions of existence and stability of positive equilibria, and persistence (nonextinction) of variables. Their work shows that, surprisingly, the large class of complex balanced mass-action systems have unique positive equilibria and admit a global Lyapunov function, which makes them locally asymptotically stable independently of reaction rate constant values. This robustness is relevant in applications, where exact values of system parameters are typically unknown. Moreover, Horn conjectured that the unique equilibria are in fact globally asymptotically stable [21], a question known as the Global Attractor Conjecture. The conjecture stayed open until recent years, when new work fueled in part by advances in systems biology led to a series of partial results. It was shown that trajectories of complex balanced systems either converge to the positive equilibrium or go to boundary equilibria [28, 29], establishing that persistence implies global stability. A series of subsequent papers showed persistence for complex balanced systems in two variables and other classes of systems, and proved the Global Attractor Conjecture in two and three variables [1, 7, 2, 3, 25, 8, 19, 19]. This work led to a very recent proof of the Global Attractor Conjecture in full generality by Craciun [9].

Much less is known about the corresponding reaction-diffusion models, although a number of recent papers have focused on extending the results above in the PDE setting. A promising venue for relating the PDE and ODE models is by way of space discretization (the method of lines). As proof of concept, the network $A+B\rightleftharpoons C$ was considered in [24] where it was shown that solutions of the discretized system converge to the solution of the PDE system as the space discretization grows finer. Solutions of the reaction-diffusion system $A+B\rightleftharpoons C$ have been shown to approach a positive spatially homogeneous distribution [27] via semigroup theory. Newer work uses entropy techniques to prove global asymptotic stability for other systems, including dimerization networks $2A\rightleftharpoons B$ [10] and monomolecular networks [16].

For general complex balanced systems it was shown that under the assumption of equal diffusion constants, their $\omega$-limit set consists of constant functions corresponding to equilibria of the space-homogeneous ODE system and that, moreover, the unique positive equilibrium is asymptotically stable. This is the analogous of the “persistence implies global stability” result from the ODE setting, albeit in the case of equal diffusion constants. Recent results by Desvillettes, Fellner and collaborators [12, 14] removed the requirement of equal diffusion constants, and showed that in the absence of boundary equilibria, the positive equilibrium of general complex balanced reaction-diffusion systems attracts all solutions with positive initial data. These papers also considered special cases of networks with boundary equilibria, where a more detailed analysis showed that positive solutions remain globally asymptotically stable. However, the general case of systems with boundary equilibria remains open, and the analysis of such systems is on a case-by-case basis.

Our paper studies two cases of complex balanced reaction networks with boundary equilibria, and shows that under mild boundedness conditions on the initial data, solutions converge asymptotically to the unique positive equilibria. Namely, we consider the three-species system $A+nB{\rightleftharpoons}B+C$ (Theorem 1.1), and the two-species system $m_{1}A+n_{1}B{\rightleftharpoons}m_{2}A+n_{2}B$ (Theorem 1.2).

In the remainder of this introductory section we set up terminology and notation, we discuss some of the techniques used here and in previous work, and we state our main theorems. Sections 2 and 3 contain the proofs of the results for the three-species and two-species systems. We conclude with a few remarks and open problems (Section 4) and with an appendix collecting a few technical results needed in the paper.

1.1 Terminology and previous results

Let us consider $0<T\leq\infty$ and the semilinear parabolic system

[TABLE]

where $\mathbf{c}:\Omega\times[0,T)\rightarrow\mathbb{R}^{n}$ is the vector of concentrations at spatial position $x\in\Omega$ (an open subset of $\mathbb{R}^{d}$) and time $t\in[0,\infty)$, $\mathcal{D}$ is a positive definite, diagonal $n\times n$ matrix, and $R:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is a vector field whose components are polynomials (determined by the chemical reactions under consideration). This system can be linear and “trivial” (at least in the sense that “enough” of its equations decouple), such as

[TABLE]

(which corresponds to the reaction $A\displaystyle\rightarrow B$ with reaction rate $k>0$), linear and nontrivial (weakly coupled) such as

[TABLE]

(which corresponds to $A\overset{k_{1}}{\underset{k_{2}}{\rightleftharpoons}}B$). However, as soon as a reaction includes two or more reactants, the system becomes nonlinear in the zero order terms (semilinear). For example, the single reaction $A+B\overset{k}{\rightarrow}C$ yields

[TABLE]

We use this last system to illustrate some terminology and notation. Here $A,\ B,$ and $C$ are the three species of the network, and $A+B$ and $C$ are its complexes. In general, complexes are formal linear combinations of species with non-negative integer coefficients, and sit on both sides of a reaction arrow. It is useful to think of complexes as vectors in a natural way, for example $A+B$ corresponds to $y=(1,1,0)$, and $C$ to $y^{\prime}=(0,0,1)$. The concentrations of $A,\ B,\ C$ are non-negative functions of time and space and are collected in the concentration vector $\mathbf{c}=(a,b,c)$. The reaction rate of a reaction is given by mass-action, and is proportional to the concentration of each reactant species. This way, the reaction $A+B\overset{k}{\rightarrow}C$ has rate $kab$. The reaction rate constant $k$ is a reaction-specific positive number. In general, the rate of a the reaction $y\overset{k}{\to}y^{\prime}$ is given by

[TABLE]

where $n$ is the number of species, and complexes $y$ and $y^{\prime}$ are viewed as vectors, as illustrated above. The reaction rate $kab$ enters with negative sign in the equations for $a_{t}$ and $b_{t}$ ($A$ and $B$ are being consumed in the reaction), and with positive sign in the equation for $c_{t}$ ($C$ is being produced). The aggregate contribution of all reaction rates are collected in the vector $R(\mathbf{c})=(-kab,-kab,kab)$. In general, this is given by

[TABLE]

where $k_{y\to y^{\prime}}$ is the rate constant of $y\to y^{\prime}$ and the summation is over all reactions $y\to y^{\prime}$ in the network. Finally, $\mathcal{D}=\mathrm{diag}\{d_{a},d_{b},d_{c}\}\in M_{3\times 3}(\mathbb{R})$ denotes the diagonal matrix of diffusion constants.

In the previous example the first two equations have the benefit of being decoupled, but that feature is lost as soon as we allow for reversibility; indeed, corresponding to $A+B\overset{k_{1}}{\underset{k_{2}}{\rightleftharpoons}}C$ we have

[TABLE]

When it comes to basic questions on the existence, uniqueness, smoothness and non-negativity of solutions (if the initial data components are nonnegative), for linear systems the answers are provided in the (by now, classical) literature (see, e.g., [26]). However, complexity adds quickly as more reactions and/or more reactants enter the system. There is no general result in the literature that guarantees long time existence, uniqueness and non-negativity of solutions, let alone smoothness and other, more delicate properties such as comparison principles and convergence to equilibrium. In this paper we discuss chemical reaction diffusion systems which have a specific structure relative to a positive equilibrium, i.e. a steady state solution with all positive components.

In general, we say that an equilibrium point $c_{0}$ is a complex balanced equilibrium if for all complexes $\bar{y}$ we have

[TABLE]

i.e., the total chemical flux that exits the complex $\bar{y}$ equals the total chemical flux that enters the complex $\bar{y}$ (for any choice of $\bar{y}$). A chemical system is called a complex balanced system if it admits a positive complex balanced equilibrium. Due to the particular polynomial nature of $R(\mathbf{c})$ in this case, the steady states are all constant vectors (i.e. independent of location, as well as of time). If at least one component of a steady state is null, then the corresponding state is said to be a boundary equilibrium. All systems arising from complex balanced CRDSs admit a “canonical” Lyapunov functional of the relative Boltzmann entropy type. Its general form (again, see, e.g., [14]), this logarithmic free relative energy functional reads

[TABLE]

where $\mathbf{c}_{\infty}:=(c_{1,\infty},...,c_{n,\infty})$ is the constant vector denoting the positive complex balanced equilibrium. The entropy dissipation functional is computed by differentiating $E$ along trajectories; that is, once all the time derivatives of concentrations are replaced by their equation specific expressions and the Neumann BC are used to integrate by parts wherever the Laplacian appears, one gets

[TABLE]

where $\rho$ is the number of reactions and $\Phi(x,y):=x\log(x/y)-x+y$. Of course, one gets exponential decay to zero for $E$ if one can prove that there exists a positive constant $\alpha$ such that

[TABLE]

Naturally, $E(t)$ should not only be identically zero when $\mathbf{c}(t)=\mathbf{c}_{\infty}$, but it should also be bounded below by some increasing function of the distance (from some norm) between $\mathbf{c}(t)$ and $\mathbf{c}_{\infty}$.

For complex balanced systems, in the spatially isotropic case ($\mathcal{D}=0$, so the PDE’s are reduced to ODE’s) there is some recent work by Craciun [9], which answers in the affirmative a long standing conjecture on the convergence to the positive equilibrium in each stoichiometric class, called the Global Attractor Conjecture. This conjecture states that regardless of the existence of boundary equilibria, trajectories starting in the positive orthant converge to the unique positive equilibrium in the corresponding stoichiometric class. In the PDE case, the most general result concerns the case where there are no boundary equilibria. Very recently, Desvillettes, Fellner and Tang [14] showed that, contingent on the existence of suitable solutions (essentially, solutions that may not be classical but they are renormalized and do satisfy a weak entropy entropy-dissipation law), one obtains exponentially fast convergence to the equilibrium which lies in the same stoichiometric class as the initial data, which is merely assumed nonnegative and integrable over some bounded, $C^{2}$ domain in $\mathbb{R}^{d}$. This is also a remarkably general result in the sense that the initial concentrations are only assumed to lie in $L^{1}(\Omega)$. This improvement (over the previous works, where $L^{\infty}$-bounds were imposed on the initial data) is achieved via the use of the Log-Sobolev inequality (see, e.g., [14]) in order to establish the entropy-entropy dissipation inequality (EEDI) (3). In all the previous works, the EEDI follows from the standard zero-average Poincaré inequality applied to the square roots of the concentration functions, combined with their uniform $L^{\infty}$-bounds (in space-time); these bounds need to be proved a priori. This uniform $L^{\infty}$-bound is key to the proofs of convergence to equilibrium in most of the works on this topic, (in fact, to our knowledge, the only exception comes when there are no boundary equilibria [14]) and the constant $\alpha$ from the EEDI (3) tends to vanish as the $L^{\infty}$-bound on the solution blows up. We note that (2) is one of the two systems studied in [10], and the authors use the uniform $L^{\infty}$-bound as available in the literature (for this particular system). In [24] the authors carry out the proof in some detail (adapted from a proof in [5]), and show that the properties of the Neumann Heat Kernel involved in it hold for the discrete Neumann Heat Kernel as well; as a consequence, one can emulate the proof in the continuous case to obtain uniform $L^{\infty}$ bounds for the discretized problem. The main idea of the proof is a bootstrapping argument in which the bounds obtained on $a$ and $b$ in terms of $c$ (from the first two equations of (2) we get $a_{t}-k_{a}\Delta a\leq k_{2}c$ and $b_{t}-k_{b}\Delta b\leq k_{2}c$) are fed into the inequality $c_{t}-k_{c}\Delta c\leq k_{1}ab$ (from the third equation of (2)) to yield an $L^{\infty}$ bound on $c$ at some time $t\geq\delta>0$ in terms of a sublinear function of the bound at times $t\leq\delta/2$. The success of this method relies on the right hand sides of the first two equations of (2) being bounded above in the positive orthant by a constant multiple of $c$. It therefore fails for systems with all multiple-species complexes (no single species in any complex) or systems where single species appear more than once in any given complex (such as $A+B{\rightleftharpoons}2C$). These bounds are crucial to the proof of consistency and, ultimately, convergence [24].

Our method to prove these uniform $L^{\infty}$ bounds seems confined to one-D, as it uses the stronger form of Poincaré’s inequality on a bounded interval (where the essential sup norm of a Sobolev function is bounded in terms of its average and the $L^{1}$ norm of its Sobolev derivative). This leads to a uniform estimate (with respect to $t$) of the $L^{2}$ norm of the solutions in cylinders of type $(t,t+1)\times(0,1)$, which, once more using $d=1$, leads to a uniform in time $L^{\infty}$ bound in the case where one of the right hand side polynomials is bounded above by a quadratic polynomial. It is an important improvement that we can deal with the at most quadratic case, since previous results (when boundary equilibria are present) only dealt with two species reactions or, if at least three species are present, the right hand side of an equation from the system is dominated by a first-degree polynomial [10], [11], [12], [13], [16] etc.

1.2 The three-species system.

A case not covered so far in the literature is $A+nB{\rightleftharpoons}B+C$ ($n\geq 2$ is an integer); this has boundary equilibria in some (not all) stoichiometric classes, translates to a $3\times 3$ system ($2\times 2$ being, in general, easier to treat via the standard maximum principle for the heat equation), and the right hand side is not bounded above (in the positive orthant) by a linear term. More precisely, the PDE system we are looking at is

[TABLE]

where $\nu$ is the (outward) normal vector to $\partial\Omega$. Notice that by rescaling time $t$, space $x$ and the concentrations $(a,b,c)$, we can always assume that the reaction rates and the domain volume are 1; to fix ideas, let us choose, without loss of generality, $n=2$ and $k_{1}=k_{2}=1$. The only important restriction we impose is the choice of spatial dimension $d=1$. This is necessary due to the uniform $L^{\infty}$ bound on the solution, which is crucial to our analysis and we can (for the time being) only justify in the $d=1$ case; that is, only the estimates in subsection 2.1 are predicated on this restriction. Therefore, the $3\times 3$ reaction-diffusion system considered here is

[TABLE]

The conserved (in time) quantities here are $\bar{a}+\bar{c}$ and $\bar{b}+\bar{c}$, where $\bar{f}$ denotes the average of the function $f$ over $\Omega$. If $b_{\infty}>0$ we obviously can only have a boundary equilibrium at $(0,b_{\infty},0)$ (i.e. $a_{\infty}=c_{\infty}=0$). The conservation of $\bar{a}+\bar{c}$ forces $a\equiv c\equiv 0$, the second equation of the system decouples into $b_{t}-d_{b}b_{xx}=0$, and $b_{\infty}=\bar{b}_{0}$. This steady state cannot be approached from any initial state for which $\bar{a}_{0}+\bar{c}_{0}>0$, so no initial data in the positive orthant will converge to it. The other nontrivial type of steady states is given by $b_{\infty}=0$ and $a_{\infty}+c_{\infty}=\bar{a}_{0}+\bar{c}_{0}>0$. If $c_{\infty}=0$, we get $b\equiv c\equiv 0$ (from the conservation of $\bar{b}+\bar{c}$), so, once again, an equilibrium of the type $(a_{\infty},0,0)$ cannot be approached from the positive orthant. We are left with the case $(a_{\infty},0,c_{\infty})$ for $a_{\infty}c_{\infty}>0$. We do not know how to prove that $\bar{a}_{0}\bar{b}_{0}\bar{c}_{0}>0$ prevents convergence to such steady states, but in this paper we will prove a weaker statement, namely:

Theorem 1.1**.**

If $a_{0},\ b_{0},c_{0}\in L^{\infty}(0,1)$ are a.e. nonnegative and such that $\bar{a}_{0}\bar{b}_{0}\bar{c}_{0}>0$ and $b_{0}\geq\delta$ a.e. in $(0,1)$ for some $\delta>0$, then the (unique) global classical solution to (4) converges asymptotically exponentially fast (at an explicit rate) to the unique positive equilibrium in its stoichiometric class.

The above theorem will be proved in Section 2.

1.3 The two-species system

Finally, in Section 3 we prove similar results on the convergence to the positive equilibrium for a two-species reversible reaction-diffusion network with boundary equilibria:

[TABLE]

Assume $m_{1}>m_{2}>0$ and $0<n_{1}<n_{2}$ with $\bar{m}=m_{1}-m_{2}$, $\bar{n}=n_{2}-n_{1}$ and $\bar{m}<\bar{n}$. The $2\times 2$ reaction-diffusion system is

[TABLE]

Theorem 1.2**.**

Let $\Omega$ be a bounded domain of $\mathbb{R}^{d}$ with a smooth boundary, for some integer $d\geq 1$. If $0<\alpha\leq a_{0}(x),b_{0}(x)\leq\beta<+\infty$ for a.e. $x$ in $\Omega$, then the (unique) global classical solution to (5) converges asymptotically exponentially (at an explicit rate) to the unique positive equilibrium in its stoichiometric class.

2 Asymptotic decay for the three-species system

We consider the entropy functional $E(a,b,c)$ and the corresponding entropy dissipation (when computed along solutions) $D(a,b,c)=-\frac{d}{dt}{E(a,b,c)}$ associated to the system:

[TABLE]

and

[TABLE]

We would also like to record (for later use) the following conservation laws

[TABLE]

for all $t\geq 0$. Note that these are simply obtained by adding equations 1 and 3 (respectively, 2 and 3) and integrating in space over $[0,1]$ by taking into account the boundary conditions. Note also that $M_{1}$ and $M_{2}$ are finite as long as $a_{0},b_{0},c_{0}\in L^{1}(0,1)$.

2.1 Local $L^{2}$ estimate (a priori estimate)

Proposition 2.1**.**

Let $(a,b,c)$ be a solution for (4) with initial condition $(a_{0},b_{0},c_{0})$ such that $a_{0}>0,\ b_{0}>0,\ c_{0}>0$ a.e. in $[0,1]$ and $a_{0}\ln{a_{0}},b_{0}\ln{b_{0}},c_{0}\ln{c_{0}}\in L^{1}(0,1)$. Then there exists a real constant $C$ such that

[TABLE]

for any $\tau>0$.

Proof.

We start with the obvious inequality (which holds for all $x\in[0,1],\ t>0$)

[TABLE]

then use Hölder’s inequality to get

[TABLE]

Obviously, the above inequalities also hold for $b$ and $c$. Next we integrate the entropy dissipation in time to obtain

[TABLE]

where we have only displayed the dependence on time of the components of the solution vector. Since the last integrand in right hand side of (7) is nonnegative, we conclude

[TABLE]

Since $x(\ln x-1)\geq-1$ for all $x\geq 0$ (at $x=0$ this holds in the limiting sense), we get

[TABLE]

which implies

[TABLE]

for any $t\geq 0$, where $d:=\min\{d_{a},d_{b},d_{c}\}$. Finally, we take into account (8) and (9) to estimate

[TABLE]

Similar inequalities hold for $b$ and $c$, therefore we have finished the proof. ∎

2.2 Uniform $L^{\infty}$ estimate

In this section we shall prove that a classical solution to (4) is bounded uniformly in time (and therefore, it also exists for all time). To achieve this, our goal is to place ourselves in the setting of Theorem 4.1 [17]. We shall refrain from transcribing the assumptions (H1)–(H3) from [17] here, as they are universally satisfied by CRDN systems with nonnegative and essentially bounded initial conditions. On the other hand, assumption (H4’) is both specific to our case and nontrivial to verify. We state it below, as it refers to a general semilinear parabolic $m\times m$ system

[TABLE]

where $x\in\Omega$, $t>0$, $u=(u_{1},...,u_{m})$. It reads:

[TABLE]

[TABLE]

[TABLE]

The following result is a version of Theorem 4.1 [17].

Theorem 2.2**.**

Suppose the initial data $u_{j,0}\in L^{\infty}(\Omega)$, $j=1,...,m$, the generic assumptions (H1)–(H3) from [17] hold. Further assume (11) holds for some $1\leq p<\infty$ and

[TABLE]

Then

[TABLE]

A word of caution is in order: in [17] the analysis is performed on the whole space $\Omega=\mathbb{R}^{d}$ and for more general elliptic operators (instead of the Laplacian). However, we argue that Theorem 2.2 holds for bounded domains $\Omega$ as well. Indeed, our elliptic operator is the Laplacian and on a bounded domain $\Omega\subset\mathbb{R}^{d}$ we have that the estimate (used to prove Lemma 2.1 in [17]) on the fundamental solution to the corresponding parabolic operator

[TABLE]

for $x,\ \xi\in\Omega,\ 0<\tau<t\leq T<\infty$ (where $c_{1}(T),\ c_{2}(T)$ are bounded for finite $T$) holds for the Heat Kernel with Neumann BC on a bounded domain as well (see, e.g, [6]) with $c_{1}(t)=ce^{t}$, and $c_{2}(t)=C$, for some constants $c,\ C>0$. Likewise, for the homogeneous problem with initial value at $T\geq 0$, we also have the estimate

[TABLE]

which is used to prove Lemma 4.1 [17]. In the proof of Lemma 2.2 [17] there is an estimate on solutions of the adjoint equation; [22] is given as a reference. We limit ourselves to noting that [22] covers the case of and provides the same estimate on bounded domains as well. These are the estimates one needs to check in order to convince oneself that Theorem 2.2 holds on a bounded domain $\Omega\subset{R}^{d}$.

If we go back to (4) and denote by

[TABLE]

we see that (11) is satisfied with $\alpha_{1,1}=1,\alpha_{2,1}=-1,\ \alpha_{2,2}=1,\ \alpha_{3,1}=0,\ \alpha_{3,2}=-1,\ \alpha_{3,3}=1$, $r=2$, $K_{1}=k_{2}/2$, $K_{2}=0$. Thus, if we can find $1\leq p<\infty$ such that (12) and the first inequality in (11) (the one bounding $r$ in terms of $p$) are satisfied, we can apply Theorem 2.2 in order to obtain the uniform $L^{\infty}$ bound. But Proposition 2.1 shows that $p=2$ does the job.

The same reference [17] guarantees that the solution is classical, unique and nonnegative. Therefore, we have proved:

Theorem 2.3**.**

If $a_{0},\ b_{0},\ c_{0}\in L^{\infty}(\Omega)$ and are a.e. nonnegative, then a unique classical solution to the system (4) exists for all time. Furthermore, the solution is uniformly (with respect to time) bounded in $L^{\infty}(\Omega)$.

2.3 $L^{1}$ convergence

Let us make the assumption $\beta=\|\frac{1}{b_{0}}\|_{L^{\infty}[0,1]}<\infty$; because the classical solution is continuous, there exists $t_{1}>0$ such that $\|\frac{1}{b(\cdot,t)}\|_{L^{\infty}[0,1]}<10\beta$ for all $t\in[0,t_{1}]$. We next divide the second equation in (4) by $-b^{2}$ and use the uniform (in time) $L^{\infty}$ boundedness of $a$ to get

[TABLE]

Using the maximum principle for the heat equation, we have that, for all $t\in[0,t_{1}]$,

[TABLE]

We can iterate this inequality in time to get

[TABLE]

for all $t>0$. Therefore, we now have an estimate on how fast $b$ can decay to zero.

There exists a unique equilibrium with all positive components for (4) and by (4) and (8) we see that it is given by $v_{\infty}:=(a_{\infty},b_{\infty},c_{\infty})$, where its components are uniquely determined by

[TABLE]

Now we introduce the relative entropy

[TABLE]

and its corresponding entropy dissipation

[TABLE]

where

[TABLE]

At this point we introduce the notation

[TABLE]

On the basis of the following identity

[TABLE]

we get

[TABLE]

The Logarithmic Sobolev Inequality

[TABLE]

(where $C_{LSI}$ only depends on the domain $[0,1]$) yields

[TABLE]

for an explicit constant $C_{2}=\min\{d_{a},d_{b},d_{c}\}\cdot C_{LSI}$. Next, we define two integrand functions:

[TABLE]

[TABLE]

and set $S(a,b,c):=S_{1}(a,b,c)+S_{2}(a,b,c).$ From (17), (19), (21) and (14) we get

[TABLE]

where

[TABLE]

and $\hat{S}$ is the convexification of $S$, i.e. the supremum of all affine functions below $S$. The last inequality above holds due to Jensen’s inequality and the unit volume of the spatial domain.

We next define the compatible class:

[TABLE]

In this class, the first two conditions are related to the conservation laws (8) while the last one follows from the decreasing relative entropy. Since we know $\hat{S}=\widehat{S_{1}+S_{2}}\geq\hat{S_{1}}+\hat{S_{2}}$, $S_{1}$ is convex and $S_{2}$ is non-negative, we have

[TABLE]

Furthermore, it is not hard to verify that

[TABLE]

It follows

[TABLE]

Let

[TABLE]

By (23) and (24) we get $C_{4}$ can only be zero if there exists a sequence $\{v_{n}\}_{n}\subset C_{M_{1},M_{2}}$ such that $v_{n}\rightarrow(a_{\infty},b_{\infty},c_{\infty})$ as $n\rightarrow\infty$. This means

[TABLE]

In order to show that the above limit inferior is positive we use the following lemma [23]:

Lemma 2.4**.**

There exists $\delta>0$ such that for all $v\in B(v_{\infty},\delta)$ (ball centered at $v_{\infty}$ and of radius $delta$) $S(v)$ is locally convex in this ball.

In particular, we get that $\hat{S}\equiv S$ in the ball centered at $(a_{\infty},b_{\infty},c_{\infty})$ with radius $\delta$. Let us now define

[TABLE]

and consider the Taylor expansion of

[TABLE]

around the unique positive equilibrium $(a_{\infty},b_{\infty},c_{\infty})$. Since $a_{\infty}b_{\infty}=c_{\infty}$, we have $D_{2}(a_{\infty},b_{\infty},c_{\infty})=\nabla{D_{2}(a_{\infty},b_{\infty},c_{\infty})}=0$ and quadratic term in the expansion is

[TABLE]

Thus,

[TABLE]

Since $v\in C_{M_{1},M_{2}}$ (which means $x+z=a_{\infty}+c_{\infty},y+z=b_{\infty}+c_{\infty}$), we get

[TABLE]

Then

[TABLE]

Also notice (by direct computation and using that $c_{\infty}=a_{\infty}b_{\infty}$) the identity $D_{2}(v)=bc_{\infty}S_{2}(v)$, which implies (in view of the above inequality)

[TABLE]

Combining the above two steps, we have

[TABLE]

Therefore, in light of (22), we obtain

[TABLE]

so,

[TABLE]

where $C_{5}=\min\{1,C_{4}\}\times\min\{\beta C_{2},a_{\infty}b_{\infty}^{2},b_{\infty}c_{\infty}\}$. Then Gronwall’s lemma yields

[TABLE]

for all $t>0$.

Now we need the following lemma [4]:

Lemma 2.5**.**

For all non-negative and measurable functions $a,b,c:[0,1]\rightarrow\mathbb{R}$ and $\int_{0}^{1}(a+c)=M_{1},\int_{0}^{1}(b+c)=M_{2}$. Then there exists a constant $C_{K}>0$ related only with domain and $M_{1},M_{2}$ such that:

[TABLE]

Therefore, we get

[TABLE]

where $C_{6}=\frac{E(a_{0},b_{0},c_{0}|a_{\infty},b_{\infty},c_{\infty})}{C_{K}}$.

The above inequality shows that the solution stays away from the boundary equilibrium; in fact, its converges to the unique positive equilibrium in the $L^{1}$ norm. In order to show that the convergence rate is, in fact, exponential, we use the above inequality to conclude that there exists a time $T_{\epsilon}$ such that $\|a(t)\|_{1},\|b(t)\|_{1},\|c(t)\|_{1}>\epsilon^{2}>0$ (for some sufficiently small $\epsilon<1$) for all $t>T_{\epsilon}$.

2.4 Entropy entropy-dissipation estimate

By using the inequality (37), we obtain

[TABLE]

where $C_{7}:=\min(4d_{a},4d_{b},4d_{c},a_{\infty}b_{\infty}^{2}+b_{\infty}c_{\infty})$. Due to (8), we have $M:=\max(M_{1},M_{2})$ such that $\overline{a(t)},\overline{b(t)},\overline{c(t)}<M$ for all $t\geq 0$. In what follows we drop the dependence on $t$ from the notation; each time we write $\overline{a}$ or the likes we mean the spatial average of $a(t)=a(\cdot,t)$. Due to (39), we see that

[TABLE]

Since $0<a_{\infty},b_{\infty},c_{\infty}<M$, we have

[TABLE]

where

[TABLE]

Next we claim that there exists a real constant $C_{9}$ such that

[TABLE]

In order to get the above estimate, we introduce the deviations from the mean, i.e. $\delta_{a}=\sqrt{a}-\overline{\sqrt{a}},\delta_{b}=\sqrt{b}-\overline{\sqrt{b}},\delta_{c}=\sqrt{c}-\overline{\sqrt{c}}$. Now we make the decomposition

[TABLE]

where $D_{L}=\{x\in[0,1]:|\delta_{a}|,|\delta_{b}|,|\delta_{c}|\leq L\}$ for a fixed constant $L$. We expand

[TABLE]

and

[TABLE]

to see that on the set $D_{L}$ one has

[TABLE]

and

[TABLE]

where R_{1}:=\big{(}\sqrt{M_{2}}+L\big{)}^{2}+\sqrt{M_{1}}\big{(}2\sqrt{M_{2}}+L\big{)} and R_{2}:=\sqrt{M_{2}}+\big{(}\sqrt{M_{2}}+L\big{)}. Thus,

[TABLE]

where $R(M_{1},M_{2},L):=\frac{4R_{1}^{2}}{a_{\infty}b_{\infty}^{2}}+\frac{4R_{2}^{2}}{b_{\infty}c_{\infty}}$.

On the set $D_{L}^{\complement}$, by using Poincaré’s inequality, we get

[TABLE]

Since

[TABLE]

we infer

[TABLE]

where

[TABLE]

Pick $K>\frac{R+1}{\min\{1,C_{P}\}}$ and combine (28) and (29) to conclude

[TABLE]

where C_{K,R}=\min\big{\{}K,K\tilde{R},\frac{1}{2},KC_{P}-R\big{\}}=\min\big{\{}K\tilde{R},\frac{1}{2}\big{\}} (because $K-R>1$). Therefore, (27) is proved with $C_{9}=\frac{C_{K,R}}{3K}$.

It remains to show that there exists a constant $C_{10}$ such that

[TABLE]

To this end, we introduce $\mu_{a},\mu_{b},\mu_{c}$ to parameterize $\sqrt{\overline{a}},\sqrt{\overline{b}},\sqrt{\overline{c}}$ with $\sqrt{\overline{a}}=\sqrt{a_{\infty}}(1+\mu_{a})$,$\sqrt{\overline{b}}=\sqrt{b_{\infty}}(1+\mu_{b})$,$\sqrt{\overline{c}}=\sqrt{c_{\infty}}(1+\mu_{c})$, where $-1\leq\mu_{a},\mu_{b},\mu_{c}<\mu_{k}$ for $\mu_{k}=\frac{\sqrt{k}}{\min\{\sqrt{a_{\infty}},\sqrt{b_{\infty}},\sqrt{c_{\infty}}\}}-1$. Since $\delta_{a}=\sqrt{a}-\overline{\sqrt{a}}$, we have

[TABLE]

where $T(a)=\frac{1}{\sqrt{\overline{a}}+\overline{\sqrt{a}}}\leq\frac{1}{\epsilon}$; this inequality follows from $\overline{a}=\|a\|_{1}>\epsilon^{2}>0$. Similarly,

[TABLE]

where $T(b)=\frac{1}{\sqrt{\overline{b}}+\overline{\sqrt{b}}},T(c)=\frac{1}{\sqrt{\overline{c}}+\overline{\sqrt{c}}}\leq\frac{1}{\epsilon}$. And since

[TABLE]

due to this lower bound on $\overline{\sqrt{b}}$, we can factor out $(\frac{\overline{\sqrt{b}}}{\sqrt{b_{\infty}}})^{2}$ and reduce (4) to the system associated with the reversible reaction $A+B{\rightleftharpoons}C$ (which does not have boundary equilibria). More precisely, we have

[TABLE]

Since \|\delta_{a}\|_{2}^{2}=\overline{a}-\big{(}\overline{\sqrt{a}}\big{)}^{2}, we get $\|\delta_{a}\|_{2}^{2}\leq k$; similarly, $\|\delta_{b}\|_{2}^{2},\|\delta_{c}\|_{2}^{2}\leq k$. Combined with $T(a),T(b),$ $T(c)\leq\frac{1}{\epsilon}$ and $0\leq 1+\mu_{a},\ 1+\mu_{b},\ 1+\mu_{c}\leq 1+\mu_{k}$, (31) gives

[TABLE]

where

[TABLE]

Poincaré’s inequality reveals

[TABLE]

for $C_{12}:=\frac{C_{P}\epsilon^{4}}{2b_{\infty}kC_{11}}$. On the other hand,

[TABLE]

so we would like to compare $[(1+\mu_{a})(1+\mu_{b})-(1+\mu_{c})]^{2}$ and $\mu_{a}^{2}+\mu_{b}^{2}+\mu_{c}^{2}$. The conservation laws (8) (applied to $\overline{a},\ \overline{b},\ \overline{c}$ and $a_{\infty}$, $b_{\infty}$, $c_{\infty}$) yield

[TABLE]

and so, $a_{\infty}\mu_{a}+c_{\infty}\mu_{c}=b_{\infty}\mu_{b}+c_{\infty}\mu_{c}=0$. Thus, unless $\mu_{a},\mu_{b},\mu_{c}$ are all zero (trivial case!), we get $\mu_{a}\mu_{c}<0$ and $\mu_{b}\mu_{c}<0$. If $\mu_{a},\mu_{b}>0$ and $0>\mu_{c}$, then

[TABLE]

Otherwise, if $\mu_{a},\mu_{b}<0$ and $0<\mu_{c}$, then

[TABLE]

Therefore, in both cases we have

[TABLE]

(Notice that when $\mu_{a}=\mu_{b}=\mu_{c}=0$, both sides of the inequality are equal to zero.) Set

[TABLE]

to conclude the proof of (30).

Proof of Theorem 1.1:

Proof.

Finally, by (25), (26), (27) and (30), we obtain

[TABLE]

In view of the above inequality and (21), we discover

[TABLE]

where

[TABLE]

In conclusion, we have proved that the solution decays exponentially to the positive equilibrium (with explicit rate). ∎

3 Asymptotic decay for the two-species system

In this section we prove Theorem 1.2.

3.1 Uniform boundedness and global existence for the two-species system

To show the uniform boundedness for classical solutions to (5) we estimate the $L^{p}$ norm and pass to the limit as $p\rightarrow\infty$. For the lower bound we have

[TABLE]

Similarly, we have

[TABLE]

Therefore we can let $p>1$ be sufficiently large so that $q=\frac{(p+1)\bar{n}}{\bar{m}}-1>1$; then

[TABLE]

This gives us

[TABLE]

[TABLE]

Since

[TABLE]

as we we let $p\rightarrow\infty$, we have

[TABLE]

Also we have a similar result on b, i.e.

[TABLE]

[TABLE]

Again, as $p\rightarrow\infty$, we have

[TABLE]

Thus, we get the uniform lower bound $\epsilon^{2}$ for $a(x,t)$ and $b(x,t)$ with $\epsilon^{2}:=\min(L_{a}^{-1},L_{b}^{-1})$. As for the upper bound, we use the same method to see that

[TABLE]

Similarly, we have

[TABLE]

Therefore, we can let $p>1$ be sufficient large to render $q=\frac{(p-1)\bar{n}}{\bar{m}}+1>1$; then

[TABLE]

This gives us

[TABLE]

[TABLE]

But

[TABLE]

so, when we let $p\rightarrow\infty$, we obtain

[TABLE]

We have a similar estimate on $b$,

[TABLE]

[TABLE]

Again, as $p\rightarrow\infty$, we get

[TABLE]

Thus, we get the uniform upper bound $\omega$ for $a(x,t)$ and $b(x,t)$ with $\omega:=\max(U_{a},U_{b})$, and this also implies the existence of a unique global classical solution.

3.2 Convergence for the two species system

We use the same entropy entropy dissipation method to obtain an explicit exponential convergence rate for the two species system in any dimension.

Again we introduce the relative entropy

[TABLE]

and its corresponding entropy dissipation

[TABLE]

Due to the following identity

[TABLE]

and the Logarithmic Sobolev Inequality (20) we have

[TABLE]

where $D_{1}=\min(d_{a},d_{b})\cdot C_{LSI}$.

From inequality (37) we get the following estimate

[TABLE]

where $D_{2}=\min(4d_{a},4d_{b},4d_{c},a_{\infty}^{m_{1}}b_{\infty}^{n_{1}}+a_{\infty}^{m_{2}}b_{\infty}^{n_{2}})$.

The obvious conservation law

[TABLE]

shows there exists $N$ such that $\overline{a},\overline{b}<N$. Thus, (39) yields

[TABLE]

Since $0<a_{\infty},b_{\infty}<N$,

[TABLE]

where

[TABLE]

Now we claim there exists a constant $D_{4}>0$ such that

[TABLE]

Now we again introduce the deviations $\delta_{a}=\sqrt{a}-\overline{\sqrt{a}},\delta_{b}=\sqrt{b}-\overline{\sqrt{b}}$ and make the decomposition

[TABLE]

where $D_{L}:=\{x\in\Omega\,:|\,\delta_{a}|,|\delta_{b}|\leq L\}$ with a fixed constant $L$. On the set $D_{L}$ we get

[TABLE]

[TABLE]

where $R_{1}$ and $R_{2}$ are finite due to the boundedness of $|\delta_{a}|,|\delta_{b}|,\overline{\sqrt{a}},\overline{\sqrt{b}}$. Then we get

[TABLE]

where $R(|\delta_{a}|,|\delta_{b}|,\overline{\sqrt{a}},\overline{\sqrt{b}})=\frac{4R_{1}^{2}}{a_{\infty}^{m_{1}}b_{\infty}^{n_{1}}}+\frac{4R_{2}^{2}}{a_{\infty}^{m_{2}}b_{\infty}^{n_{2}}}$ is finite (depends on the choice of $L$ and $N$).

On the set $D_{L}^{\complement}$, by using Poincaré’s inequality, we get

[TABLE]

Since

[TABLE]

we infer

[TABLE]

where

[TABLE]

We combine the above two parts, pick $K>\frac{R+1}{\min\{1,C_{P}\}}$ and have the following

[TABLE]

where $C_{K,R}:=\min\{K\tilde{R},\frac{1}{2}\}$. We can fix $L>0$ and get the corresponding $R$, then pick sufficiently large $K$ (e.g. $K>R+1$) such that we obtain (35) with $D_{4}=\frac{C_{K,R}}{3K}$.

It remains to show that there exists a constant $D_{5}$ such that

[TABLE]

We again introduce $\mu_{a},\mu_{b}$ to parameterize $\sqrt{\overline{a}},\sqrt{\overline{b}}$ with $\sqrt{\overline{a}}=\sqrt{a_{\infty}}(1+\mu_{a})$,$\sqrt{\overline{b}}=\sqrt{b_{\infty}}(1+\mu_{b})$, where $\mu_{\epsilon}\leq\mu_{a},\mu_{b}<\mu_{\omega}$ with $\mu_{\epsilon}=\frac{\epsilon}{\max\{\sqrt{a_{\infty}},\sqrt{b_{\infty}}\}}-1$ and $\mu_{\omega}=\frac{\sqrt{\omega}}{\min\{\sqrt{a_{\infty}},\sqrt{b_{\infty}}\}}-1$. We have $\overline{\sqrt{a}}=-\frac{\|\delta_{a}\|_{2}^{2}}{\sqrt{\overline{a}}+\overline{\sqrt{a}}}+\sqrt{\overline{a}}=\sqrt{\overline{a}}-T(a)\|\delta_{a}\|_{2}^{2}$ , where $T(a)=\frac{1}{\sqrt{\overline{a}}+\overline{\sqrt{a}}}$. Similarly, $\overline{\sqrt{b}}=\sqrt{\overline{b}}-T(b)\|\delta_{b}\|_{2}^{2}$ , where $T(b)=\frac{1}{\sqrt{\overline{b}}+\overline{\sqrt{b}}}$. Both $T(a),T(b)$ have uniform (in time) upper and lower bounds. Since we get the lower bound for $\overline{\sqrt{a}},\overline{\sqrt{b}}$, we can factor out \bigg{(}\frac{{\overline{\sqrt{a}}}^{m_{2}}{\overline{\sqrt{b}}}^{n_{1}}}{\sqrt{a_{\infty}}^{m_{2}}\sqrt{b_{\infty}}^{n_{1}}}\bigg{)}^{2} to reduce the original system (5) to the system associated with the reaction $\bar{m}A{\rightleftharpoons}\bar{n}B$ (which does not have boundary equilibria). Then we get

[TABLE]

We evaluate

[TABLE]

where $S=\max{(|S_{1}|,|S_{2}|)}$.

We have $\|\delta_{a}\|_{2}^{2}=\overline{a}-(\overline{\sqrt{a}})^{2}\leq\omega$, similarly $\|\delta_{b}\|_{2}^{2}\leq\omega$ and $T(a),T(b)\leq\frac{1}{\epsilon}$, $\mu_{a},\mu_{b}<\mu_{\omega}$; so $S$ is uniformly bounded and

[TABLE]

where $D_{6}=4\frac{\epsilon^{2(m_{2}+n_{1})}}{a_{\infty}^{m_{2}}b_{\infty}^{n_{1}}}\|S\|_{\infty}$. Poincaré’s inequality yields

[TABLE]

where $D_{7}=\frac{\epsilon^{2(m_{2}+n_{1})}}{2a_{\infty}^{m_{2}}b_{\infty}^{n_{1}}}D_{6}^{-1}$. On the other hand,

[TABLE]

so we need to compare $[(1+\mu_{a})^{\bar{m}}-(1+\mu_{b})^{\bar{n}}]^{2}$ with $\mu_{a}^{2}+\mu_{b}^{2}$. The conservation law gives

[TABLE]

so, as long as $\mu_{a},\mu_{b}$ are nonzero, $\mu_{a}$ and $\mu_{b}$ must have different signs. If $\mu_{a}>0>\mu_{b}$, then

[TABLE]

Otherwise,

[TABLE]

so in both cases we have

[TABLE]

Notice that when $\mu_{a}=\mu_{b}=0$, both sides of the inequality are equal to zero.

Proof of Theorem 1.2:

Proof.

Set $D_{5}=\frac{D_{7}}{\max(a_{\infty},b_{\infty})}$ to see that (36) holds, and then combine (33), (34), (35) and (36) to reveal

[TABLE]

In view of (32), we get

[TABLE]

for $D_{8}=\min(\frac{D_{2}D_{4}D_{5}}{D_{3}},D_{1})$, which finally proves that the solution decays exponentially to the positive equilibrium at an explicit rate. ∎

4 Remarks and open problems

We believe this approach works to prove and quantify the decay rate to the complex balanced equilibrium for more general systems of the type

[TABLE]

and even more complex systems such as

[TABLE]

where $k\geq 1$ and $l,\ m\geq 2$ are integers.

In order to obtain the uniform essential bound on the densities we used the $L^{\infty}$ version of Poincaré’s inequality, which is only available in 1D. Can we employ more refined techniques in order to prove the essential bounds in higher dimensions? These are some questions we plan to address in future work.

5 Appendix

Lemma 5.1**.**

For any $x,\ y>0$ we have

[TABLE]

Proof.

The case $x=y$ is trivial. If $x>y>0$ we use the Jensen inequality for the convex function $f(s):=[y+s(x-y)]^{-1}$ to get

[TABLE]

By using this inequality and $x+y>2\sqrt{xy}$, we conclude

[TABLE]

Or suppose we have $y>x>0$ and set $g(u)=e^{u}$, $u(s)=a+s(b-a)$, for $b>a$. Jensen’s inequality shows

[TABLE]

Let $b=\ln y,\ a=\ln x$ to deduce

[TABLE]

which implies

[TABLE]

It follows

[TABLE]

∎

The second important tool is:

Lemma 5.2**.**

For each fixed $y>0$,

[TABLE]

Proof.

Suppose $x>y>0$ and set $k=\sqrt{\frac{x}{y}}>1$. We have

[TABLE]

Since we have $k>1$, we use (38) to obtain

[TABLE]

Combine this with $\sqrt{x}-\sqrt{y}\geq 0$ and $y>0$ to get

[TABLE]

We can use the same way to show this inequality is correct when $y>x>0$. ∎

6 Acknowledgments

Gheorghe Craciun and Jiaxin Jin acknowledge support from the National Science Foundation under grants DMS–1412643 and DMS–1816238. Casian Pantea was partially supported by National Science Foundation grant DMS–1517577. Adrian Tudorascu acknowledges support from the National Science Foundation under grant DMS–1600272, and from the MSRI Research Professorship in the program “Hamiltonian systems, from topology to applications through analysis” from the Fall of 2018.