New universal Lyapunov functions for non-linear reaction networks

TL;DR

This paper introduces a new broad family of universal Lyapunov functions for nonlinear reaction networks, extending classical concepts and depending only on equilibrium states and network structure, not on specific kinetic parameters.

Contribution

The authors constructed a new family of universal Lyapunov functions applicable to any nonlinear reaction network with detailed or complex balance, independent of kinetic parameters.

Findings

New Lyapunov functions depend only on equilibrium and network structure.

Applicable to all nonlinear reaction networks with detailed or complex balance.

Constructed using partial equilibria and convex envelopes.

Abstract

In 1961, R\'enyi discovered a rich family of non-classical Lyapunov functions for kinetics of the Markov chains, or, what is the same, for the linear kinetic equations. This family was parameterised by convex functions on the positive semi-axis. After works of Csisz\'ar and Morimoto, these functions became widely known as $f$-divergences or the Csisz\'ar--Morimoto divergences. These Lyapunov functions are universal in the following sense: they depend only on the state of equilibrium, not on the kinetic parameters themselves. Despite many years of research, no such wide family of universal Lyapunov functions has been found for nonlinear reaction networks. For general non-linear networks with detailed or complex balance, the classical thermodynamics potentials remain the only universal Lyapunov functions. We constructed a rich family of new universal Lyapunov functions for {\em any non-linear reaction network} with detailed or complex balance. These functions are parameterised by compact subsets of the projective space. They are universal in the same sense: they depend only on the state of equilibrium and on the network structure, but not on the kinetic parameters themselves. The main elements and operations in the construction of the new Lyapunov functions are partial equilibria of reactions and convex envelopes of families of functions.

Equations

$$H(f)=\int f(x,v)\ln (f(x,v)){\mathrm{d}}^{3}v{\mathrm{d}}^{3}x$$

$$H=\sum _{i=1}^n{c}_i\left(\ln \left(\frac{c_i}{c_i^{\mathrm{eq}}}\right)-1\right),$$

$$H_f(c)=\sum _i{c}_i^{\mathrm{eq}}f\left(\frac{c_i}{c_i^{\mathrm{eq}}}\right),$$

$$\sum _{i=1}^n{\alpha }_{ri}{A}_i\to \sum _{j=1}^n{\beta }_{rj}{A}_j(r=1,\dots ,m),$$

$$w_r=k_r\prod _{i=1}^n{c}_i^{{\alpha }_{ri}},$$

$$\frac{\mathrm{d}c}{\mathrm{d}t}=\sum _{r=1}^m{\gamma }_r{w}_r$$

$$\frac{\mathrm{d}c}{\mathrm{d}t}\in Span\{{\gamma }_r|r=1,\dots ,m\}.$$

$$\text{mbox}Foralli,r\text{mbox}either{\alpha }_{ri}=0\text{mbox}or{\alpha }_{ri}\ge 1,$$

$$(Span\{{\gamma }_r|r=1,\dots ,m\}+c)\cap {\mathbb{R}}_{\ge 0}^n,$$

$$\sum _{i=1}^n{\alpha }_{ri}{A}_i\rightleftharpoons \sum _{j=1}^n{\beta }_{rj}{A}_j(r=1,\dots ,m),$$

$$w_r^{+}=k_r^{+}\prod _{i=1}^n{c}_i^{{\alpha }_{ri}},w_r^{-}=k_r^{-}\prod _{i=1}^n{c}_i^{{\beta }_{ri}},$$

$$w_r^{+}(c_i^{\mathrm{eq}})=w_r^{-}(c_i^{\mathrm{eq}})\text{mbox}forallr.$$

$$k_r^{+}\prod _{i=1}^n(c_i^{\mathrm{eq}}{)}^{{\alpha }_{ri}}=k_r^{-}\prod _{i=1}^n(c_i^{\mathrm{eq}}{)}^{{\beta }_{ri}}=w_r^{\mathrm{eq}}.$$

$$k_r^{+}=w_r^{\mathrm{eq}}{\left(\prod _{i=1}^n(c_i^{\mathrm{eq}}{)}^{{\alpha }_{ri}}\right)}^{-1},k_r^{-}=w_r^{\mathrm{eq}}{\left(\prod _{i=1}^n(c_i^{\mathrm{eq}}{)}^{{\beta }_{ri}}\right)}^{-1}.$$

$$w_r^{+}=w_r^{\mathrm{eq}}\prod _{i=1}^n{\left(\frac{c_i}{c_i^{\mathrm{eq}}}\right)}^{{\alpha }_{ri}}{w}_r^{-}=w_r^{\mathrm{eq}}\prod _{i=1}^n{\left(\frac{c_i}{c_i^{\mathrm{eq}}}\right)}^{{\beta }_{ri}},$$

$$\frac{\mathrm{d}c}{\mathrm{d}t}=\sum _{r=1}^m{\gamma }_r{w}_r^{\mathrm{eq}}\left[\prod _{i=1}^n{\left(\frac{c_i}{c_i^{\mathrm{eq}}}\right)}^{{\alpha }_{ri}}-\prod _{i=1}^n{\left(\frac{c_i}{c_i^{\mathrm{eq}}}\right)}^{{\beta }_{ri}}\right].$$

$$\frac{\partial H}{\partial c_i}=\ln \left(\frac{c_i}{c_i^{\mathrm{eq}}}\right)$$

$$\frac{\mathrm{d}H}{\mathrm{d}t}=-\sum _{r=1}^m(w_r^{+}-w_r^{-})(\ln w_r^{+}-\ln w_r^{-})\le 0,$$

$$k_r=w_r^{\mathrm{eq}}{\left(\prod _{i=1}^n(c_i^{\mathrm{eq}}{)}^{{\alpha }_{ri}}\right)}^{-1}$$

$$\frac{\mathrm{d}c}{\mathrm{d}t}=\sum _{r=1}^m{\gamma }_r{w}_r^{\mathrm{eq}}\prod _{i=1}^n{\left(\frac{c_i}{c_i^{\mathrm{eq}}}\right)}^{{\alpha }_{ri}}.$$

$$\frac{\mathrm{d}H}{\mathrm{d}t}=\sum _i\frac{\partial H}{\partial c_i}\frac{\mathrm{d}{c}_i}{\mathrm{d}t}=\sum _r({\gamma }_r,\nabla H)w_r^{\mathrm{eq}}\exp ({\alpha }_r,\nabla H),$$

$$\theta (\lambda )=\sum _{\rho }{w}_r^{\mathrm{eq}}\exp [(\nabla H,(\lambda {\alpha }_r+(1-\lambda ){\beta }_r))]$$

$$\frac{\mathrm{d}H}{\mathrm{d}t}=-{\frac{\mathrm{d}\theta (\lambda )}{\mathrm{d}\lambda }|}_{\lambda =1}$$

$$\sum _r{w}_r^{\mathrm{eq}}\prod _{i=1}^n{\left(\frac{c_i}{c_i^{\mathrm{eq}}}\right)}^{{\alpha }_{ri}}=\sum _r{w}_r^{\mathrm{eq}}\prod _{i=1}^n{\left(\frac{c_i}{c_i^{\mathrm{eq}}}\right)}^{{\beta }_{ri}}.$$

$$R_j^{+}=\{r|{\alpha }_r=y_j\},R_j^{-}=\{r|{\beta }_r=y_j\}$$

$$\sum _{j=1}^q\prod _{i=1}^n{\left(\frac{c_i}{c_i^{\mathrm{eq}}}\right)}^{y_{ji}}\left[\sum _{r\in R_j^{+}}{w}_r^{\mathrm{eq}}-\sum _{r\in R_j^{-}}{w}_r^{\mathrm{eq}}\right]=0$$

