The entropy method for reaction-diffusion systems without detailed balance: first order chemical reaction networks

In this paper, the applicability of the entropy method for the trend towards equilibrium for reaction-diffusion systems arising from first order chemical reaction networks is studied. In particular, we present a suitable entropy structure for weakly reversible reaction networks without detail balance condition. We show by deriving an entropy-entropy dissipation estimate that for any weakly reversible network each solution trajectory converges exponentially fast to the unique positive equilibrium with computable rates. This convergence is shown to be true even in cases when the diffusion coefficients all but one species are zero. For non-weakly reversible networks consisting of source, transmission and target components, it is shown that species belonging to a source or transmission component decay to zero exponentially fast while species belonging to a target component converge to the corresponding positive equilibria, which are determined by the dynamics of the target component and the mass injected from other components. The results of this work, in some sense, complete the picture of trend to equilibrium for first order chemical reaction networks.


Introduction and Main results
This paper investigates the applicability of the entropy method and proves the convergence to equilibrium for reaction-diffusion systems, which do not satisfy a detailed balance condition.
The mathematical theory of (spatially homogeneous) chemical reaction networks goes back to the pioneer works of e.g. Horn, Jackson, Feinberg and the Volperts, see [Fei79,Fei87,FH,Hor72,Hor74,HJ72,Vol,VVV] and the references therein. The aim is to study the dynamical system behaviour of reaction networks independently of the values of the reaction rates. It is conjectured since the early of 1970s that in a complex balanced system, the trajectories of the corresponding dynamical system always converge to a positive equilibrium. This conjecture was given the name Global Attractor Conjecture by Craciun et al. [CDSS]. The conjecture in its full generality is -up to our knowledge -still unsolved so far, despite many attempts have been made by mathematicians to attack this problem.
From the many previous works concerning the large time behaviour of chemical reaction networks, the majority of the existing results considers the spatially homogeneous ODE setting. The PDE setting in terms of reactiondiffusion systems is less studied. Also detailed quantitative statements like, e.g. rates of convergence to equilibrium, constitute frequently open questions even in the ODE setting.
Our general aim is to prove quantitative results on the large-time behaviour of chemical reaction networks modelled by reaction-diffusion systems. In the present work, we study reaction-diffusion systems arising from first order chemical reaction networks and show that all solution trajectories converge exponentially to corresponding equilibria with explicitly computable rates.
Our approach applies the so called entropy method. Going back to ideas of Boltzmann and Grad, the fundamental idea of the entropy method is to quantify the monotone decay of a suitable entropy (e.g. a convex Lyapunov) functional in terms of a functional inequality connecting the time-derivative of the entropy, the so called entropy dissipation functional, back to the entropy functional itself, i.e. to derive a so called entropy entropy-dissipation (EED) inequality. Such an EED inequality can only hold provided that all conservation laws are taken into account. After having established an EED inequality and applying it to global solutions of a dissipative evolutionary problem, a direct Gronwall argument implies convergence to equilibrium in relative entropy with rates and constants, which can be made explicit.
By being based on functional inequalities (rather than on direct estimates on the solutions), a major advantage of the entropy method is its robustness with respect to model variations and generalisations. Moreover, the entropy method is per se a fully nonlinear approach.
The fundamental idea of the entropy method originates from the pioneer works of kinetic theory and from names like Boltzmann and Grad in order to investigate the trend to equilibrium of e.g. models of gas kinetics.
A systematic effort in developing the entropy method for dissipative evolution equations started not until much later, see e.g. the seminal works [Tos,TV,CJMTU,AMTU,DV01] and the references therein for scalar (nonlinear) diffusion or Fokker-Planck equations, and in particular the paper of Desvillettes and Villani concerning the trend to equilibrium for the spatial inhomogeneous Boltzmann equation [DV05]. The derivation of EED inequalities for scalar evolution equations is typically based on the Barky-Emery strategy (see e.g. [CJMTU, AMTU]), which seems to fail (or be too involved) to apply to systems.
The great challenge of the entropy method for systems is, therefore, to be able to derive an entropy entropy-dissipation inequality, which summarises (in the sense of measuring with a convex entropy functional) the entire dissipative behaviour of solutions to a (possibly nonlinear) dynamical system to which the EED inequality shall be applied to. Preliminary results based on a (nonexplicit) compactness-contradiction argument in 2D were obtained e.g. in [Grö, GGH, GH] in the context semiconductor drift-diffusion models.
The first proof of an EED inequality with explicitly computable constants and rates for specific nonlinear reaction-diffusion systems was shown in [DF06] and followed by e.g. [DF07,DF08,DF15,FLT,MHM]. The application of these EED inequalities to global solutions of the corresponding reactiondiffusion systems proves (together with Csiszár-Kullback-Pinsker type inequalities) the explicit convergence to equilibrium for these reaction-diffusion systems.
We emphasise that all these previous results on entropy methods for systems assumed a detailed balance condition and, thus, features the free energy functional as a natural convex Lyapunov functional.
A main novelty of the paper lies in demonstrating how the entropy method can be generalised to first order reaction networks without detailed balance equilibria. In particular we shall consider first weakly reversible networks and secondly even more general composite systems consisting of source, transmission and target components (see below for the precise definitions).
We feel that it is important to point out that while there are certainly many classical approaches by which linear reaction-diffusion systems can be successfully dealt with, our task at hand is to clarify the entropic structure and the applicability of the entropy method for linear reaction networks as a first step before being able to turn to nonlinear problems in the future.
The goal of this present work is to prove the explicit convergence to equilibrium for the complex balanced (and more general) systems corresponding to first order reaction networks. To be more precise, we study a first order reaction network of the form a ji a ij Figure 1. A first-order chemical reaction network where S i , i = 1, 2, . . . , N, are different chemical substances (or species) and a ij , a ji ≥ 0 are constant reaction rates. In particular, a ij denotes the reaction rates from the species S j to S i . First order reaction networks appear in many classical models, see e.g. [Smo, Rot]. More recently, first order catalytic reactions are used to model transcription and translation of genes in [TVO]. The evolution of the surface morphology during epitaxial growth involves the nucleation and growth of atomic islands, and these processes may be described by first order adsorption and desorption reactions coupled with diffusion along the surface. A first order reaction network can also be used to describe the reversible transitions between various conformational states of proteins (see e.g. [MGetal]). RNA also exists in several conformations, and the transitions between various folding states follow first order kinetics (see [BRetal]).
In the present paper, we investigate the entropy method and the large time behaviour of reaction-diffusion systems modelling first order reaction networks with mass action kinetics. More precisely, we shall consider the reaction network N in the context of reaction-diffusion equations and assume that for all i = 1, 2, . . . , N the substances S i are described by spatial-temporal concentrations u i (x, t) at position x ∈ Ω and time t ≥ 0. Here, Ω shall denote a bounded domain Ω ⊂ R n with sufficiently smooth boundary ∂Ω (that is ∂Ω ∈ C 2+α to avoid all difficulties with boundary regularity, although the below methods should equally work under weaker assumptions) and the outer unit normal ν(x) for all x ∈ ∂Ω. Due to the rescaling x → |Ω| 1/n x, we can moreover consider (without loss of generality) domains with normalised volume, i.e. |Ω| = 1.
In addition, we assume that each substance S i diffuses with a diffusion rate d i ≥ 0 for all i = 1, 2, . . . , N. Finally, we shall assume mass action law kinetics as model for the reaction rates, which leads to the following linear reactiondiffusion system: where X(x, t) = [u 1 (x, t), u 2 (x, t), . . . , u N (x, t)] T denotes the vector of concentrations subject to non-negative initial conditions X 0 (x) = [u 1,0 (x) ≥ 0, u 2,0 (x) ≥ 0, . . . , u N,0 (x) ≥ 0] T , D = diag(d 1 , d 2 , . . . , d N ) denotes the diagonal diffusion matrix and the reaction matrix A = (a ij ) ∈ R N ×N satisfies the following conditions: The conditions (1.2) on the reaction matrix A imply in particular that the vector (1, 1, . . . , 1) T constitutes a left-eigenvector corresponding to the eigenvalue zero. Thus, the homogeneous Neumann boundary conditions guarantee that solutions to (1.1) admit the following conservation of total mass : where M > 0 is the initial total mass, which we shall assume positive.