$$\prod _{i=1}^n{\left(\frac{c_i}{c_i^{\mathrm{eq}}}\right)}^{y_{ji}}$$

$$\sum _{r\in R_j^{+}}{w}_r^{\mathrm{eq}}-\sum _{r\in R_j^{-}}{w}_r^{\mathrm{eq}}=0$$

$$\frac{\mathrm{d}c}{\mathrm{d}t}={\gamma }_r\left[\prod _{i=1}^n{\left(\frac{c_i}{c_i^{\mathrm{eq}}}\right)}^{{\alpha }_{ri}}-\prod _{i=1}^n{\left(\frac{c_i}{c_i^{\mathrm{eq}}}\right)}^{{\beta }_{ri}}\right].$$

$$\begin{array}{rl}\frac{\mathrm{d}G(c,c^{\mathrm{eq}})}{\mathrm{d}t}& =\sum _{r=1}^m{w}_r^{\mathrm{eq}}({\nabla }_{c}G(c,c^{\mathrm{eq}}),{\gamma }_r)\left[\prod _{i=1}^n{\left(\frac{c_i}{c_i^{\mathrm{eq}}}\right)}^{{\alpha }_{ri}}-\prod _{i=1}^n{\left(\frac{c_i}{c_i^{\mathrm{eq}}}\right)}^{{\beta }_{ri}}\right]\\ & =\sum _{r=1}^m{w}_r^{\mathrm{eq}}{D}_{r}G(c,c^{\mathrm{eq}}),\end{array}$$

Full text

Despite many years of research, no such wide family of universal Lyapunov functions has been found for nonlinear reaction networks. For general non-linear networks with detailed or complex balance, the classical thermodynamics potentials remain the only universal Lyapunov functions.

We constructed a rich family of new universal Lyapunov functions for any non-linear reaction network with detailed or complex balance. These functions are parameterized by compact subsets of the projective space. They are universal in the same sense: they depend only on the state of equilibrium and on the network structure, but not on the kinetic parameters themselves.

The main elements and operations in the construction of the new Lyapunov functions are partial equilibria of reactions and convex envelopes of families of functions.

1 Introduction

The classical Lyapunov functions in kinetics are closely related to the concepts of entropy and free energy. The first example was provided by Boltzmann in 1872 [1]. He proved that in an isolated system the functional

[TABLE]

does not increase in time due to kinetic equation (the Boltzmann equation). Here, $f(x,v)$ is the distribution density of particles in the phase space, $x$ is position in space, $v$ is velocity of a particle. Boltzmann’s proof [1] used the principle of detailed balance. Several years later he found more general conditions for $H$-theorem [2] and invented what we call now semidetailed balance or cyclic balance or complex balance (for modern theory of chemical reaction networks we refer to tutorial [3]).

For perfect chemical mixtures with components $A_{1},\ldots,A_{n}$ in isothermal isochoric conditions (fixed volume) the analogue of Boltzmann’s $H$-function is:

[TABLE]

where $c_{i}\geq 0$ is the concentration of $A_{i}$ and $c_{i}^{\rm eq}>0$ is an equilibrium concentration of $A_{i}$ (under the standard convention that $x\ln x=0$ for $x=0$). The idealization of perfect mixtures is applicable to rarefied gases or to reactions of small admixtures in solutions.

The $H$-function (2) is a Lyapunov function for all mass action law systems with detailed or complex balanced equilibrium (see e.g. [3])). We call this property ‘universality’.

In information theory, the function $H$ appears as a measure of relative information (in the distribution $c_{i}$ with respect to the distribution $c_{i}^{\rm eq}$) and analogue of the $H$-theorem states that random manipulations with data decrease the relative information with respect to the equilibrium that does not change under manipulations [4, 5, 6].

It is not much surprising that the $H$ function (2) is essentially the only universal Lyapunov function for all imaginable perfect kinetic systems with detailed balance. Nevertheless, if we restrict the choice of the reaction mechanism then the class of Lyapunov functions, which are conditionally independent of reaction rate constants for a given detailed balanced or complex balanced equilibrium, can be extended. We call such Lyapunov functions conditionally universal (for a given reaction mechanism).

In 1961, Rényi discovered a class of conditionally universal Lyapunov functions for Markov chains [7]. After works [8, 9] these functions were studied by many authors under the name $f$-divergences or Csiszár–Morimoto divergences. It is known that any universal Lyapunov functions for Markov chains has the form of $f$-divergence [10, 11, 12] or is a monotonic function of such a divergence. The continuous time Markov kinetic equation coincide with the kinetic equations for linear (monomolecular) reactions of perfect systems, and $f$-divergences can be considered as a direct generalization of (2):

[TABLE]

where $f$ is a convex function on the positive semi-axis.

For $f(x)=x\ln x$ and after adding a constant term proportional to $\sum_{i}c_{i}$ we get the classical formula (2). (Recall that $\sum_{i}c_{i}$ does not change in linear kinetics.)

Existence of a very rich family of conditionally universal Lyapunov functions for the linear reaction mechanisms makes us guess that there should be many conditionally universal Lyapunov functions for any given nonlinear reaction mechanism as well. In this paper, we construct new conditionally universal Lyapunov functions for any given reaction mechanism using partial equilibria of all single reactions and detailed or complex balance conditions.

In the next Sec. 2, we give the necessary formal definitions and introduce notations for mass action law systems. The necessary and sufficient conditions that a convex function is a Lyapunov function for all reaction networks with given reaction mechanism and equilibrium point under detailed or complex balance assumption are proven in Sec.3. The construction of a new family of conditionally universal Lyapunov functions for any reaction network is presented and the main result, Theorem 2, is proven for mass action law systems in Sec. 4. In Sec. 5 we outline the possible generalizations and applications of the results. In Conclusion the main results of the work are summarized and an open question is formulated.

2 Prerequisites: mass action law and classical Lyapunov functions

In this section, we formally introduce mass action law and equations of chemical kinetics. For more detailed introduction, including thermodynamical backgrounds, detailed kinetics, applied kinetics, and mathematical aspect of kinetics, we refer to the modern book [13]. Tutorial [3] gives the mathematical introduction in dynamics of chemical reaction networks. Formalism of chemical kinetics with special attention to heterogeneous catalysis is discussed in detail in the monograph [14].

2.1 Mass action law

Consider a closed system with $n$ chemical species $A_{1},\ldots,A_{n}$, participating in a complex reaction network. The reaction network is represented in the form of the system of stoichiometric equations of elementary reactions (called also reaction mechanism):

[TABLE]

where $\alpha_{ri}\geq 0$, $\beta_{rj}\geq 0$ are the stoichiometric coefficients, $r=1,\ldots,m$, $i,j=1,\ldots,n$, $m$ is the number of elementary reactions, $n$ is the number of components. In this representation, the direct and reverse elementary reactions are considered separately.

The stoichiometric vector $\gamma_{r}$ of the elementary reaction is $\gamma_{r}=(\gamma_{ri})$, $\gamma_{ri}=\beta_{ri}-\alpha_{ri}$ (‘gain minus loss’). The gain vector is $\alpha_{r}=(\alpha_{ri})$ and the loss vector is $\beta_{r}=\beta_{ri}$.

Elementary reactions of the form $A_{i}\to A_{j}$ are called linear or monomolecular reactions.

According to the mass action law, the reaction rate for the elementary reactions (4) are

[TABLE]

where $k_{r}\geq 0$ is the reaction rate constant and the standard convention is used: for any $x\geq 0$, $x^{0}=1$.

The kinetic equations for a perfect system in isochoric isothermal conditions have the form

[TABLE]

The stoichiometric subspace is $Span\{\gamma_{r}|r=1,\ldots,m\}$. For any set of values $k_{r}\geq 0$ and $c_{i}\geq 0$ the time derivative of $c$ belongs to the stoichiometric subspace

[TABLE]

Therefore, each linear functional $b(c)$ that annuls $Span\{\gamma_{r}|r=1,\ldots,m\}$ (i.e. $b(\gamma_{r})=0$ for all $r$ is the conservation law: ${\mathrm{d}}b(c)/{\mathrm{d}}t=0$ according to (6). Such functionals are called stoichiometric conservation law.

We always assume that there exists a strictly positive stoichiometric conservation law that is such a vector $b=(b_{i})$, $b_{i}>0$ that $\sum_{i}b_{i}\gamma_{ri}=0$ for all $r$. In chemical kinetics, this may be the conservation of mass or of total number of atoms, for example. Due to this assumption, every stoichiometric vector $\gamma$ has both positive and negative components. This assumption has many important consequences. For example, if

[TABLE]

then the right hand side of (6) is a Lipschitz function in any reaction polyhedron

[TABLE]

where $c$ is an arbitrary vector with non-negative coordinates and $\mathbb{R}^{n}_{\geq 0}$ is the cone (orthant) of such vectors. Therefore, solution of (6) exists and is unique for every initial state with non-negative coordinates. Moreover, it remains non-negative and, therefore, any reaction polyhedron (8) is positively invariant: solutions do not leave it in positive time. Algorithms for detailed analysis of the structure of reaction polyhedra are developed in [15]. It involves the classical double-description problem [16, 17]. Reaction polyhedra for various examples of chemical reactions are presented in [18, 19, 20]. The boundedness of positive semi-trajectories has strong stability implications for reaction networks [21].

In many applications, the stoichiometric coefficients are non-negative integers. For global existence, uniqueness of solutions, and preservation of positivity the condition (7) and existence of positive stoichiometric conservation law are sufficient. Moreover, in this work we do not need uniqueness and global existence of solutions and analyze just the local conditions ${\mathrm{d}}G/{\mathrm{d}}t\leq 0$ for various functions $G$ at strictly positive points $c$. Therefore, a weaker condition is needed: instead of (7) we assume below that $\alpha_{ri}\geq 0$, $\beta_{ri}\geq 0$.

2.2 Detailed balance and classical Lyapunov function

First proof of Boltzmann’s $H$-theorem [1] and its analogues for chemical kinetics were based on detailed balance . Boltzmann’s argument were analyzed by Tolman [22]. In general situation, this principle means that at equilibrium, each elementary process is equilibrated with its reverse process.

To formulate this principle for chemical kinetics, we combine direct and reverse reactions together and rewrite the reaction mechanism in the form

[TABLE]

(formally, the transformation from (4) to (9) is always possible if we allow some of $k_{r}^{-}=0$ when there is no reverse reaction in the initial mechanism).

We use notation $k_{r}^{+}$ and $k_{r}^{-}$ for reaction rate constants of direct and reverse reactions, and $w_{r}^{+}$, $w_{r}^{-}$ for reaction rates of these reactions:

[TABLE]

For the total rate of the couple of direct and reverse reaction we use $w_{r}=w_{r}^{+}-w_{r}^{-}$. In these notations the kinetic equations are the same (6) (with different value of $m$).

The principle of detailed balance for the mass action law is: there exists a positive equilibrium $c_{i}^{\rm eq}>0$ with detailed balance,

[TABLE]

If the set of the stoichiometric vectors $\{\gamma_{r}\}$ is linearly dependent then this condition implies algebraic relations between reaction rate constants. Each elementary reaction is equilibrated at the point of detailed balance, $c^{\rm eq}$. For systems that obey the mass action law, this means that

[TABLE]

It is convenient to use the detailed balance relations (12) and introduce a set of independent parameters instead of the reaction rate constants: the equilibrium fluxes $w_{r}^{\rm eq}$ and the equilibrium concentrations $c_{i}^{\rm eq}$. The reaction rate constants have a simple and explicit expression through the equilibrium flows and equilibrium concentrations:

[TABLE]

In this parameterization, the mass action law (10) with detailed balance condition takes the form:

[TABLE]

In this parameterization, the kinetic equations for the mass action law with detailed balance are

[TABLE]

The classical $H$-theorem for isochoric isothermal perfect chemical systems with detailed balance can be produced now by simple straightforward calculations. Notice that

[TABLE]

and calculate the time derivative of $H$ according to the system of kinetic equations with detailed balance condition (15):

[TABLE]

and ${{\mathrm{d}}H}/{{\mathrm{d}}t}=0$ for points of detailed balance an only for such points. Therefore, existence of a positive point of detailed balance for the mass action law kinetic equation implies that $H$ (2) does not increase in time and all the non-negative steady states of kinetic equations are detailed balance equilibria.

Definition of $H$ and $c_{i}^{\rm eq}$ requires additional comments. Equilibrium concentrations of chemical mixture depend on the conserved quantities those do not change in the course of chemical reactions. For example, they depend on the atomic balances. Without fixing these values, the vector of positive equilibria $c_{i}^{\rm eq}$ ($i=1,\ldots,n$) is defined ambiguously. Any vector of positive equilibria $c_{i}^{\rm eq}>0$ can be used in the definition of $H$ (2). After that, the $H$-function (2) can be used for all values of $c_{i}\geq 0$ for all possible values of conserved quantities. It is a simple exercise to show that for the mass action law kinetics with detailed balance the difference between $H$-functions with different choices of equilibria $c_{i}^{\rm eq}$ does not change in time.

The $H$-function (2) was utilised by Zeldovich in his proof of uniqueness of positive equilibrium for given values of conserved quantities (1938, reprinted in 1996 [23]). It was recognised as the main instrument for analysis of stability of perfect kinetic systems in 1960s-1970s [24, 25].

2.3 Complex balance and Lyapunov functions

In 1887, Lorentz stated that the collisions of polyatomic molecules are irreversible and, therefore, Boltzmann’s $H$-theorem is not applicable to the polyatomic media [26]. Boltzmann found the solution immediately and invented what we call now semidetailed balance or cyclic balance or complex balance [2].

Now, it is proven that the Lorentz objections were wrong and the detailed balance conditions hold for polyatomic molecules [27]. Nevertheless, this discussion was seminal. The complex balance is a popular assumption in chemical kinetics beyond the detailed balance [28, 29]. The comparative analysis of detailed and complex balance assumption in practice of modeling of chemical reaction networks was presented in work [30]. It is demonstrated how the generalized mass action law with complex balance appears as a macroscopic limit of the microscopic Markov kinetics [31]. The formal structures of complex balance are also useful for analysis of systems with time delays [32].

Return to the ‘irreversible’ representation of the reaction mechanism (4), where direct and reverse reaction are considered separately. Let $c^{\rm eq}$ be a positive vector of concentrations (it will be a point of complex balance). Similarly to (13) introduce new parameters

[TABLE]

$w_{r}^{\rm eq}$ is the rate of the $r$th reaction at the point $c^{\rm eq}$. The kinetic equations with these parameters are

[TABLE]

Let us calculate the time derivative of the function $H$ (2) by virtue of the system (19).

[TABLE]

where $\nabla H$ is given by (16)

An auxiliary function $\theta(\lambda)$ of one variable $\lambda\in[0,1]$ is convenient for analysis of ${\mathrm{d}}H/{\mathrm{d}}t$ (see [15, 31]):

[TABLE]

With this function, ${{\mathrm{d}}H}/{{\mathrm{d}}t}$ defined by (20) has a very simple form:

[TABLE]

The condition $\theta(1)=\theta(0)$ is sufficient for the inequality $\theta^{\prime}(1)\leq 0$, because $\theta(\lambda)$ is a convex function. Hence, if $\theta(1)=\theta(0)$ then ${{\mathrm{d}}H}/{{\mathrm{d}}t}\leq 0$. The explicit form of the condition $\theta(1)=\theta(0)$ is

[TABLE]

Let us consider the family of all gain and loss vectors $\{\alpha_{r},\beta_{r}\}$ ($r=1,\ldots,m$). Usually, some of these $2m$ vectors coincide. Assume that there are $q$ different vectors among them. Let $y_{1},\ldots,y_{q}$ be these vectors. For each $j=1,\ldots,q$ we take

[TABLE]

We can rewrite (23) in the form

[TABLE]

The monomials

[TABLE]

are linearly independent functions of $c$ for every finite set of vectors $y_{j}$. Therefore for any $j=1,\ldots,q$

[TABLE]

This complex balance condition is equivalent to the condition $\theta(0)=\theta(1)$ and is sufficient for the inequality ${{\mathrm{d}}H}/{{\mathrm{d}}t}\leq 0$. Therefore, $H$ is a Lyapunov function for systems with complex balance.

3 General $H$ theorem for perfect systems with detailed and complex balance

The systems with detailed or complex balance have a classical Lyapunov function (2) but we are interested in construction of new Lyapunov functions. In this section, we find necessary and sufficient conditions that a convex function is a conditionally universal Lyapunov function for the reaction network with given reaction mechanism and equilibrium point under detailed or complex balance assumption.

Let us start from the systems with detailed balance. By definition, a function $G(c,c^{\rm eq})$ is a conditionally universal Lyapunov function for this reaction network if it is a Lyapunov function for system (15) for every set of non-negative values of the equilibrium fluxes $w_{r}^{\rm eq}$.

The following lemma is a simple consequence of the form of the kinetic equations (15). Let $G(c,c^{\rm eq})$ be continuous function and a convex function of $c$ for all $c^{\rm eq}$.

Lemma 1**.**

Time derivative of a function $G(c,c^{\rm eq})$ by virtue of system (15) is non-negative, ${\mathrm{d}}G(c,c^{\rm eq})/{\mathrm{d}}t\leq 0$, for all positive values of $w_{r}^{\rm eq}$, $c$, and $c^{\rm eq}$ if and only if the derivatives of this function by virtue of the following systems are non-negative for all $r$, all positive $c$ and given $c^{\rm eq}$:

[TABLE]

Proof.

For smooth $G(c,c^{\rm eq})$, calculate the time derivative of this function by virtue of the system (15):

[TABLE]

where $(\nabla_{c}G(c,c^{\rm eq}),\gamma_{r})$ is the standard inner product (derivative of $G(c,c^{\rm eq})$ in the direction $\gamma_{r}$) and $D_{r}G(c,c^{\rm eq})$ is the time derivative of $G(c,c^{\rm eq})$ by virtue of the system (27).

The coefficients $w_{r}^{\rm eq}$ are independent non-negative variables. Therefore, the time derivative of $G(c,c^{\rm eq})$ by virtue of the original system (15) is a conical combination of $D_{r}G(c,c^{\rm eq})$. Non-negativity of all conical combinations of $D_{r}G(c,c^{\rm eq})$ means that each term is non-negative. Thus, non-negativity of ${{\mathrm{d}}G(c,c^{\rm eq})}/{{\mathrm{d}}t}$ for the reaction network with any values of equilibrium fluxes is equivalent to non-negativity of ${{\mathrm{d}}G(c,c^{\rm eq})}/{{\mathrm{d}}t}$ for each single reaction subsystem of the network (i.e. to the inequality $D_{r}G(c,c^{\rm eq})\geq 0$ for all $r$). For continuous convex functions $G(c,c^{\rm eq})$ we have to use in (32) the subgradients of $G(c,c^{\rm eq})$ instead of the gradients. (Recall that according to A.D. Alexandrov theorem continuous convex functions are almost everywhere twice differentiable and, therefore, the criterion (32) with classical gradients is valid for them almost everywhere.) The time derivative in (28) becomes an interval of non-positive numbers and the proof remains the same. For more detail about continuous but non-smooth Lyapunov functions we refer to [33, 34] ∎

Lemma 1 allows us to reduce a complex validation of inequalities with $m+2n$ variables $w_{r}^{\rm eq}$, $c_{i}^{\rm eq}$, and $c_{i}$ to a series of $m$ simpler inequalities with $2n$ variables. Moreover, non-negativity of $D_{r}G(c,c^{\rm eq})$ means that the function $G(c,c^{\rm eq})$ does not increase with time along solutions of one-dimensional kinetic equations (27): in this system, ${{\mathrm{d}}c}/{{\mathrm{d}}t}$ is proportional to vector $\gamma_{r}$ and for any positive solution $c(t)$ ($t>0$) the difference $c(t_{1})-c(t_{0})$ is always proportional to $\gamma_{r}$, $c(t_{1})-c(t_{0})=\xi\gamma_{r}$ with some scalar multiplier $\xi$ ($t_{0,1}>0$). This one-dimensional dynamics admits a strongly convex Lyapunov functions $H$ (2) with time derivative $-(w_{r}^{+}-w_{r}^{-})(\ln w_{r}^{+}-\ln w_{r}^{-})$ (17).

For each positive concentration vector $c$, the interval

[TABLE]

is positively invariant with respect to (27). Restriction of $H$ (2) on this interval is a strongly convex function. The derivative of this function has the logarithmic singularity at the ends of the interval. $H$ has the unique minimizer on $I$. It is the unique positive equilibrium point of (27) on $I$ (and, by definition, a partial equilibrium of the complete system (15)). These observations allow us to formulate the following criterion for the conditionally universal convex Lyapunov functions of the reaction kinetics (15) for reaction networks with detailed balance.

Let $G(c,c^{\rm eq})$ be continuous function and a convex function of $c\in\mathbb{R}_{>0}^{n}$ for all $c^{\rm eq}$.

Theorem 1**.**

$G(c,c^{\rm eq})$* is a monotonically non-increasing function of time on the positive solutions of the kinetic equations (15) for all non-negative values of equilibrium fluxes if and only if for every positive concentration vector $c$ and every stoichiometric vector $\gamma_{r}$ of the reaction mechanism the minimizer of $H(c,c^{\rm eq})$ on the interval $I$ (29) is, at the same time, a minimizer of $G(c,c^{\rm eq})$ on this interval:*

[TABLE]

Proof.

According to Lemma 1, it is necessary and sufficient to prove this theorem for the one-step reaction ($m=1$). Consider restriction of the one-step kinetic equation (27) on the interval $I$ (29). On this interval, the system has one equilibrium (the partial equilibrium $c_{\gamma_{r}}^{*}$). It is stable, and the restriction of $H$ on this interval is the Lyapunov function of the system. The point $c_{\gamma_{r}}^{*}$ is the minimizer of $H$ on $I$.

Assume that $c_{\gamma_{r}}^{*}$ is a minimizer of $G(c,c^{\rm eq})$ for $c\in I$. Then convexity of $G(c)$ implies that $G$ is monotonically non-increasing function of time on the solution $c(t)$ of (27) on $I$. ($G(c(t))$ decreases monotonically to the minimal value when the point $c(t)$ approaches its minimizer $c_{\gamma_{r}}^{*}$.)

Assume now that $G(c(t))$ does not increase in time due to dynamic of (27) on $I$. This dynamics lead $c(t)$ to the unique equilibrium $c_{\gamma_{r}}^{*}$. This equilibrium should be a minimizer of $G$ on $I$. Indeed, if at some point $y\in I$ the function $G$ takes smaller value then in $c_{\gamma_{r}}^{*}$, then in the motion from $y$ to $c_{\gamma_{r}}^{*}$ the value of $G$ should increase, which contradicts the assumption. ∎

Thus, to check that a convex function is a conditionally universal function for the reaction network with detailed balance, it is sufficient to check that its minimizers in the direction of the stoichiometric vectors of the reaction mechanism include the minimizers of $H$ (i.e. one-step partial equilibria). Of course, this is a much simpler task than analysis of the signs of ${\mathrm{d}}G/{\mathrm{d}}t$ for all states and all values of parameters. Nevertheless, this simple check gives necessary and sufficient conditions for a function to be a conditionally universal Lyapunov function for the kinetic equations with a given reaction mechanism.

A cone of possible velocities is a convenient tool for analysis of conditionally universal Lyapunov functions for reaction networks. Consider all systems with detailed balance, a given reaction mechanism and a positive equilibrium $c^{\rm eq}$. According to the mass action law kinetic equations with detailed balance (15) the time derivatives ${\mathrm{d}}c/{\mathrm{d}}t$ at given point $c$ form a cone $\mathbb{Q}_{\rm DB}(c)$:

[TABLE]

where $cone$ stands for the conic hull and $sign(x)=1$ if $x>0$, $-1$ if $x<0$ and 0 if $x=0$.

Remark 1** (Differential criterion of conditionally universal Lyapunov functions).**

A smooth function $G(c,c^{\rm eq})$ is a conditionally universal Lyapunov functions for a given reaction mechanism and positive equilibrium $c^{\rm eq}$ if

[TABLE]

For continuous convex functions we have to use in (32) the subgradients of $G(c,c^{\rm eq})$ instead of gradients (see also the proof of Lemma 1).

The systems of complex balance are defined by linear relationships (26) between non-negative values $w_{r}^{\rm eq}$. Therefore, the time derivatives ${\mathrm{d}}c/{\mathrm{d}}t$ at given point $c$ for all systems with complex balanced equilibrium $c^{\rm eq}$ also form a cone, $\mathbb{Q}_{\rm CB}(c)$.

It is obvious that $\mathbb{Q}_{\rm DB}\subseteq\mathbb{Q}_{\rm CB}$. Surprisingly, these cones coincide [35, Theorem 2 (Local equivalence of detailed and complex balance)]. This means that for every mass action law system with the complex balanced positive equilibrium and any concentration vector $c$ there exists a mass action law system with detailed balance and the same positive equilibrium such that the velocity vectors ${\mathrm{d}}c/{\mathrm{d}}t$ at point $c$ for these systems coincide. If some reactions of the complex balance system are irreversible, a reverse reaction should be added. Its rate constant will be zero for the original complex balance system and non-zero for the detailed balance system. This theorem was proven in [35] even for more general kinetic law, the generalized mass action law. A bit earlier, such theorem was proven for Markov processes [36]. Continuous-time Markov kinetics with a finite number of states and a given positive equilibrium was studied. This class of systems is significantly wider than the systems with detailed balance. Nevertheless, for an arbitrary probability distribution $P$ and a general system there exists a system with detailed balance and the same equilibrium that has the same velocity ${\mathrm{d}}P/{\mathrm{d}}t$ at point $P$. The results are extended to nonlinear systems with the generalized mass action law.

Remark 2** (Coincidence of conditionally universal Lyapunov functions for detailed and complex balance systems).**

According to [35, Theorem 2], a function $G(c,c^{\rm eq}$ is a Lyapunov function for all systems (6) with given reaction mechanism and complex balanced equilibrium $c^{\rm eq}$ if and only if it is the universal Lapunov functions for all systems (15) with detailed balance, the same equilibrium and the same mechanism (supplemented by the reverse reactions, if necessary).

For some reaction mechanisms there exist Lyapunov functions without any relation to detailed balance or complex balance. For example, assume that all the elementary reactions have the form

[TABLE]

(only one $\alpha_{ri}$ can be non-zero; direct and reverse reactions are considered separately and some reactions can be irreversible). If there exists a positive balance $\sum_{i}m_{i}c_{i}=M=const$ then for any two solutions of the kinetic equations $c^{1}(t)$, $c^{2}(t)$ with the same value of $M$ the weighted $l_{1}$ distance between them $\sum_{i}m_{i}|c^{1}_{i}(t)-c^{2}_{i}(t)|$ monotonically decreases [37].

Convergent dynamics with quadratic Lyapunov norms $\|x\|^{2}=(x,Px)$, where $P$ is a symmetric positive definite matrix, was studied by Demidovich in 1960s and widely used [38]. Systems (33) give us example of a class of non-linear convergent systems in weighted $l_{1}$ norm. Some other examples of reaction mechanisms with such convergence property were produced in [39, 40] on the basis of monotonicity idea. All these selected mechanisms are rather simple. They have convergence property for any values of reaction rate constants. The monotonicity idea is useful for selection of reaction kinetic equations with stable dynamics without detailed or complex balance properties [39, 41].

On the contrary, in this paper we consider reaction networks with an arbitrary (presumably, nonlinear) reaction mechanism but with specific restrictions on the reaction rate constants. They should obey the principle of detailed balance. (Later on we explain why the same results are valid for systems with the so-called complex balance.) For such systems, there exist thermodynamic Lyapunov functions. For perfect systems under isothermal isochoric conditions the explicit form of this function is presented by (2). The situation with linear kinetics was similar when Rényi revealed $f$-divergences for Markov chains. The decrease of relative entropy (information) in time was well-known but there were no other Lyapunov functions until Rényi work [7]. Below we construct a wide family of additional Lyapunov functions for any nonlinear reaction network, obeying the mass action law and the principle of detailed balance.

4 Partial equilibria and new Lyapunov functions for mass action law

Let vector $\gamma$ have both positive and negative components. For every vector of concentrations $c$ we define the corresponding partial equilibrium in direction $\gamma$ as

[TABLE]

This partial equilibrium $c_{\gamma}^{*}(c)$ is the minimizer of $H$ on the interval

[TABLE]

This interval is bounded. For a positive point $c$ the minimizer $c_{\gamma}^{*}(c)$ is also positive. This is an elementary consequence of the logarithmic singularity of $(c\ln c)^{\prime}$ at zero. Here and below, argmin is the set of points where the function gets its minimum. The functions $H(c)$ is strongly convex on each bounded set because its Hessian has the form

[TABLE]

where $\delta_{ij}$ is the Kronecker delta. Therefore, each the argmin set in (34) consists of one point.

For monomolecular and for bimolecular reactions there are simple analytic expression for partial equilibria. Consider a monomolecuar reaction $A_{i}\rightleftharpoons A_{j}$. The non-zero components of the stoichiometric vector $\gamma$ are: $\gamma_{i}=-1$, $\gamma_{j}=1$. For a given vector $c$, the partial equilibrium is given by the equation $k^{+}c_{i}^{*}(c)=k^{-}c_{j}^{*}(c)$ under condition that $c^{*}(c)=c+x\gamma$. Simple algebra gives:

[TABLE]

Other components of $c^{*}(c)$ coincide with those of $c$. The sum $c_{i}+c_{j}=b$ does not change in the reaction $A_{i}\rightleftharpoons A_{j}$. Rewrite (35) using this ‘partial balance’ $b$:

[TABLE]

The ‘rate constants’ in (35) and (36) can be defined through a positive equilibrium point $c^{\rm eq}$: for the linear reaction, $A_{i}\rightleftharpoons A_{j}$, $k^{+}c_{i}^{\rm eq}=k^{-}c_{j}^{\rm eq}$ and we can take, for example,

[TABLE]

(we use the normalization condition $k^{+}+k^{-}=1$ to select one solution from the continuum of proportional sets of constants). The expression for the partial equilibrium for the linear reaction (36) is

[TABLE]

For a bimolecular reaction $A_{i}+A_{j}\rightleftharpoons A_{k}$ the non-zero components of the stoichiometric vector $\gamma$ are: $\gamma_{i}=\gamma_{j}=-1$, $\gamma_{k}=1$. Two independent ‘partial balances’ that do not change in the reaction are:

[TABLE]

The partial equilibrium $c^{*}(c)$ is the positive solution of the equation $k^{+}c_{i}^{*}(c)c_{j}^{*}(c)=k^{-}c_{k}^{*}(c)$ under condition that $c^{*}(c)=c+x\gamma$. After solving of quadratic equation for $x$ we get:

[TABLE]

The signs in front of square root are selected to provide positivity of $c^{*}(c)$.

We can rewrite (38) using any positive equilibrium $c^{\rm eq}$. Indeed, at this point, the detailed balance gives $k^{+}c^{\rm eq}_{i}c^{\rm eq}_{j}=k^{-}c^{\rm eq}_{k}$ and

[TABLE]

Therefore,

[TABLE]

Similar formulas can be easily obtained for the bimolecular reactions $A_{i}+A_{j}\rightleftharpoons A_{k}+A_{l}$ and $A_{i}+A_{j}\rightleftharpoons 2A_{k}$. They require nothing more than the quadratic formula and the detailed balance condition.

Partial equilibria appear in non-equilibrium thermodynamics from the very beginning. Already Jaynes considered conditional maximization of entropy as a basic method of equilibrium and non-equilibrium statistical physics [42]. According to Grmela, the equilibrium and nonequilibrium thermodynamics as well as the equilibrium and nonequilibrium statistical mechanics should be considered as particular representations of the Dynamical Maximum Entropy Principle [43]. Conditional maximization of entropy gives thermodynamic basis to thermodynamic of driven systems [44]. Partial equilibria (or quasiequilibria - conditional maximizers of entropy for more general conditions) are often considered as constrained equilibria as a result of introducing external or internal constraints. For example, the models of extended irreversible thermodynamics can be produced by the conditional maximization of entropy subject to various dynamical constrains and hypotheses about slow and fast variables [45]. Methods of invariant manifolds produces dynamical correction to these conditional maximum entropy models [46].

In this work, we use a very simple case of partial equilibria: these are equilibria with constrains under which only one reversible elementary reaction is on, while other reactions are frozen. To be more precise, we use the term ‘$H$-function’ instead of ‘entropy’.

The partial equilibrium $H$-function is associated with the partial equilibria:

[TABLE]

We do not assume here any time separation and constrained dynamics and use partial equilibria in a completely different way. Dynamics of $c(t)$ remains undeformed and foll;owes the original kinetic equations. For each stoichiometric vector $\gamma$ we consider the partial equilibrium $c^{*}_{\gamma}(c(t))$ as a projection of the genuine dynamics on the hypersurface of partial equilibria. This set of projections is a ‘shadow’ of kinetic curves. The new Lyapunov function for the genuine kinetics is combined from the entropies of these projections. For every compact set $\Gamma$ of vectors $\gamma$ with both positive and negative components, we define

[TABLE]

This definition has a simple explanation: for a given initial positive concentration vector $c$ and each stoichiometric vector $\gamma\in\Gamma$ we find the partial equilibrium of the one-step system ($m=1$) with this stoichiometric vector. This partial equilibrium, $c_{\gamma}^{*}(c)$, is the projection of the initial vector $c$ parallel to the vector $\gamma$ onto the hypersurface of partial equilibria defined by the equation

[TABLE]

From all these projections ($\gamma\in\Gamma$) we select the most non-equilibrium state, i.e., the state with the maximal value of $H=H(c_{\gamma}^{*}(c))$ (40). This maximal value is the new function $H_{\Gamma}^{*}(c)$ (41). It is defined by the set $\Gamma$ that should include all the stoichiometric vectors of the reaction mechanism.

It is necessary to stress that for every reaction mechanism there exists continuum of compact sets $\Gamma$, which include the stoichiometric vectors of this mechanism. The partial equilibrium is the same for the vector $\gamma$ and $x\gamma$ for any $x\neq 0$. Therefore, rigorously speaking, the functions $H_{\Gamma}^{*}(c)$ should be indexed by subsets of the projective space (the space of one-dimensional subspaces), and not by sets of vectors.

The hypersurface of partial equilibria for a given vector $\gamma$ is the hyperplane orthogonal to $\gamma$ in the entropic inner product. The level sets of the partially equilibrium function $H_{\gamma}^{*}(c)$ are cylinders with the spheric base and axis parallel to $\gamma$.

Fig. 1 represents a very simplified example of the partial equilibrium $c_{\gamma}^{*}(c)$ and the Lyapunov function $H_{\Gamma}^{*}(c)$ construction. The system with three components of the same molecular weight is presented in the triangle $c_{1}+c_{2}+c_{3}=const$ drawn in barycentric coordinates. The reaction mechanism consists of three reactions $A_{1}\rightleftharpoons A_{2}$,$A_{2}\rightleftharpoons A_{3}$, and $2A_{1}\rightleftharpoons A_{2}+A_{3}$. The equilibrium $c^{\rm eq}$ is assumed in the center of the triangle ($c^{\rm eq}_{1}=c^{\rm eq}_{2}=c^{\rm eq}_{3}$). The partial equilibria of the first two reactions form the straight lines in the triangle, the medians, while the partial equilibria of the non-linear reaction form a parable $c_{2}c_{3}/c_{1}^{2}=1$. These three lines intersect in the equilibrium due to detailed balance.

Explicit expressions for the partial equilibria $c^{*}_{\gamma}(c)$ are:

For the reaction $A_{1}\rightleftharpoons A_{2}$

[TABLE]

• 2.

For the reaction $A_{2}\rightleftharpoons A_{3}$

[TABLE]

• 3.

For the reaction $2A_{1}\rightleftharpoons A_{2}+A_{3}$

[TABLE]

where $b_{1}=c_{1}+c_{2}+c_{3}$ and $b_{2}=c_{2}-c_{3}$ are the independent ‘partial balances’ for this reaction.

Fig. 1a shows the partial equilibria for an arbitrarily selected point $c$. In Fig. 1b, one level set of $H_{\Gamma}^{*}$ is presented, where $\Gamma=\{\gamma_{1},\gamma_{2},\gamma_{3}\}$:

[TABLE]

$c^{*}_{\gamma_{i}}$ $i=1,2,3$ are given by (42), (43), (44) and $H$-function is the classical one (2).

The sublevel set of $H_{\Gamma}^{*}(c)$ is the intersection of strips with sides parallel to the stoichiometric vectors $\gamma_{i}$. These strips are sublevel sets for the partial equilibrium entropies $H_{\gamma_{i}}^{*}(c)$. In higher dimensions, the level sets of partial equilibrium $H$-function $H_{\gamma}^{*}(c)$ are cylindrical hypersurfaces with the generatrix parallel to $\gamma$. The base (or the directrix) of this cylindric surface is the level set of $H$ on the surface of partial equilibrium. In Fig. 1b, the ‘surfaces’ of partial equilibria are lines, the level sets of $H$ on these lines are couples of points. These points are highlighted. Note that the sublevel areas for the new function $H^{*}_{\Gamma}$ in Fig. 1 are convex polygons, whereas for the classical $H$-function they have smooth border.

Assume that the principle of detailed balance holds and the vector of equilibrium concentrations $c_{i}^{\rm eq}$ in the definition of $H$-function (2) is the point of detailed balance.

Theorem 2**.**

Let all the stoichiometric vectors of the reaction mechanism (4) belong to $\Gamma$ then $H_{\Gamma}^{*}(c)$ (40) is monotonically non-increasing function along the solutions of the kinetic equations (15) for all positive $c$ and all non-negative values of equilibrium fluxes $w^{\rm eq}$.

Proof.

To use Theorem 1 we have to prove two statements:

The function $H_{\Gamma}^{*}(c)$ is convex in $\mathbb{R}_{>0}^{n}$; 2. 2.

For each $\gamma\in\Gamma$ and a positive vector $c$ the minimizer of $H$ on the interval $(c+\gamma\mathbb{R})\cap\mathbb{R}_{>0}^{n}$ is, at the same time, a minimizer of $H_{\Gamma}^{*}$ on this interval.

We prove convexity of the function $H_{\Gamma}^{*}(c)$ in $\mathbb{R}_{>0}^{n}$ in two steps.

Convexity of $H_{\gamma}^{*}(c)$ for one-element sets $\Gamma=\{\gamma\}$.

• 2.

Convexity of $H_{\Gamma}^{*}(c)=\max_{\gamma\in\Gamma}H_{\gamma}^{*}(c)$.

Let us prove convexity of $H_{\Gamma}^{*}(c)$ for one-element sets $\Gamma$, that is, we will prove convexity of the partial equilibrium $H$-function $H_{\gamma}^{*}(c)$ (40). By definition, convexity of $H_{\gamma}^{*}(c)$ in $\mathbb{R}_{>0}^{n}$ means that for each two positive concentration vectors $c^{1}$ and $c^{2}$ and a number $\lambda\in[0,1]$ the inequality holds:

[TABLE]

First, notice that due to the convexity of $H$

[TABLE]

Secondly, $H(c_{\gamma}^{*}(c^{1,2}))=H_{\gamma}^{*}(c^{1,2})$ by the definition of the partial equilibrium $H$-function (40). Therefore, the previous inequality can be rewritten as

[TABLE]

Finally,

[TABLE]

hence,

[TABLE]

and

[TABLE]

because the last value is the minimum of $H$ on the interval

[TABLE]

Convexity of the partial equilibrium $H$-function $H_{\gamma}^{*}(c)$ is proven.

Convexity of $H_{\Gamma}^{*}(c)$ follows from the convexity of the partial equilibrium $H$-function $H_{\gamma}^{*}(c)$, from the definition of $H_{\Gamma}^{*}(c)$ as the maximum of $H_{\gamma}^{*}(c)$ ($\gamma\in\Gamma$), and from the following fact from convex analysis: Maximum of a set of convex functions is again convex. The shortest proof is based on the definition of a convex function as a function with convex epigraph [47] and follows from the observation that the epigraph of the maximum of a family of functions is the intersection of their epigraphs.

Let us analyse the minimizers of $H_{\Gamma}^{*}$ on the interval $I=(c+\gamma\mathbb{R})\cap\mathbb{R}_{>0}^{n}$ for a positive concentration vector $c$. Select $\gamma\in\Gamma$. The minimizer of $H$ on the interval $I$ is the partial equilibrium $c_{\gamma}^{*}(c)$, by the definition (34). Function $H_{\gamma}^{*}(c)$ is constant on the interval $I=(c+\gamma\mathbb{R})\cap\mathbb{R}_{>0}^{n}$ and $H_{\gamma}^{*}(c)=H(c_{\gamma}^{*}(c))$ on $I$. Therefore, $c_{\gamma}^{*}(c)$ is a minimizer of $H_{\gamma}^{*}(c)$ on $I$ (trivially, as all other points of $I$ do).

Notice, that for all $\gamma^{\prime}\in\Gamma$

[TABLE]

and the equality is strong if $c_{\gamma}^{*}(c)\neq c_{\gamma^{\prime}}^{*}(c_{\gamma}^{*}(c))$. Indeed, $c_{\gamma^{\prime}}^{*}(c_{\gamma}^{*}(c))$ is the unique minimizer of $H$ on the interval $(c_{\gamma}^{*}(c)+\gamma^{\prime}\mathbb{R})\cap\mathbb{R}_{>0}^{n}$. If this minimizer does not coincide with $c_{\gamma}^{*}(c)$ then

[TABLE]

According to inequality (45),

[TABLE]

For any family of convex functions $\Phi$ on the interval $I$ the following statement holds. Let $f\in\Phi$. If $y$ is a minimizer of $f(x)$ on $I$, $F(x)=\max\{\phi(x)|\phi\in\Phi\}$, and $F(y)=f(y)$ then $y$ is a minimizer of $F(x)$ on $I$. Indeed, for any $z\in I$ $f(z)\geq f(y)$ because $y$ is a minimizer of $f$. At the same time, $F(z)\geq f(z)$ by the definition of $F$. Hence, for any $z\in I$ $F(z)\geq F(y)$ and $y$ is a minimizer of $F(x)$ on $I$.

Let us take $f=H_{\gamma}^{*}(c)$ and $\Phi=\{H_{\gamma^{\prime}}^{*}(c)|\gamma^{\prime}\in\Gamma\}$. Then $F=H_{\Gamma}^{*}(c)$. Select $y=c_{\gamma}^{*}(c)$. Notice that $y$ is a minimizer of $f$ and $F(y)=f(y)$ according to (46). Therefore, $c_{\gamma}^{*}(c)$ is a minimizer of $H_{\Gamma}^{*}(c)$ on $I$. By combining this result with the proven convexity of $H_{\Gamma}^{*}(c)$ and applying Theorem 1, we prove that $H_{\Gamma}^{*}(c)$ is a monotonically non-increasing function on solutions of kinetic equations. ∎

Corollary 1**.**

Assume that the complex balance condition holds and the vector of equilibrium concentrations $c_{i}^{\rm eq}$ in the definition of $H$-function (2) is the point of complex balance.

Let all the stoichiometric vectors of the reaction mechanism (4) belong to $\Gamma$ then the time derivative $H_{\Gamma}^{*}(c)$ (41) by virtue (15) is non-positive for all positive $c$ and all non-negative values of equilibrium fluxes $w^{\rm eq}$: ${\mathrm{d}}H_{\Gamma}^{*}(c)/{\mathrm{d}}t\leq 0$.

Proof.

This corollary follows immediately from Theorem 2 and from [35, Theorem 2 (Local equivalence of detailed and complex balance)] because $\mathbb{Q}_{\rm DB}(c)=\mathbb{Q}_{\rm CB}(c)$ for all positive $c$. ∎

When $c(t)$ approaches an equilibrium, $c^{\rm eq}$, all the partial equilibria also converge to $c^{\rm eq}$. Locally, near the equilibrium $c^{\rm eq}$, this picture can be simplified. Let us use quadratic approximation to $H(c)$. The entropic inner product is defined as $\langle x|y\rangle=(x,(D^{2}H)_{c^{\rm eq}}y)$, where $(D^{2}H)_{c^{\rm eq}}$ is the second differential of $H$ at the equilibrium (Hessian matrix) and $(\ ,\ )$ is the standard inner product.

Remark 3** (Universal Lyapunov functions for linearized kinetic equations).**

In the entropic inner product, the quadratic approximation to the classical $H$-function near equilibrium $c^{\rm eq}$ is

[TABLE]

For this approximation, elementary linear algebra gives

[TABLE]

Formula (47) gives universal Lyapunov functions for linearized kinetic equations with the given reaction mechanism and detailed balance conditions.

Thus, the new family of universal Lyapunov functions for chemical kinetic equations with detailed or complex balance is constructed.

5 Outline of possible generalizations and applications

In this section, we briefly outline the possible generalizations and applications of the main results. This is basically a ‘to do’ section and detailed analysis of all generalizations and applications is beyond the scope of the paper.

Generalized mass action law

Analysis of the proven statements demonstrate that two properties of $H$ are used: strict convexity in directions $\gamma$ in positive orthant $\mathbb{R}_{>0}$ and logarithmic singularity of the $\partial H/\partial c_{i}$ when $c_{i}\to 0$. One more requirement is the generalized mass action laws for the reaction rates:

[TABLE]

where $\exp(\nabla H,\alpha_{r})$ is the ‘Boltzmann factor’, and $\phi_{r}\geq 0$ is a non-negative quantity, ‘kinetic factor’ [31]. Detailed balance and complex balance conditions for the generalized mass action law have the same form as for perfect systems. For detailed balance such kinetics was introduced by Feinberg [48] (Marselin–De Donder kinetics, see also [3]). Grmela studied properties of geometry of nonlinear non-equilibrium thermodynamics for the generalized mass action law [49]. Detailed balance in the limit of non-equilibrium reactions was analyzed in [50]. General analysis of non-classical entropies and their relations to the second law of thermodynamics was presented in work [35]. Complex balance conditions for macroscopic kinetics was proven for systems with Markov microscopic description in [31]. The basic constructions of partial equilibria and functions $H^{*}_{\Gamma}$ can be extended to the generalized mass action law.

Non-isochoric and non-isothermal conditions

The Lyapunov function (2) and kinetic equations (6) with mass action low (5) are valid for isochoric isothermal conditions (constant volume and temperature). For other conditions they should be modified. Such a modification is rather simple and does not add any substantial change but the equation and Lyapunov functions have different analytic form (see, for example, [14, 51]). The analysis of Lyapunov functions for these conditions can be provided in the formalism of the generalized mass action law. A general form for the description of non-isothermal reactions in closed chemical systems in terms of the Marcelin–De Donder kinetics and explicit form of the $H$-function for the systems with detailed balance under various conditions is presented in [52].

Spatially distributed systems

Transport in Boltzmann’s equation is conservative (the free flight) and dissipative terms are local (the collision integral). Generalization of our approach to the collision integral has to be done. In general, models of complex transport processes can include both dissipative and conservative terms. There are many attempts to create thermodynamic theory of such processes. In the GENERIC approach, conservative and dissipative components are explicitly separated with some commutativity conditions between them [53]. Generalization of proposed construction onto dissipative components of transport processes seems to be a challenging task. If the dissipative part is described in the language of quasi-chemical formalism with a finite number of steps [54] then this generalization is more straightforward.

Possible applications

What can add the new Lyapunov functions to research tools? Of course, more Lyapunov functions are better than less. The non-classical Lyapunov functions for linear systems are widely used for various estimates and information analysis in the situations, where linear Markov chains could serve as adequate models of information transformation. Just for example of various applications, we can refer to works [55, 56, 57]. Universal Lyapunov functions are instruments for evaluation of possible dynamics when the reaction rate constants are unknown or highly uncertain. Without any knowledge of the reaction mechanism we use the thermodynamic potentials for evaluation of the attainable sets of chemical reactions (the theory and algorithms are presented in [15, 19, 58], some industrial applications are discussed in [18, 59, 60]). Convexity allows us to transform the $n$-dimensional problems about attainability and attainable sets into an analysis of one-dimensional continua and discrete objects, thermodynamic trees [19]. When the reaction mechanism is known, we can use this information about mechanism for sharper estimations. A new class of estimates is needed and the new Lyapunov functions give a collection of instruments for such estimates. The use of many Lyapunov functions in the analysis of attainable sets makes these estimates more narrow and close to reality.

6 Conclusion

6.1 New Lyapunov functions: the main result

A rich family of universal Lyapunov functions for any linear or non-linear reaction network with detailed or complex balance is presented. Consider a mass action law system with a given reaction mechanism, the set of stoichiometric vectors $\Upsilon$ and detailed or complex balance. For construction of the new Lyapunov functions for this system the following operations are needed.

Normalize vectors from $\Upsilon$. Let $\Upsilon_{0}$ be the set of normalized stoichiometric vectors. 2. 2.

Select a finite set of vectors $\Gamma$ such that $\Upsilon_{0}\subset\Gamma$ and each vector from $\Gamma$ has both positive and negative components. 3. 3.

For each $\gamma\in\Gamma$ calculate

[TABLE] 4. 4.

Find

[TABLE]

Thus, we have to perform several smooth one-dimensional convex minimizations and then select the maximum of these minima.

Theoretically, $H^{*}_{\Gamma}$ is a Lyapunov function for any compact $\Gamma$. Such a function can be approximated (from below) by the functions calculated for finite $\varepsilon$-networks on $\Gamma$.

Functions $H^{*}_{\Gamma}$ are not presented as elementary functions even if $H$ has explicit analytical form. If the reaction mechanism consists of mono- and bimolecular reactions then the partial equilibria can be represented by radicals and $H^{*}_{\gamma}(c)$ can be written as an elementary function but in general case the less implicit form with min and max operations is necessary. If we substitute $H$ by its quadratic approximation near equilibrium then functions $H^{*}_{\Gamma}$ can be represented more explicitly (47) and the difference between $H$ and $H^{*}_{\Gamma}$ becomes obvious.

The mystery about the fundamental difference between the rich family of Lyapunov functions for linear networks and a very limited collection of Lyapunov functions for non-linear networks in thermodynamic conditions is resolved: there is no such crucial difference anymore.

6.2 An open problem

There remains an obvious difference between explicit analytic expression of $f$-divergences (3) and not so obvious construction of Lyapunov functions for general networks using partial equilibria of non-linear reactions. For linear reactions, the partial equilibria have very simple analytic expression, for bimolecular reactions they are given explicitly using quadratic formula, but for trimolecular reactions the analytic formulas become too bulky.

Alt least one important question is still open. The new Lyapunov functions $H_{\Gamma}^{*}$ are, at the same time, universal Lyapunov functions for linear kinetics, if the stoichiometric vectors of the linear reaction mechanism $A_{i}\rightleftharpoons A_{j}$ are included in $\Gamma$. Due to the results of [10, 11], such a function should be, essentially, a $f$-divergence (3) $H_{f}(c)$, or, more precisely, it should be a monotonic function of $H_{f}(c)+\lambda_{i}c_{i}$ for some constant $\lambda$. Nevertheless, now we know nothing about these $f$-divergences except their existence. Constructive transformation of $H_{\Gamma}^{*}$ into $f$-divergence is desirable because an explicit form (3) brings some benefits for analysis.

Acknowledgments

The work was supported by the University of Leicester and the Ministry of Science and Higher Education of the Russian Federation (Project No. 14.Y26.31.0022).