Definition 1.1 (Detailed Balance Equilibrium).
A positive state X ∞ = (u 1,∞ , u 2,∞ , . . . , u N,∞ ) > 0 is called a detailed balance equilibrium for the reaction network N if a positive reaction rate constant a ij > 0 for i = j implies also a positive reversed reaction rate constant a ji > 0 and that the forward and backward reaction rates balance in equilibrium, i.e.
The reaction network N is called to satisfy detailed balance condition if it admits a detailed balance equilibrium.

Definition 1.2 (Complex Balance Equilibrium).
A positive state X ∞ = (u 1,∞ , u 2,∞ , . . . , u N,∞ ) > 0 is called a complex balance equilibrium for the reaction network N if for all k = 1, 2, . . . , N, the total inflow into the substance S k balances in equilibrium the total out-flow from S k to all other substances S i , i.e.
The reaction network N is called complex balanced if it admits a complex balance equilibrium.
Example 1.1 (Complex and detailed balance equilibria).
It is easy to see that detailed balance equilibria are also complex balance equilibria while the reverse does not hold in general, even in reversible networks. For example, consider the reversible reaction network in Figure 2, where all a 12 a 31 a 13 a 23 a 32 Figure 2. A reversible network reaction rates constants a ij > 0 are assumed positive and the network is thus fully reversible. The corresponding reaction-diffusion system with homogeneous Neumann boundary conditions 31 )u 1 + a 12 u 2 + a 13 u 3 , ∂ t u 2 − d 2 ∆u 2 = a 21 u 1 − (a 12 + a 32 )u 2 + a 23 u 3 , ∂ t u 3 − d 3 ∆u 3 = a 31 u 1 + a 32 u 2 − (a 13 + a 23 )u 3 , ∂ ν u 1 = ∂ ν u 2 = ∂ ν u 3 = 0.
For X ∞ to be a detailed balance equilibrium, it is however additionally necessary that      a 12 u 2,∞ = a 21 u 1,∞ , a 32 u 2,∞ = a 23 u 3,∞ , a 13 u 3,∞ = a 31 u 1,∞ , (1.6) which obviously implies (1.5). Yet the equations (1.6) can only have a solution if a 12 · a 23 · a 31 a 21 · a 32 · a 13 = 1, (1.7) holds; in other words if the product of the reaction rate constants multiplied in the clockwise sense of the above reaction network graph equals the product of the reaction rate constants multiplied in the counterclockwise sense. The condition (1.7) is thus necessary and sufficient for system (1.4) to admit a detailed balance equilibrium.
Remark 1.1 (General definition of detailed and complex balance). The concepts of detailed balance and complex balance are also defined for general higher order chemical reaction networks, see e.g. [Hor72]. For simplicity, we choose here the definition corresponding to the first order network N . Roughly speaking, a state X ∞ is called a complex balanced equilibrium if at equilibrium the in-flow to each specie S i is equal to the out-flow from S i .

Remark 1.2 (Detailed balance and reversability).
It follows from Definition (1.1) that if N satisfies the detailed balance condition, then it is also reversible in the sense that for any reaction S i → S j also the reverse reaction S j → S i takes place.
The set of complex balanced systems is much larger than the one of detailed balance systems. Horn already gave necessary and sufficient conditions for a network to satisfy the complex balance condition in [Hor72]. For convenience of the reader, we shall present in the following the associated definitions based on directed graphs as representations of reaction networks.
Given the reaction network N , one defines a corresponding directed graph G by considering the substances S i , i = 1, 2, . . . , N, as the N nodes of G, which are connected for all i = j = 1, 2, . . . , N by an edge with starting node S i and finishing node S j if and only if the reaction S i a ji − → S j occurs with a positive reaction rate constant a ji > 0.

Definition 1.3 (Weak reversability).
A directed graph G is called weakly reversible if for any edge S i → S j with i = j, there exists a path S j ≡ S j 1 → S j 2 → . . . → S jr ≡ S i where S j 1 , S j 2 , . . . , S jr are other nodes of G.
A reaction network is call weakly reversible if the corresponding graph is weakly reversible.

Definition 1.4 (Strongly connected components).
A subset H ⊂ G is called a strongly connected component of a graph G if for any two nodes S i , S j in H, we can find a path from S i to S j of the form It's easy to see that if G is weakly reversible then G consists of disjoint strongly connected components, i.e. there are no edges connecting two nodes of different strongly connected components.

Remark 1.3 (Separability of strongly connected components in first order networks).
For reaction networks of first order, each node represents exactly one substance. Thus, a first order reaction network with multiple disjoint strongly connected components can be split into the corresponding subsystems, which can be studied independently. For higher order networks, this is not necessarily true since one substance might need to be represented by different nodes.
Because of Remark 1.3, we will only consider connected (yet not necessarily strongly connected) reaction networks throughout this paper.
The particular case of weakly reversible networks plays an important role in our analysis. The following definition provides a characterisation of the weak reversibility in terms of reaction matrices: Definition 1.5 (Decomposable and indecomposable reaction matrices). A reaction matrix A is called decomposable if there exists a permutation matrix P such that The following lemma shows that weakly reversible reaction networks subject to positive initial mass have a unique positive equilibrium. Lemma 1.2 (Indecomposable reaction matrices and complex balance equilibria). Assume that the initial total mass M as defined in the conservation law (1.3) is positive.
Then, the first order reaction network N is complex balanced if and only if A is indecomposable. Moreover, the complex balance equilibrium X ∞ = (u 1,∞ , u 2,∞ , . . . , u N,∞ ) of system (1.1) satisfies (1.9) Proof. Proposition 2.5 Remark 1.4 (Complex balanced higher order systems are necessarily weakly reversible). For higher order reaction network, it holds only true that systems with complex balance equilibrium are necessarily weakly reversible. Thus, weakly reversible systems constitute the more general class of reaction networks.
Remark 1.5. The equilibrium X ∞ in (1.9) is spatially homogeneous. Thus, it coincides with the equilibrium for the corresponding spatially homogeneous ODE system X t = AX of the reaction network given in Figure 1. In [And] or [SiMa], the authors proved that X(t) −→ X ∞ as t −→ +∞. However, the method used in this paper cannot be directly applied to prove the convergence to equilibrium for PDE system (1.1). In addition, the rate of convergence to equilibrium is unknown even in ODE cases.
The first main result of this paper concerns the trend to equilibrium for weakly reversible reaction networks of the form displayed in Figure 1. In the following, we will apply the entropy method to prove explicit exponential convergence of solutions of system (1.1) to the unique equilibrium.
As mentioned above, all previous results of explicit EED inequalities (see e.g. [DF06,DF07,DF08,DF15,FLT,MHM]) considered reaction-diffusion systems satisfying a detailed balance condition. Thus, the free energy functional of the form provided a natural Boltzmann-type entropy functional. We refer e.g. to [BD] for a formal derivation of a prototypical reaction-diffusion system with quadratic mass action law reaction rates from an underlying kinetic systems of reacting Boltzmann equations in terms of a macroscopic limit. However, considering the above system (1.4), a direct computation shows that the free energy functional is only a Lyapunov functional if and only if the detailed balance condition (1.7) holds. Thus, for all reaction rates for which the above system (1.4) is a complex balanced system, the free energy functional fails to provide a suitable entropy structure to be exploited via the entropy method. In fact, one can check that given suitable rate constants a ij and initial data, any convex functional of the form (1.10) with Φ i : [0, ∞) → R being convex C 1 functions with Φ(0) = 0, might transiently increase in time. In particular, this is also true for the quadratic free energy functional Φ i ∼ u 2 i , which would be the most natural convex entropy functional to apply to linear systems.
An algebraic reason for this failure can be seen in a certain incompatibility of the eigenstructures of diffusion and reaction in (1.1). To illustrate the problem, consider a diagonalisable reaction matrix A: In the special ODE case (i.e. D = 0) in (1.1), a simple transformation onto the eigenbasis of A →Ã = diag(0, . . . , 0, λ < 0, . . .) uncovers the mass conservation laws as well as the exponential decay of all non-conserved eigenmodes.
In the PDE case, however, this transformation yields a non-diagonal diffusion matrixD. Thus, by testing the transformed reaction-diffusion system (1.1) withX T from the left, integration over Ω will not result in a quadratic Lyapunov functional of the form (1.10) since the quadratic form ∇X TD ∇X is indefinite in general.
On the other hand, by observing that the diffusion matrix D in (1.1) is already in diagonal form and by testing (1.1) with X T from the left, the positive definitness of the quadratic form ∇X T D∇X ≥ 0 (after integration by parts) is trivial, yet the reaction quadratic form X T AX is indefinite except when A describes a detailed balance system.
An algebraic solution and characterisation of this difficulties arising for complex balance systems is currently work in preparation in [FPT]. The key idea is to look for a positive, diagonal multiplier matrix B such that the multiplied diffusion quadratic form ∇X T BD∇X ≥ 0 is positive definite and the multiplied reaction quadratic form X T BAX ≤ 0 is simultaneously negative semi-definite. In [FPT], we are able to show that such a positive, diagonal multiplier matrix B exists and B has entries given by the inverses of certain co-factors.
In order to interpret the algebraic considerations of [FPT], we also observed at the hand of example systems, that the multiplication with the positive, diagonal multiplier matrix B suggests to consider the following quadratic relative entropy functional (1.11) The first key result of this paper (see Lemma 2.6 below) proves that the relative entropy (1.11) decays monotone in time according to the following explicit form (1.12) Remark 1.6 (Relative entropies necessary for complex balance systems?).
For detailed balance systems and their natural entropy (free energy) functional it is straightforward to calculate the corresponding entropy dissipation functional. For further details, readers are referred e.g. to [MHM].
For complex balanced systems, however, the dissipative structure (1.11) and (1.12) is the result of non-trivial calculations (see Lemma 2.6 below) and seems to be -up to our knowledge -a novel results, cf. Remark 1.8. The entropy structure seems to suggest that for complex balanced systems an admissible entropy functional needs to also include (information on) the complex balanced equilibrium X ∞ .
This point, however, is still unclear since the algebraic considerations in [FPT] seem also to suggest that the above relative entropy with respect to the equilibrium might not be unique and it is an open problem if alternative entropy structures exist, which do not contain the complex balanced equilibrium X ∞ .
Remark 1.7 (Constructive entropy structure). The above entropy (1.11) and entropy-dissipation (1.12) functionals provide an explicit entropy structure exhibiting how weakly reversible linear reaction networks dissipate the entropy (1.11) and, as a consequence, converge exponentially to a unique positive equilibrium. Together with the below results on general first order reaction networks, the results of the present paper provide a constructive and more general improvement of the abstract previous results in [DFM], where linear reaction-diffusion type systems were considered under suitable, yet non-constructive assumptions on diffusion and reaction matrices.

Remark 1.8 (Lyapunov functionals for ODE systems).
For complex ODE systems, L 1 -type Lyapunov functionals are most commonly used in the study of the large-time-behaviour. For reaction-diffusion systems, however, L 1 -functionals are not useful for the entropy method and proving explicit convergence to equilibrium, since they do not measure the spatial diffusion in an exploitable way.
We also remark, that while relative entropy functionals of the form were known to constitute a monotone decaying Lyapunov functional for complex balanced ODE reaction networks (see e.g. [HJ72,Gop,SiMa]), up to our knowledge and somewhat surprisingly, no explicit expression of the entropy dissipation −dV /dt in complex balanced systems has been derived so far. We also refer the reader to e.g. [MiSi] for the stability of some mass action law reaction-diffusion systems, where the author used techniques of ω-limit sets along with the monotonicity of L 1 -type Lyapunov functional.
Our results in this paper are significantly stronger in the sense that we show, by using the entropy method, the exponential convergence to equilibrium with computable rates.
In addition and in comparison to ω-limit techniques, the entropy method has also the major advantage of relying on functional inequalities rather than on specific estimates of solutions to a given system. Having such functional entropy entropy-dissipation inequalities once and for all established makes the entropy method robust with respect to model variations and generalisations.
As example, it is the intrinsic robustness of the entropy method, which makes it possibly to also apply to non weakly reversible reaction networks, see Theorems 1.5 and 1.6 below.
With the help of the explicit form of entropy dissipation (1.12), we are able to show in Lemma 2.7 below an entropy-entropy dissipation inequality of the form where λ > 0 is an explicitly computable constant. By means of the EED inequality (1.14), we state the first main result of this paper to be proved in Section 2 below: Theorem 1.3 (Exponential equilibration of weakly reversible linear reaction networks). Assume that the reaction network in Figure 1 is weakly reversible, the diffusion coefficients d i are positive for all i = 1, 2, . . . , N, and the initial mass M is positive.
Then, the unique global solution to initial-boundary problem (1.1) converges exponentially to the unique positive equilibrium where the constant λ > 0 can be computed explicitly.
The assumption on the positivity of all diffusion coefficients in Theorem 1.3 is not necessary as such. As already shown in e.g. [DF07,FLT], the combined effect of diffusion of a specie and its weakly reversible reaction with other (possibly non-diffusive) species will lead to a indirect "diffusion-effect" on the latter specie. This indirect diffusion-effect can also be measured in terms of functional inequalities. Hence the exponential convergence to equilibrium still holds for systems with partial degenerate diffusion.
The exponential convergence for weakly reversible systems (1.1) with degenerate diffusion is stated in the following Theorem 1.4 to be proved in Section 2 below: Theorem 1.4 (Equilibration of linear networks with degenerate diffusion). Assume that the reaction network in Figure 1 is weakly reversible and the initial mass M is positive. Moreover, assume that at least one diffusion coefficient d i is positive for some i = 1, 2, . . . , N.
Then, the solution to (1.1) converges exponentially fast to the unique positive equilibrium Remark 1.9 (Same results of linear ODE reaction networks).
We remark that our approach can of course be adapted to equally apply to linear ODE reaction networks by eliminating the terms and calculations concerning spatial diffusion. Thus, all the results of this paper hold equally for such linear ODE systems.
As the second main result of this manuscript, we shall derive an entropy approach and prove convergence to equilibrium for reaction networks as in Figure 1, for which the weak reversibility assumption does not hold. For first order reaction networks, this implies that the system is not complex balanced, or in other words, that equilibria are not necessarily positive.
Due to the lack of positivity of equilibria, it follows immediately that the relative entropy used for weakly reversible systems is not directly applicable. In the following we proposed a modified entropy approach. At first, it is necessary to understand the structure of non weakly reversible reaction networks.
We state here the necessary terminology and the main ideas. Since for non weakly reversible networks, the associated directed graph G is connected but not strongly connected, G consists of r strongly connected components, which we denote by C 1 , C 2 , . . . , C r . Then, we can construct a direct acyclic graph G C , i.e. G C is a directed graph with no directed cycles as follows: a) G C has as nodes the r strongly connected components C 1 , C 2 , . . . , C r , b) for two nodes C i and C j of G C , if there exists a reaction C i ∋ S k a ℓk − − → S ℓ ∈ C j with a ℓk > 0, then there exists also the edge C i → C j on G C .
Due to the structure of G C , its nodes, or equivalently the strongly connected components of G, can be labeled as one of the following three types: where C 1 is a source component, C 2 is a transmission component and C 3 , C 4 are target components.
By definition, each of the three types of strongly connected components is subject to a different dynamic, which can be written as follows: Let C i be a strongly connected component and denote by X i the concentrations within C i . Moreover, denote by A i the reaction matrix formed by all reactions within the component C i . Then, we have • for a source component C i : are the in/out-flow of the transmission component C i . In the dynamics of transmission and target components, the in-flow F in i depends only on species which do not belong to C i , so that F in i can be treated as an external source for the system for C i . However, it may happen that F in i contains inflow from species whose behaviour is not a-priori known.
For acyclic graphs G C , however, it is possible to avoid these difficulties, since the topological order of acyclic graphs allows to re-order the r strongly connected components C 1 , C 2 , . . . , C r in such a way that for every edge C i → C j of G C it holds that i < j! This permits to study the dynamics of all components C i sequentially according to the topological order and, when at times considering a transmission component (or later a target component) C i , the required in-flow F in i contains only species whose behaviour is already known.
Due to the structure of the network, it is expected that species belonging to source or transmission components are subsequently losing mass such that the concentrations decay to zero in the large-time behaviour as time goes to infinity. In contrast, the species belonging to a target component converge to an equilibrium state, which is determined by the reactions within this component and by the mass "injected" from other components.
Since the source and transmission components do not converge to positive equilibria, the relative entropy method used for weakly reversible networks is directly not applicable. Instead, for each component C i we will modify the entropy method by introducing an artificial equilibrium state with normalised mass, which balances the reaction within C i . The artificial equilibrium will allow us to consider a quadratic functional, which is similar to the relative entropy in weakly reversible networks and which can be proved to decay exponentially to zero. This result is stated in the following Theorem: Theorem 1.5 (Exponential decay to zero of source and transmission components). For each C i being a source or a transmission component, there exist explicit constants C > 0 and λ i > 0 such that, for any specie S ℓ ∈ C i , the concentration u ℓ of S ℓ decays exponentially to zero, i.e.
For a target component C i , due to the in-flow F in i , the total mass of C i is not conserved but increasing. Hence, C i does not possesses an equilibrium as weakly reversible networks, which is explicitly given in terms of the reaction rates and the conserved initial total mass. However, since each target component is strongly connected and thus a weakly reversible reaction network with mass influx, there still exists a unique, positive equilibrium of C i denoted by X i,∞ , which balances the reactions within C i and has a total mass, which is the sum of the total initial mass of C i and the total "injected mass" from the other components via the in-flow F in i (see Lemma 3.4). We emphasis that in general the injected mass is not given explicitly but depends on the time evolution of all the influencing species higher up with respect to the topological order of the graph G C .
Since the equilibrium X i,∞ is positive, we can use again a relative entropy functional to prove the convergence of the species belonging to a target component to their corresponding equilibrium states.
Theorem 1.6 (Exponential convergence for target components). For each C i = {S i 1 , S i 2 , . . . , S i L } being a target component, there exists a unique positive equilibrium state X i,∞ = (u i 1 ,∞ , . . . , u i L ,∞ ) and for each specie S i ℓ belonging to C i , the concentration u i ℓ of S i ℓ converges exponentially to the corresponding equilibrium value for all t > 0, for some explicit constants C, λ i > 0.
Outline: The rest of the paper is organised as follows: In Section 2, we present the entropy method for weakly reversible networks and prove exponential convergence to the positive equilibrium. Non weakly reversible networks will be investigated in the Section 3. By using the structure of the underlying graphs, we are able completely resolve the large-time behaviour of all species belonging to such first order networks. Notation: We shall use the shortcut f = Ω f (x) dx, whenever |Ω| = 1, and · for the usual norm in L 2 (Ω), i.e.
By summing up all equations of system (2.1), the conservation of the total mass (1.3) follows immediately for solutions to (2.1) for all t > 0: where the initial mass M is assumed positive.

Definition 2.1 (Equilibrium).
A state X ∞ = (u 1,∞ , u 2,∞ , . . . , u N,∞ ) ∈ R N is called an equilibrium for (2.1) if X ∞ balances the reactions and satisfies the mass conservation law, that is (2.4) The existence of equilibria to (2.1) is nontrivial due to the singularity of the reaction matrix A and strongly related to the weak reversibility of the network. To be more precise, we will show that the system (2.1) admits a unique positive equilibrium if and only if the network is weakly reversible. In order to do that, we need to study the structure of the reaction matrix of the network, see e.g. [FPT].
Definition 2.2 (Decomposable reaction matrices). The reaction matrix A is called decomposable if there exists a permutation matrix P such that

5)
where B and D are square matrices. The reaction matrix A is called indecomposable if such a decomposition does not exist, i.e. if the dimension of either B or D is zero.
In the following Lemma 2.2 and 2.3, we state two criteria for a reaction matrix to be decomposable. The proof of Lemma 2.2 relies only on the graph structure while the proof of Lemma 2.3 relies on linear algebra. For the sake of readability we postpone the proofs to the Appendix 4.

Lemma 2.2 (Graph criterion for indecomposability).
A reaction matrix A is indecomposable if and only if the underlying graph G is strongly connected.
The following lemma provides an algebraic criterion to check the decomposability of a reaction matrix A. Lemma 2.3 (Algebraic criterion of decomposability, see [FPT]). For all i, j = 1, . . . , N, we denote by A ij ∈ R (N −1)×(N −1) the submatrix of A, where the i-th row and the j-th column are deleted from A. Moreover, we define the minor ρ ij := det(A ij ) for all i, j = 1, 2, . . . , N.
Then, a reaction matrix A is decomposable if and only if ρ ii = 0 for some i ∈ {1, 2, . . . , N}. Thus, for an indecomposable reaction matrix A, we have ρ ii = 0 for all i = 1, 2, . . . , N and moreover rank(A) = N − 1 and (1, . . . , 1) T is the only left-eigenvector to the eigenvalue zero, i.e. there is a unique conservation law for a weakly reversible reaction network.
Proof. The proof of the Lemma follows from direct calculations using the fact that all column sums of the reaction matrix A are zero.
Using the Lemmata 2.3 and 2.4, we are now able to show the existence of a unique positive equilibrium for (2.1).
Proposition 2.5 (Unique positive equilibria for weakly reversible networks).
The system (2.1) admits a unique positive equilibrium if and only if the network N is weakly reversible.
Proof. Sufficiency: Suppose that the network N is weakly reversible. Then, by Lemma 2.2, the reaction matrix A is indecomposable and the underlying graph G is strongly connected. Hence, from Lemma 2.3, we have that ρ ii = 0 for all i = 1, 2, . . . , N.
Then, by using the property that the sum of each column is zero, we can eliminate the last equation and rewrite the system as, by recalling that A N N is the (N − 1) × (N − 1) matrix obtained from A by removing the N-th row and N-th column and that det(A N N ) = ρ N N . Denote by D j , 1 ≤ j ≤ N − 1, the matrix obtained from A N N by replacing the j-th column by (−a 1N , −a 2N , . . . , −a N −1,N ) T . By the standard Cramer's rule, we have On the other hand, we have where we applied Lemma 2.4 at the last step. Then, from (2.6) we have u j,∞ = ρ jj ρ N N u N,∞ , for all j = 1, 2, . . . , N − 1. (2.7) From (2.7) and the mass conservation for all j = 1, 2, . . . , N thanks to sgn(ρ ii ) = (−1) N −1 for all i = 1, 2, . . . , N.
Necessity: Assume that N is not weakly reversible, then thanks to Lemma 2.2, the reaction matrix A is decomposable. By applying Lemma 2.3, there exists 1 ≤ j ≤ N such that ρ jj = 0. Similarly to the proof of the sufficient condition, we can write the equilibrium condition AX ∞ = 0 as Since det(A jj ) = ρ jj = 0, the system (2.8) either has no solution or infinitely many solutions. This contradicts the fact that the equilibrium X ∞ = (u 1,∞ , . . . , u N,∞ ) is unique. Therefore, N must be weakly reversible.
In the following, we will use the entropy method to study the trend to equilibrium. More precisely, for two trajectories X = (u 1 , u 2 , . . . , u N ) and Y = (v 1 , v 2 , . . . , v N ) to (2.1), where Y (t) has non-zero components for all times t > 0, we consider the following quadratic relative entropy functional (2.9) The following key Lemma 2.6 provides an explicit expression of the entropy dissipation associated to (2.9): Lemma 2.6 (Relative entropy dissipation functional). Assume that v i (t) = 0 for all i = 1, 2, . . . , N and t > 0. Then, we have Proof. For convenience we recall that a ij u j , for all i = 1, 2, . . . , N, a ij v j , for all i = 1, 2, . . . , N.
In order to simplify the following calculations, we introduce the difference to the equilibrium W := (w 1 , w 2 , . . . , w N ) = (u 1 − u 1,∞ , u 2 − u 2,∞ , . . . , u N − u N,∞ ) = X − X ∞ , and remark that thanks to the linearity of the system, the difference W is the solution to (2.1) subject to the shifted initial data for all x ∈ Ω.
Note that the total initial mass corresponding to W is zero, i.e.
and that W conserves the zero mass By using the relative entropy dissipation functional derived in Lemma 2.6, we have The following Lemma about entropy-entropy dissipation estimate is the main key to prove the convergence to equilibrium for (2.1).

Lemma 2.7 (Entropy-Entropy Dissipation Estimate).
There exists an explicit constant λ > 0 such that Proof. We divide the proof in several steps: Step 1. (Additivity of the relative entropy w.r.t. spatial averages) Straightforward calculation leads to where we denote W = (w 1 , w 2 , . . . , w N ) and we recall that w i = Ω w i dx for i = 1, . . . , N due to |Ω| = 1.
Step 2. (Entropy dissipation due to diffusion) By using Poincaré's inequality Step 3. (Entropy dissipation due to reactions) From (2.17) and (2.19), it remains to control By using Jensen's inequality we have (2.20) It then remains to prove that for some γ > 0. Note that if both reactions S i → S j and S j → S i do not appear in the reaction network, then we have a ij = a ji = 0 and thus a ij u j,∞ + a ji u i,∞ = 0.

Hence, the expression
may not contain all pairs (i, j) with i = j. However, the weak reversibility of the network allows to make all pairs (i, j) with i = j appear in the following sense: There exists a constant ξ > 0 such that (2.22) Indeed, assume that a ij = a ji = 0 for some i = j. Due to the weak reversibility of the network, there exists a path from S i to S j as follows a jr j r−1 − −−− → S jr ≡ S j with r ≥ 3 and a j k j k−1 > 0 for all k = 2, 3, . . . , r. Thus, there exists a σ > 0 such that Since there can only be finitely many pairs (i, j) with a ij = a ji = 0, we can repeat this procedure to finally get ( with the constraint of the conserved zero total mass Because of (2.25), (2.26) Therefore, we can estimate for various constants C N i,j=1;i<j In conclusion, we have proved (2.24), which in combination with (2.22) implies (2.21) and thus completes the proof of this Lemma.
Theorem 2.8 (Convergence to Equilibrium). Assume that the reaction network in Figure 1 is weakly reversible, the diffusion coefficients d i > 0 for all i = 1, 2, . . . , N, and the initial mass M > 0 is positive.
Then, the unique global solution to (2.1) converges to the unique positive equilibrium X ∞ in the following sense: where the constant λ > 0 can be computed explicitly.
Proof. From Lemma 2.7 we have d dt By Gronwall's inequality, and the proof is complete.
We now turn to the case of degenerate diffusion, where some of the diffusion coefficients d i can be zero. In the proof of the Theorem 2.8, we have used non-degenerate diffusion of all species in order to control distance of the concentrations to their spatial averages (see estimate (2.19)). This procedure must thus be adapted in the case of degenerate diffusion.
It was already proven in [DF07,FLT,MHM] that even if some diffusion coefficients vanish, one can still show exponential convergence to equilibrium provided reversible reactions. The technique used in these mentioned references is based on the fact that diffusion of one specie, which is connected through a reversible reaction with another specie, induces a indirect kind of "diffusion effect" to the latter specie.
We will prove that this principle is still valid for weakly reversible reaction networks as considered in this section.
Theorem 2.9 (Convergence to Equilibrium with Degenerate Diffusion). Assume that the reaction network in Figure 1 is weakly reversible, the initial mass is positive, and at least one diffusion coefficient d i is positive.
Then, the solution to (2.1) converges exponentially to equilibrium via the following estimate for some explicit rate λ ′ > 0.
Proof. We aim for a similar entropy-entropy dissipation inequality as stated in Lemma (2.7), i.e. we want to find a constant λ ′ > 0 such that (2.28) Due to the degenerate diffusion, the diffusion part of D(W |X ∞ ) is insufficent to control E(W − W |X ∞ ) as in (2.19), since some of diffusion coefficients can be zero. This difficulty can be resolved by quantifying the fact that diffusion of one specie is transferred to another species when connected via a weakly reversible reaction path. Without loss of generality, we assume that d 1 > 0 and estimate D(W |X ∞ ) by By arguments similar to (2.22) and (2.23), we have To control E(W −W |X ∞ ), we use the following estimate for all i = 2, 3, . . . , N: for some β > 0. Indeed, thanks to Poincaré's inequality ∇f 2 ≥ C P f − f 2 , we estimate for various sufficiently small constants C (2.32) Now, thanks to (2.30) and (2.31) (2.33) Thus (2.28) is proved and the proof is complete.
Remark 2.1. The estimate (2.31) is usually interpreted as follows: the sum of the dissipation due to the diffusion of w 1 and the dissipation caused by the reaction between w 1 and w i are bounded below by (2.32), which is essentially a diffusion dissipation term of the specie w i (after having applied Poincaré's inequality). In this sense, a "diffusion effect" has been transferred onto w i . We remark that while the presented proof for the linear case is straightforward, the proof of an analogous estimate to (2.31) in nonlinear cases turns out to be quite tricky. Readers are referred to [DF07] or [FLT,Lemma 3.6] for more details.

Non-weakly reversible networks
In this section, we consider reaction networks N which are not weakly reversible and, thus, correspond to a connected yet not strongly connected directed graph G. We will show that in the large time behaviour, each specie tends exponentially fast either to zero or to a positive equilibrium value depending on its position in the graph representing the network.
For weakly reversible reaction-diffusion networks (corresponding to strongly connected graphs), it was proven in Section 2 that each specie converges exponentially fast to a unique, positive equilibrium value, which is given explicitly in terms of the reaction rates and the conserved initial total mass.
For non weakly reversible reaction networks, however, we will show that while the equilibria are still unique and attained exponentially fast, the equilibrium values are in general no longer explicitly given but depend on the position in the graph in general and on the history of the concentrations of the influencing species in particular.
Moreover, since non weakly reversible reaction networks (2.1) may no longer have positive equilibria, the relative entropy method used in Section 2 is not directly applicable. Nevertheless, we will see that the relative entropy and the ideas of the entropy method still play the essential role our analysis of non-weakly reversible networks.
As the large time behaviour of the species depend on their position within the network, we need to first state some important properties of the underlying graph G. The following Lemmas 3.1 and 3.2 are well known in graph theory. We refer the reader to the book [BJG08] for a reference.
Lemma 3.1 (Strongly connected components form acyclic graphs G C ). Let G be a directed graph which is connected but not strongly connected such that the graph G contains at least r ≥ 2 strongly connect components, which we shall denote by C 1 , C 2 , . . . , C r . Thus, we can define a directed graph G C of strongly connected components as follows -G C has as nodes the r strongly connected components C 1 , C 2 , . . . , C r , -for two nodes C i and C j of G C , if there exists a reaction C i ∋ S k a ℓk − − → S ℓ ∈ C j with a ℓk > 0, then we define a directed edge C i → C j of G C . Then, the directed graph G C is acyclic, that is G C does not contain any cycles.
Proof. The proof can be found in e.g. [BJG08,Chapter 1] and shows that if G C would contain a cycle than this cycle should have been contained in a strongly connected component in the first place.
Lemma 3.2 (Topological order of acyclic graphs, [BJG08,Chapter 1]). There exists a reordering of the nodes of G C in such a way that for all direct edges C i → C j we always have i < j.
From now on, we will always consider topologically ordered graphs G C . For each i = 1, 2, . . . , N, we denote by N i the number of species belonging to Each component C i belongs to one of the following three types: there does not exist an edge C i ∋ S k → S j ∈ C i , • Transmission component: If C i is neither a source component nor a target component, then C i is called a transmission component.
The above classification of strongly connected components greatly improves the notation of the corresponding dynamics, which quantifies the behaviour of the species belonging to the three types of components. In the following, we denote by . . . , u L[i] ) T the concentration vector of the species belonging to C i . The evolution of the species belonging to a component C i depends on the type of C i : (i) For a source component C i , the system for X i is of the form where the diffusion matrix D i is (3.4) and the out flow matrix is defined as where the lower summation index L[i] + 1 follows for the topological order of the graph G C .
Roughly speaking, f L[i−1]+k is the sum of all the reaction rates from the specie S L[i−1]+k to species outside of C i . It may happen that f L[i−1]+k = 0 for some k = 1, 2, . . . , N, but there exists at least one k 0 such that f L[i−1]+k 0 > 0 since C i is a source component.
(ii) If C i is a transmission component, the system for X i writes as where the diffusion matrix D i , the reaction matrix A i and the out flow matrixF out i are defined as above in (3.3), (3.4) and (3.5), respectively. The in-flow vector F in i is defined by We remark that by studying all components C i within the topological order of G C , the dynamics of the previous components C 1 , C 2 , . . . , C i−1 is already known at the time we analyse the component C i . Thus, in system (3.6) the in-flow vector F in i can be considered as a given external in-flow.
x ∈ Ω, where the reaction matrix A i and the in-flow F in i are defined in the same way as above in (3.4) and (3.7).
By modifying the relative entropy method in Section 2, we obtain the first main result of this section: Theorem 3.3 (Exponential decay to zero for source and transmission components). Let C i be a source component or a transmission component.
Then, there exist constants C > 0 and λ i > 0 such that, for all S j ∈ C i , the concentration u j of S j decays exponentially to zero, i.e.
Proof. Since the ongoing outflow vanishes the mass of all source components and subsequently all transmission components, the corresponding equilibrium values are expected to be zero and the relative entropy method used for weakly reversible networks is not directly applicable here. We instead introduce a concept of "artificial equilibrium states with normalised mass" for these components, which allows to derive a quadratic entropy-like functional, which can be proved to decay exponentially. Due to their different dynamics, we have to distinguish the two cases: C i is a source component and C i is a transmission component.
In order to simplify the notation, we shall denote Then, the concentration vector X i and the reaction matrix A i can be rewritten as Case 1: C i is a source component.
We recall the corresponding system from (3.2) x ∈ Ω. (3.9) We now introduce an artificial equilibrium state X i,∞ = (v 1,∞ , v 2,∞ . . . , v N i ,∞ ) T with normalised mass to (3.9), which is defined as the solution of the system (3.10) It follows from Lemma 2.5 that there exists a unique positive solution X i,∞ to (3.10). Here we notice that X i,∞ balances all reactions within C i while the total mass contained in X i,∞ is normalised to one.
In the following we will study the evolution of the quadratic entropy-like functional By similar calculations as in Lemma 2.6, we obtain the time derivative of this quadratic functional (3.12) We remark that since C i is a source component, there exists an index k 0 ∈ {1, 2, . . . , N i } such that the out-flow f L[i−1]+k 0 > 0 is positive. Then, an estimate similar to (2.22) gives for various constants C for a constant λ i > 0. It follows that which proves the statement of the Theorem in the case that C i is a source component.
Case 2: C i is a transmission component. By recalling that the components C i are topologically ordered, we can assume without loss of generality that u ℓ , with ℓ = 1, 2, . . . , L[i − 1], obeys the following exponential decay u ℓ (t) 2 ≤ Ce −λ * t , ℓ = 1, 2, . . . , L[i − 1], for all t > 0. (3.14) for some λ * > 0. We also recall the system for C i , (3.15) with F in i is defined as (3.7). Denote by X i,∞ = (v 1,∞ , . . . , v N i ,∞ ) T the artificial equilibrium state of (3.15), which is the unique positive solution to (3.16) Again, we can compute the time derivative of to be estimated. Thanks to the decay (3.14) of u ℓ , we can estimate Then, with the help of (3.19), we estimate and similarly to (3.13), we obtain for some λ > 0, (3.20) From (3.20), we can use the classic Gronwall lemma to have which ends the proof in the case that C i is a transmission component.
For a target component, we need to define its corresponding equilibrium state. This equilibrium state balances the reactions within the component and has as total mass the sum of the initial total mass of the target component plus the total "injected mass" from the other components. In general, the injected mass will not be given explicitly but depend on the time evolution of the influences species prior to C i in terms of the topological order.
Lemma 3.4 (Equilibrium state of target components). For each target component C i , there exists a unique positive equilibrium state Proof. By Theorem 3.3 we have u ℓ (t) 2 ≤ Ce −λ * t for some λ * > 0. Thus, Jensen's inequality yields and the right hand side of the second equation in (3.21) is finite. Therefore, the existence of a unique X i,∞ satisfying (3.21) follows from Proposition 2.5.
Theorem 3.5 (Convergence to equilibrium states for target components). Assume that C i is a target component and X i,∞ is the equilibrium state of C i defined by Lemma 3.4. Then, there exist C > 0 and λ i > 0 such that for k = 1, 2, . . . , N i , we have Proof. We recall first the system for a target component C i , a L[i−1]+ℓ,k u k .
Note that the total mass of C i is not conserved but increases in time due to the in-flow vector F in i . To compute the total mass of C i at a time t > 0, we sum up all the equations of (3.23) then integrating over Ω, thanks to the homogeneous Neumann boundary condition and the fact that (1, . . . , 1) T is a left eigenvector with eigenvalue zero of A i since A i is a reaction matrix. Thus, we have ( 3.24) We define v k (t) = u L[i−1]+k − v k,∞ the distance from u L[i−1]+k to its corresponding equilibrium state for all k = 1, 2, . . . , N i . It implies that (v k ) k=1,...,N i solves the system (3.23) subject to the initial data v k,0 = u L[i−1]+k,0 − v k,∞ for all k = 1, 2, . . . , N i . We define V i = (v 1 , v 2 , . . . , v N i ) and consider the relative entropy-like functional (3.25) By using again arguments of Lemma 2.6, we calculate the entropy dissipation where we have used u ℓ (t) 2 ≤ Ce −λ * t for some λ * > 0 and for all ℓ = 1, . . . , L[i − 1] in the last estimate. By inserting (3.27) into (3.26), we obtain (3.28) It follows from (3.28) that (3.29) To control E 2 , we use arguments similar to Step 3 in the proof of Lemma 2.7. First, by using (3.24), we have the total mass of (v k ) 1≤k≤N i is computed as, (3.31) By using (2.22) and (3.31), we estimate for ε > 0 is sufficiently small. It follows from (3.30) that where we have used the fact that u ℓ (t) 2 ≤ Ce −λ * t for all ℓ = 1, 2, . . . , L[i−1]. Hence, (3.32) implies that (3.33) Combining (3.33) and (3.29) yields (3.34) Therefore, by applying a classic Gronwall lemma, Finally, by recalling (3.25), we obtain where λ i = min{c, λ * }. This completes the proof of the Theorem.

Appendix
In this Appendix, we give the proofs of the Lemmas 2.2 and 2.3.
Proof of Lemma 2.2. Necessity: We prove by contradiction. Assume that the graph G is not strongly connected. Then, it follows that we can divide the set of vertices {S i } N i=1 into two sets C 1 = {S i 1 , S i 2 , . . . , S i k } and C 2 = {S i k+1 , S i k+2 , . . . , S i N } with 1 ≤ k ≤ N − 1 such that there is no edge starting at a vertex in C 2 and finishing at a vertex in C 1 . Let σ be the permutation satisfying σ(i j ) = j for all j = 1, 2, . . . , N. Note that σ is a re-numbering of the vertices of the graph. We denote by G σ the graph after re-numbering by σ and by A σ the corresponding reaction matrix of G σ . Then, it is obvious that A σ has the following block-diagonal form where B σ ∈ R k×k and D σ ∈ R (N −k)×(N −k) . On the other hand, we have the following relation A σ = P σ AP T σ (4.2) between A and A σ , where P σ is the permutation matrix corresponding to σ. Thus (4.1) and (4.2) imply that A is decomposable, which leads to a contradiction.
Sufficiency: Assume that the matrix A is decomposable, i.e. there is a permutation matrix P such that with B ∈ R k×k and D ∈ R (N −k)×(N −k) . Then, by re-numbering the vertices of G by the permutation σ corresponding to P , we obtain a graph G σ which has the following property: There is no edge starting at a vertex in {k + 1, k + 2, . . . , N} and finishing at a vertex in {1, 2, . . . , k}. This means that G σ is not strongly connected and so G is also not strongly connected.
To prove Lemma 2.3, we need the following elementary Lemma.