Universal Gorban’s Entropies: Geometric Case Study

Recently, A.N. Gorban presented a rich family of universal Lyapunov functions for any linear or non-linear reaction network with detailed or complex balance. Two main elements of the construction algorithm are partial equilibria of reactions and convex envelopes of families of functions. These new functions aimed to resolve “the mystery” about the difference between the rich family of Lyapunov functions (f-divergences) for linear kinetics and a limited collection of Lyapunov functions for non-linear networks in thermodynamic conditions. The lack of examples did not allow to evaluate the difference between Gorban’s entropies and the classical Boltzmann–Gibbs–Shannon entropy despite obvious difference in their construction. In this paper, Gorban’s results are briefly reviewed, and these functions are analysed and compared for several mechanisms of chemical reactions. The level sets and dynamics along the kinetic trajectories are analysed. The most pronounced difference between the new and classical thermodynamic Lyapunov functions was found far from the partial equilibria, whereas when some fast elementary reactions became close to equilibrium then this difference decreased and vanished in partial equilibria.


Classical Entropic Lyapunov Functions for General Kinetics
The classical example of the Lyapunov functional in kinetics was provided by Boltzmann in 1872 [1] (twenty years before the famous Lyapunov thesis): where f (x, v) is an one-particle distribution function in space (x) and velocity (v). The analogue of this functional for chemical reaction was known already for Gibbs [2]: where c i ≥ 0 is the concentration of the ith component A i and c eq i > 0 is an equilibrium concentration of A i (under the standard convention that x ln x = 0 for x = 0). This is the thermodynamic potential for systems under constant temperature and volume (up to a constant factor).
In 1938, Zeldovich [3] used convexity of function (2) and logarithmic singularity of its derivatives at zeros for his proof of uniqueness of positive chemical equilibrium for given values of linear balances. In the 1960s, this approach was applied for many systems under different conditions and became standard [4].
For systems with detailed balance, the time derivative of H is the sum (or integral, for continua of elementary processes) of the terms: where w + and w − are the rates of the direct and reverse elementary process, respectively, and the term (w + − w − ) ln(w + /w − ) ≥ 0 is the entropy production in an elementary process. Boltzmann used principle of detailed balance in the proof of his H-theorem in 1872, but in 1887 he invented a remarkable generalization of his theorem (after criticisms by Lorentz) [5]. His new sufficient condition for H-theorem, the cyclic balance or the semidetailed balance, was several times rediscovered later on. In chemical kinetics, it is called 'the complex balance' [6]. For linear kinetics, the generalisation from detailed balance to complex balance is equivalent to the generalisation from the reversible Markov chains to general Markov chains (with positive equilibrium). For non-linear kinetics this condition seems to be more restrictive (it will be discussed below in more detail).
Shannon proved an analogue of the H-theorem for general random manipulation with information (for Markov chains, essentially). This is the information processing lemma [7].
The classical Lyapunov functions (1) and (2) have an important property, universality: they do not depend directly on the collision and reaction mechanisms and kinetic constants but on the equilibrium distributions (concentration and the detailed or complex balance condition in general form [8]). This universality can be considered as a manifestation of the universality of thermodynamics that does not depend on the microscopic details directly.

General Lyapunov Functions for Linear Kinetics
In 1960, an extremely rich family of Lyapunov functions was discovered for general Markov chains. Rényi [9] proved that the following functions ( f -divergences) are the Lyapunov functions for general linear kinetics (Markov chains) with positive equilibrium c where f is an arbitrary convex function on the positive semi-axis.
Moreover, H f are not just Lyapunov functions but divergences: is monotonically non-increasing function of time t for any two kinetic curves c 1 (t) and c 2 (t) with the same value of ∑ i c i . This discovery attracted less immediate attention than the Rényi entropy proposed in the same paper. (Here, P here is a vector of probability distribution with coordinates p i .) Nevertheless, a bit later the theory of f -divergences was recognised as an important instrument of information theory and kinetics [10,11]. In 2003, P.A. Gorban proved in that all universal Lyapunov functions for Markov kinetics can be produced by monotonic transformations of f -divergences [12]. In 2009, Amari [13] got a similar result.

Conditionally Universal Lyapunov Functions for General Kinetics
The f -divergences (4) are universal Lyapunov functions as they do not depend on kinetic constants directly but on the equilibrium only. Nevertheless, their universality is weaker than the universality of the classical thermodynamic Lyapunov functions like (2) because the classical thermodynamic potentials change monotonically in time for any reaction mechanism, linear or non-linear, under conditions of detailed or complex balance, whereas f -divergences are defined for linear kinetics only: for the sets of elementary processes like A i A j , where A i are the components (or states). If a function changes monotonically in time for a given reaction mechanism under conditions of detailed or complex balance, then we call it a conditionally universal Lyapunov function for this reaction mechanism [8].
For linear reaction mechanisms a rich family of conditional Lyapunov functions (4) is known since 1960 [9]. Nevertheless, there were no general constructions of conditionally universal Lyapunov functions for non-linear reaction mechanisms till the series of works [8,14,15], where new conditionally universal Lyapunov functions were constructed for an arbitrary reaction mechanism under detailed or complex balance condition. These functions differ from the classical thermodynamic potentials, by construction. Nevertheless, it could be important to analyse how different they are. For this purpose, in this paper we compare the level sets of these functions and their changes over time for several typical chemical reaction examples.

Structure of the Paper
The basic notions and generalised mass action law equations are systematically introduced in Section 2. The time derivative of the thermodynamic Lyapunov functions is calculated explicitly for systems with detailed balance. The conditionally universal Lyapunov functions are characterised in this section implicitly, through their geometric properties. An extension of the general results to systems with complex balance is given in Section 2.5. The explicit construction and algorithm for calculation of Gorban's Lyapunov functions are described in Section 3. Section 4 is devoted to the case studies and comparative analysis of the level sets and dynamics of classical and Gorban's Lyapunov functions for several reaction kinetic systems. The results and outlooks are summarized in Conclusion.

Generalised Mass Action Law
The construction of the Generalised Mass Action Law (GMAL) kinetic equations uses several basic elements:

•
The list of components that is a finite set of symbols A 1 , . . . , A m ; The list of elementary reactions (the reaction mechanism) that is a finite set of the stoichiometric equations where ρ = 1, . . . , m is the reaction number and the stoichiometric coefficients α ρi , β ρi are nonnegative real numbers; • A dimensionless free entropy S(N) that is a concave function in R n ≥0 .
We use the following notations: α ρ , β ρ are the vectors with coordinates α ρi , β ρi , respectively, γ ρ = β ρ − α ρ is the stoichiometric vector of the reaction (5) (the 'gain minus loss' vector). The definition of the dimensionless free entropy function for a physico-chemical system depends on the conditions. For isolated systems it is just the thermodynamic entropy divided by the gas constant R. For isothermal isochoric conditions S = −F/(RT), where F is the Helmholtz free energy, T is the temperature, for isothermal isobaric conditions S = −G/(RT), where G is the Gibbs energy (free enthalpy), etc. [16][17][18]. Introduction in the theory of thermodynamic potentials including free entropies (Massieu-Plank functions) is given by Callen [19].
For the general GMAL construction, S is just a concave function. For the sake of generality, the value S = −∞ is also allowed. The function H = −S is assumed to be a closed convex function, this means that the its sublevel set {N ∈ R n ≥0 |H(N) ≤ a} is a closed set for any real a. It is also assumed that H takes finite values on a convex domain U ⊂ R n ≥0 with non-empty interior. H is twice differentiable almost everywhere in U (A.D. Alexandrov theorem [20,21]). Following Boltzmann's tradition, we will use further the H-function H = −S.
A non-negative quantity, the reaction rate r ρ is defined by GMAL almost everywhere in U for every elementary reaction (5) [15,22] (compare to the thermodynamic GMAL presentations of reaction rates in earlier works [16,[23][24][25]): where the kinetic factor ϕ ρ ≥ 0 is an intensive quantity.
Here and below, all the equalities and inequalities with gradients of H are considered 'almost everywhere' in U if the convex function H is not everywhere continuously differentiable.
For the perfect isothermal isochoric mixtures H-function has the form where c i = N i /V and c * i = const. For such systems, GMAL (6) becomes the standard mass action law: The corresponding GMAL kinetic equation is where V > 0 is a positive extensive variable (volume). It can also change with time and its dynamic is defined by the equation of state and by the conditions of the process. The structure of kinetic Equation (9) and GMAL formula for reaction rates (6) allow the elegant expression for dH/dt. Let an auxiliary function of real variable θ(λ) be given by the following expression for a given composition vector N [16,22,26]: Function θ(λ) is convex. With this function, dH/dt has a very simple form: Convexity of θ(λ) implies the following sufficient condition of non-positivity dH/dt.
General kinetic Equation (9) with GMAL reaction rate (6) can describe arbitrarily complex dynamics and approximate any dynamical system in U even for perfect mixtures and constant kinetic factors [27]. The specific thermodynamic properties of kinetic equations are based on special relations between kinetic factors ϕ ρ that are detailed balance and complex balance.

Detailed Balance
The principle of detailed balance is a special symmetry between direct and reverse elementary reactions caused by the so-called microreversibilty (invariance of the equations of microscopic dynamics with respect to time reversal). In the GMAL formalism, the principle of detailed balance has a simple form: kinetic factors of direct and reverse elementary reaction coincide. In such situations, it is convenient to rearrange the list of elementary reactions (5), join the reactions with their reverse reactions in a shorter list of pairs of reactions: If the reverse reaction does not exist in the original reaction mechanism (5) then we can, nevertheless, add the reverse reaction formally, with zero kinetic factor. For the reaction mechanism in the reversible form (12), we use the superscripts + and − for the reaction rates and kinetic factors of the direct and reverse reactions, respectively: The rate r ρ is defined as the difference r ρ = r + ρ − r − ρ and the kinetic equations have the same form (9). The detailed balance condition is: Under this condition, a symmetry relation holds: θ(λ) = θ(1 − λ). Therefore, θ(1) = θ(0) for every composition vector N and according to Proposition 1, dH/dt ≤ 0. Direct calculation of dH/dt by virtue of the system of kinetic equations under the detailed balance condition gives the classical result (compare to (3)): dH dt Because of this property, H(N) is called the thermodynamic Lyapunov function. The detailed analysis of entropy production in nonequilibrium systems was provided recently in [28]. Grmela considered the equilibrium and nonequilibrium thermodynamics as representations of the Dynamical Maximum Entropy Principle [29].

Conditionally Universal Lyapunov Functions and Their Geometric Characterisation
In this Subsection, we consider systems of kinetic Equation (9) with the given thermodynamic Lyapunov function H, reaction rates presented by GMAL (6), and detailed balance (14) for a given reaction mechanism (12). According to inequality (15), H is a Lyapunov function for such a system for any reaction mechanism. This means that H is a universal Lyapunov function for chemical kinetics. If the reaction mechanism is fixed then the conditionally universal Lyapunov functions are introduced.

Definition 1 ([15]).
A convex function F(N) in U is a conditionally universal Lyapunov function for kinetic Equation (9), given H and reaction mechanism (12) if for any values of kinetic factors, which satisfy the detailed balance conditions (14).
For each elementary reaction ∑ i α ρi A i ∑ i β ρi A i from the reaction mechanism given by the stoichiometric Equation (12) and any X ∈ U we define an interval of a straight line Definition 2 (Partial equilibria criterion for GMAL). A convex function F(N) on U satisfies the partial equilibria criterion with a given thermodynamic Lyapunov function H and reversible reaction mechanism given by stoichiometric Equation (12) if argmin for all X ∈ U, ρ = 1, . . . , m.

Theorem 1.
[General H-theorem]A convex function F(N) on U is a conditionally universal Lyapunov function for kinetic Equation (9), given H and reaction mechanism (12) if it satisfies the partial equilibria criterion (Definition 2).

Complex Balance
Let us return to the general form of the reaction mechanism without coupling direct and reverse reactions (5). The complex balance condition means that θ(1) ≡ θ(0) for all values of the gradient vectors from R n . More formally, it means that for all vectors µ ∈ R n with coordinates µ i . Functions exp(y, µ) of vector µ ∈ R n are linear independent for any finite set of y ∈ R n . Therefore, the identity (18) can be split in the several linear conditions on the coefficients ϕ ρ . Assume that there are q different vectors y 1 , . . . , y q among {α ρ , β ρ } (ρ = 1, . . . , m). The identity (18) is equivalent to q conditions: Formal sums ∑ y i A i from stoichiometric equations are called complexes, so conditions (19) are called the complex balance conditions [6]. In physics, the terms cyclic balance conditions or semidetailed balance conditions are also used. These conditions were derived from the Markov processes of microkinetics under two asymptotic assumptions: (i) the asymptotic intermediates are in fast equilibrium with the main components and (ii) the concentration of asymptotic intermediates is small (the Michaelis-Menten-Stueckelberg theorem [22]). If each complex ∑ y i A i is once and only once the left hand part of the stoichiometric equation from the reaction mechanism (5) and once the right hand part, for the reverse reaction equation, then the complex balance conditions literally coincide with the detailed balance conditions.

Cone Theorem and H-Theorems for Complex Balancing Systems
For analysis of conditionally universal Lyapunov functions, a notion of cone of possible velocities is useful [8,15,16]. This cone is defined for a cone of kinetic equations and a given composition vector N. It consists of all possible values of the velocity vector dN/dt at this point for equations from selected cone. For example, the systems with detailed balance for a given reaction mechanism and function H form the convex cone in the space of vector fields on the composition space. The corresponding cone of possible velocities is where cone stands for the conical hull and the piecewise-constant functions sgn(r ρ (N)) do not depend on (positive) values of kinetic factors ϕ ρ under assumption of detailed balance. Indeed, Consider the complex balance systems with a given reaction mechanism and function H. They are given by linear conditions (19) and form a convex cone of vector fields in the composition space. For a given composition vector N, the cone of all values of dN/dt is a convex cone in R n . We denote this cone by Q CB (N).
For a given function H, consider a reversible reaction mechanism (12) and kinetics with detailed balance. Calculate Q DB (N). Decouple the direct and reverse reactions, consider kinetics with complex balance (for the same reaction mechanism). Calculate Q CB (N). These cones coincide: [8,14,15]). For the same set of elementary reactions, This means that the possible directions of motion for the kinetic systems with detailed and for systems with complex balance at one point coincide. The difference between these two classes of systems appears if we consider several points or kinetic curves, not pointwise.
Time derivative of a function F(N) by virtue of kinetic equations is computed pointwise, therefore, an obvious consequence of Theorem 2 is: Corollary 1. If a function F(N) is a conditionally universal Lyapunov function for systems with detailed balance, thermodynamic Lyapunov function H and reaction mechanism (12) then it is a conditionally universal Lyapunov function for the systems with complex balance, the same H and the list of elementary reactions.
So, any construction of conditionally universal Lyapunov functions for systems with detailed balance can be easily generalised for systems with complex balance.

Gorban's Lyapunov Functions H Γ
Direct application of the general H-theorem (Theorem 1) gives the following construction of conditionally universal Lyapunov functions for GMAL kinetic equations with detailed and complex balance [8,14,15]. Consider a GMAL system with the reaction mechanism (12), the convex thermodynamic Lyapunov function H and the detailed or complex balance. Let Γ ⊂ R n be a finite set of non-zero vectors, which includes all the stoichiometric vectors γ ρ . Assume, additionally, that the function H is strictly convex on each non-empty interval U ∩ (N + Rγ) (γ ∈ Γ) and achieves its minimum on this interval in an internal point (this point of minimum is unique due to strict convexity of H in direction γ). This property trivially holds for the H function for perfect systems under isothermal isochoric conditions (7) as well as for perfect systems under all other classical conditions (for example, for isothermal isobaric systems or for isolated isochoric systems if γ has both positive and negative coordinates [17]).
Two main operations in the construction of the conditionally universal Lapunov function H Γ (N) are [8]: Thus, for calculation of H Γ (N) we have to solve several 1D convex minimization problems and select the maximum of these minima. These functions are indexed by finite set Γ. The construction of H Γ (N) does not depend on the length of the vectors γ ∈ Γ. Therefore, for theoretical purposes it makes sense to consider normalised vectors or, even better, the elements of the projective space (i.e., one-dimensional subspaces of R n ). For calculations, such a normalisation is not necessary.
The quasi-equilibrium entropies (21) and partial equilibria in direction γ are standard and very old tools for description of fast equilibria and partial equilibrium approximations. For example, the classical work of Michaelis and Menten used assumption of fast equilibration of 'compounds' with stable reagents [30]. For detailed discussion of this approximation we refer to [22], for application to thermodynamics of driven systems see [24], more physical and chemical applications, from Boltzmann's equation to chemical kinetics, and general theory are presented in the book [31]. Partial equilibria are used in the construction of Gorban's universal Lyapunov functions (22) in a completely different way. They do not substitute the genuine kinetic trajectory as the partial equilibrium approximations, but rather follow the non-perturbed motion as the ensemble of its projections on the surfaces of partial equilibria ('partial equilibrium shadows'). For the calculation of Gorban's function, the closest shadow is selected. (It is the closest shadow in the entropic divergence (22)). Which shadow is closest can be changed in the course of motion.
Such ensembles of quasi-equilibrium projections were used in 1979 [32] for construction of attainability regions for chemical kinetic equations with a given reaction mechanism (this problem is close to the problem of conditionally universal Lyapunov functions). Later on, this geometric approach was used in various applications [16,33] and reappeared recently in the theory of toric differential inclusions of chemical kinetics [34].
The 'ensemble of equilibrium subsystems' has been intensively used during almost 40 years as an effective tool for mathematical analysis of complex catalytic reactions and has given rise to many useful methods reviewed in a recent book [35].

Case Studies
"A picture is worth a thousand words." Examples are needed to evaluate the difference between Gorban's entropies and the classical Boltzmann-Gibbs-Shannon entropy. The difference in their construction is obvious but we need to evaluate the difference between these functions values and between their changes in dynamics. In this section, we analyse the level sets of these functions and their dynamic changes along kinetic trajectories. Several reaction mechanisms have been selected for benchmarking:

•
Linear isomerisation of three components (Section 4.2) • Nonlinear isomerisation reaction (Section 4.3) • Water Gas Shift (WGS) reaction (Section 4.4) or in abstract notations • Hydrogen Chloride (HCl) reaction (Section 4.5) For these reaction mechanisms, we selected various cortéges of reaction rate constants: with detailed balance, with complex balance, more or less stiff, etc. The goal was to demonstrate various aspects of similarity and difference between the classical thermodynamic Lyapunov functions and Gorban's functions.
All the systems below were considered in perfect gases and under isothermal isochorich conditions. Therefore, the volume V was constant and there was no need to use two sets of variables, amounts N i and concentrations c i . We used the concentrations c i with the vector of concentrations c and the classical thermodynamic Lyapunov function for these conditions H(c) (2).
The first subsection below contains explicit formulae for points of partial equilibrium for five types of elementary reactions which are used in case studies. The following four subsections present four case studies for four different reaction systems.

Partial Equilibria for Several Typical Reactions
The calculation of Gorban's Lyapunov function H Γ requires finding the points where γ is a stoichiometric vector or any other vector with at least one positive and at least one negative element. There is no general formula for the explicit search for such points, but for some typical cases an explicit solution can be found analytically.
Since Boltzmann's H is a strictly convex function in R n >0 with c i log c i singularities at the borders, the minimizer in the direction γ is a positive vector c + γχ, where dH(c + γχ)/dχ = 0: A partial equilibrium in direction γ satisfies the following equation: Equation (23) is very similar to the usual condition of detailed balance but we have to emphasise that it does not include any reaction rate constant, does not assume the reversibility of any reaction or microreversibility and just describes the minimisers of H in the given direction. It can be considered as the thermodynamic equilibrium condition for the elementary reaction with the stoichiometric vector γ and can differ from the kinetic equilibrium condition if the detailed balance is not assumed. Possibility of such a difference in general kinetics is sometimes called the 'Wegscheider paradox' [35] to celebrate the work of Wegscheider [36].
Let us consider the isomerisation reaction A 1 A 2 . The corresponding stoichiometric vector is γ = (−1, 1). For this vector, there is one stoichiometric conservation law c 1 + c 2 = b, where b is a positive constant. Equation (23) for this vector has the form: The root of this polynomial is and the point of partial equilibrium is Let us consider the reaction of dissociation A 1 2A 2 . The corresponding stoichiometric vector is γ = (−1, 2). For this vector, there is one stoichiometric conservation law 2c 1 + c 2 = b, where b is a positive constant. Equation (23) for this vector has the form The roots of this polynomial are The sign of the root can be determined from the condition of non-negativity of concentrations. The point of partial equilibrium is Let us consider the reaction A 1 + A 2 A 3 . The corresponding stoichiometric vector is γ = (−1, −1, 1). For this vector, there are two stoichiometric conservation laws, c 2 − c 1 = b 1 and c 1 + c 3 = b 2 , where b 1 and b 2 are positive constants. Equation (23) for this vector has the form The point of partial equilibrium is where k = c The point of partial equilibrium is where k = 4 c eq 1 c eq 2 c eq 3 2 − 1.
For the reaction A 1 + A 2 A 3 + A 4 we have a stoichiometric vector γ = (−1, −1, 1, 1). For this vector, there are three stoichiometric conservation laws The point of partial equilibrium is where An analytical representation of partial equilibria can also be found for many other reactions. In this subsection, we have presented only all the reactions that are used in case studies.

Linear Kinetics
Let us consider the isomerisation cycle There is one conservation law for this system: The line of partial equilibrium for each of the three stoichiometric vectors is defined by (24). For example, for the first reaction, this partial equilibrium line is The lines of partial equilibrium and partial equilibrium points for a given point c are presented in Figure 1. The level sets for Boltzmann's H function and Gorban's H Γ function are presented in Figure 2. It is important to emphasise that these level sets are independent of kinetic constants and are completely determined by the equilibrium for Boltzmann's H function and by the equilibrium and set of stoichiometric vectors Γ for Gorban's H Γ function.  The kinetic equations for the system (29) are: For system (29) with detailed balance the conditions for the reaction rate constants are The system can be completely parametrised by three equilibrium concentrations c eq i and three reaction rate constants, for example, by the constants k + 1 , k + 2 , k + 3 . To obtain the complex balance condition it is necessary to list all the different stoichiometric vectors α ρ and β ρ : The conditions of complex balance are We can see that the complex balance conditions for this system are equivalent to the condition of stationarity of the point c eq and are not equivalent to the detailed balance conditions. The first two equations in (31) are linearly independent, but the third equation is linearly dependent on the first two equations because the sum of these three equations is equivalent to the trivial equality 0 = 0. As a result, this system can be parametrised by three equilibrium concentrations c eq i and four reaction rate constants, for example, by the constants k + 1 , k + 2 , k + 3 , k + −3 . This means that system with complex balance has one additional degree of freedom. Since system (29) can have a complex balance equilibrium, which is not a point of detailed balance, a set of parameters with a stable focus in equilibrium instead of a stable node is not a priori forbidden. To illustrate the possible behaviour of system (29), we selected the parameters presented in Table 1. The results of simulation of system (30) with the parameters listed in Table 1 are partially presented in Figure 3. All other figures can be found online in [37]. For a system with detailed balance, the reaction rate constants k + 1 , k + 2 , k + 3 presented in the Table 1 were used. We can see the different behaviour of the two H functions. For the system with detailed balance, equal equilibrium concentrations and equal reaction rate constants of direct reactions (Set S1.2, Figure 3a) there is no apparent difference between H and H Γ and H Γ = H γ 1 = H γ 3 all the time. A system with a set of parameters S1.1 demonstrates the difference between H and H Γ : there is the switch from H γ 2 to H γ 3 (see Figure 3b). Figure 3c also demonstrates the difference between H and H Γ and the switch from H γ 1 to H γ 2 . Figure 3c demonstrates weak nonmonotonicity of H γ 3 near the time of 4 s. Figure 3d presents a system with detailed balance and a set of parameters S3.2 and demonstrates fast movement from the initial point to the patrial equilibrium of the third reaction (approximately 0.9 s, defined by switch from H γ 2 to H γ 3 ) and then slowly tends to equilibrium.
It can be concluded that the simplest linear isomerisation cycle demonstrates the coincidence of behaviour of H and H Γ for certain set of parameters (see Figure 3a) and the differences between these two Lyapunov functions for other parameters. In the case when the equilibrium is a stable focus (see Figure 3c) there are an infinite number of switches between H γ i , but the high rate of convergence does not allow this effect to be graphically illustrated. In the case of a see a stable node as equilibrium (see Figure 3b,d) we can observe only a finite number (usually one or two) of switches.   Figure 3. The left column presents the trajectories of the system (30) in the phase plane, the middle column contains graphs of H γ i and H Γ versus time, and the right column depicts graphs of Boltzmann's H and Gorban's H Γ versus time. Each row present system with one set of parameters: (a) Set S1.2 with detailed balance, (b) Set S1.1 without detailed balance, (c) Set S1.2 without detailed balance, and (d) Set S3.2 with detailed balance.

Nonlinear Isomerisation Reaction
Let us consider isomerisation reaction There is one conservation law for this system: The lines of partial equilibrium for the first two stoichiometric vectors are defined by (24). For example, for the first reaction, this partial equilibrium line is For the last reaction the partial equilibrium is defined by (27): The lines of partial equilibrium and partial equilibrium points for a given point c are presented in Figure 4. The level sets for Boltzmann's H function and Gorban's H Γ function are presented in Figure 5. It is important to emphasise that these level sets are independent of kinetic constants and are completely determined by the equilibrium for Boltzmann's H function (the same level sets in Figures 2a-c and 5a-c) and by the equilibrium and set of stoichiometric vectors Γ for Gorban's H Γ function (different level sets in Figures 2d-f and 5d-f).  The kinetic equations for the system (32) are: For system (32) with detailed balance the conditions for the reaction rate constants are: The system can be completely parametrised by three equilibrium concentrations c eq i and three reaction rate constants, for example, by the constants k + 1 , k + 2 , k + 3 . To obtain the complex balance condition it is necessary to list all the different stoichiometric vectors α ρ and β ρ : The conditions of complex balance are We can see that the first, third and fourth complex balance conditions for this system are equivalent to detailed balance conditions. This means that for system (33) the detailed and complex balances are the same. For simulation of system (33), we selected three equilibria and four sets of reaction rate constants, presented in Table 2. Table 2. Set of equilibrium concentrations and set of reaction rate constants for simulation of system (33).  Part of the simulation results of system (33) with the parameters listed in Table 2 are presented in Figure 6. All other figures can be found in [37].   We can see the different behaviour of the two H functions. Figure 6a presents the results for a system with equal equilibrium concentrations and equal reaction rate constants of direct reactions. In contrast to the behavior of the linear system (Figure 3a), there is a difference between H and H Γ and H Γ switches from H γ 2 to H γ 1 . All three models in Figure 6 demonstrate the non-monotonicity of H γ 3 and the difference between H and H Γ . The system in Figure 6b demonstrates the switch H Γ from H γ 2 to H γ 1 . The system in Figure 6c demonstrates the switch H Γ from H γ 2 to H γ 1 and then to H γ 3 . We also see that in this case the trajectory intersects the partial equilibrium of the first reaction and then is attracted back to this partial equilibrium. Opposite to system in Figure 3c the equilibrium of this system is a stable node but not a stable focus.
We can conclude that the nonlinear isomerisation reaction demonstrates the difference in the behaviour of H and H Γ for almost all set of parameters. Since the complex balance condition for this system is always is equivalent to the detailed balance condition, the equilibrium point always is a stable node and the number of switches between H γ i is finite and usually equal to one.

Water Gas Shift Reaction
In this subsection, we consider the famous Water Gas Shift reaction (WGS) [38]. More precisely, we consider the redox mechanism proposed by [39] and described in details in [40][41][42]. In the first part of this subsection, we consider the WGS reaction with arbitrary chosen kinetic parameters. To avoid confusion, we call this reaction 'abstract WGS'. In the last part of this subsection we consider the real WGS reaction with all the parameters defined for this reaction. The redox mechanism includes six substances: H 2 O, H 2 , CO, CO 2 , red, Ox. For the abstract WGS model, we use the following substances: There are two reactions in the WGS mechanism: The abstract WGS mechanism include the following reactions: Systems (35) and (36) have four stoichiometric conservation laws: For the WGS reaction these conservations laws mean the conservation of hydrogen, carbon, oxygen and catalyst (accelerator). For the abstract WGS reaction we use the same names of conservation laws. For the simulation we choose the following balance values: hydrogen balance b H = 1, carbon balance b C = 1, oxygen balance b O = b H + b C = 2, and catalyst balance b A = 0.5. These values of balances correspond to one of the standard modes of WGS reaction [40]: "1:1 molar feed ratio [H 2 O/CO]" without hydrogen and carbon dioxide in the initial composition. The line of partial equilibrium for both stoichiometric vectors is defined by (28). For example, for the first reaction, the line of partial equilibrium is The lines of partial equilibrium and the level sets for Boltzmann's H function and Gorban's H Γ function are shown in Figure 7. It is important to emphasise that these level sets are independent of kinetic constants and are completely determined by the equilibrium for Boltzmann's H function and by the equilibrium and set of stoichiometric vectors Γ for Gorban's H Γ function.  The kinetic equations for the system (36) are: For system (36) with detailed balance, the conditions for the reaction rate constants are: The system can be completely parametrised by six equilibrium concentrations c eq i and two reaction rate constants, for example, by the constants k + 1 , k + 2 . To obtain the complex balance condition it is necessary to list all the different stoichiometric vectors α ρ and β ρ : (1, 0, 0, 0, 1, 0), 1, 0, 0, 0, 1), There are two pairs of identical equalities: the first equality coincides with the second one, and the third equality coincides with the fourth one. Moreover, the first and the third equalities are equivalent to the detailed balance conditions. This means that there is no difference between the detailed and complex balance conditions for system (38). For simulation, we use equilibrium concentrations c eq = (0.25, 0.25, 0.5, 0.5, 0.25, 0.25) and reaction rate constants of direct reactions k + = (1, 1). The simulation results are presented in Figure 8. This figure clearly shows the difference between H and H Γ functions and the switching from the H Γ = H γ 2 to H Γ = H γ 1 during dynamics. Now we consider the real WGS reaction with the list of substances H 2 O, H 2 , CO, CO 2 , red, Ox and reactions (35). All conservation laws are the same as for the abstract model, since these two models have the same structure. The coincidence of the complex balance condition with the detailed balance condition also takes place for WGS reaction. The WGS reaction parameters were found for the condition described in [40] "a 1:1 molar feed ratio [H 2 O/CO] and 220 • C the conversion reaches 70%. The equilibrium conversion for these conditions is calculated as 87%". Additional parameters of reactor are described in [40]: ". . . catalyst loading: 1.0 g; . . . GHSV: 6100 h −1 . Size of reactor is 1/2 inch in diameter and 12 inch long." From this information we can identify required values. Let us consider the case b H = 1. Then from the equality of concentrations of H 2 O and CO and absence of all other gases in the original composition we can find The time of a gas movement trough the reactor can be calculated as From the conversion 70% we can require c 3 (t r ) = 0.3b C . From the equilibrium concentration of CO we can find c eq 3 = 0.13b C . From known values of b H , b C , b O , c eq 3 and degree of conversion at time t r we can find b A , c eq 5 , k + 1 , k + 2 by solving optimisation problem The found parameters used in simulation are: the reaction rate constants of direct reactions k + = (80. 53, 146.31) and the equilibrium point c eq = (0.0073, 0.9927, 0.13, 0.87, 0.0015, 0.1227). The lines of partial equilibrium and the level sets for Boltzmann's H-function and Gorban's H Γ function are presented in Figure 9. The results of simulation are presented in Figure 10. We can see that for the real WGS reaction, the equilibrium is very close to the boundary. As a result, the line of partial equilibrium of the first reaction also almost coincides with two sides of boundary of the reaction polygon. The trajectory very quickly achieved the vicinity of the partial equilibrium line of the first reaction and then moved along this line to equilibrium. The time of achieving of the vicinity of the partial equilibrium line of the first reaction could be easily evaluated by switching H Γ from H γ 2 to H γ 1 and was approximately 3 microseconds. It was a very short time compared to 0.59 seconds of the total process time in the reactor. The difference between H and H Γ is obvious for a very short initial time interval. By the way, the behaviour of the abstract system (36) qualitatively coincides with the behaviour of the real WGS system (35).

Hydrogen Chloride Reaction
In this subsection we consider the reaction of hydrogen chloride (HCl) production [43][44][45]. This reaction mechanism includes five substances H 2 In the first part of this subsection we consider reactions with arbitrary chosen reaction rate constants. To avoid confusion we call this reaction 'abstract HCl reaction'. For this reaction we used the following substances A 1 , A 2 , A 3 , A 4 , A 5 and reactions There are two conservation laws in the mechanism (41): 2c 1 + c 2 + c 5 = b H is the hydrogen conservation law and 2c 3 + c 4 + c 5 = b Cl is the chlorine conservation law. This means that there are only three independent variables in the system (41). For this study we selected the variables A 1 , A 3 , A 5 (H 2 , Cl 2 , HCl for system (40)) as independent and all the figures are presented in this space. The reaction polyhedron for this system can be found from the condition that all concentrations are nonnegative. The partial equilibrium surfaces of the first two reactions are defined by (25). For example, for the first reaction, the partial equilibrium is where k = c eq 2 ) 2 c eq 1 .
For the last two reactions, the surfaces of partial equilibrium are defined by (28). For example, for the third reaction, this surface i where surfaces for system (41) with the kinetic curve (the trajectory) are presented for the equilibrium c eq = (0.2, 0.2, 0.25, 0.1, 0.4) and the reaction rate constants of direct reactions k + = (5, 10, 2, 1) in Figure 12. It should be emphasised that the surfaces of partial equilibrium for the first two reactions only look like planes, but in fact they have square root type nonlinearity (see (42) for the surface of partial equilibrium of the first reaction). At the beginning of motion, the trajectory quickly (at about 0.03 s) achieved the partial equilibrium surface of the second reaction, then along this surface the trajectory reached (at about 0.18 s) the intersection of the partial equilibrium surfaces of the first two reactions and then moved along this intersection to the equilibrium (approximately 15 s for the tolerance level 0.0001).
(a) (b) Figure 11. The level sets for system (41) The graphs of H and H Γ are presented in Figure 13. There is a difference between H and H Γ at the initial stage of the reaction (during the first 0.02 s from the approximately 15 s of the full process). We also can observe two switches of H Γ : from H γ 4 to H γ 3 at the first few microseconds and then from H γ 3 to H γ 2 in about 5 milliseconds after the start of the process.  Now we consider the real HCl reaction (40). For the simulation, we used the information from [43][44][45]: the equilibrium point was c eq = (0.198, 0.004, 0.1995, 0.001, 0.6) and the reaction rate constants of the direct reactions were k + = (10 16 , 10 16 , 1.7 × 10 11 , 1.59 × 10 8 ). This reaction system is very stiff and the equilibrium point is almost on the edge between the vertices (0, 0, 0, 0, 1) and (0.5, 0, 0.5, 0, 0). This means that the graphs of the level sets are uninformative and we omit them. Images of the level sets for system (40) can be found in [37]. The partial equilibrium surfaces for this system are presented in Figure 14. It should be emphasized that the partial equilibrium surfaces for the first two reactions only look like planes, but actually have a square root nonlinearity (see (42) for the partial equilibrium surface of the first reaction).  Figure 14. Partial equilibrium surfaces and the trajectory for system (40).
The graphs of H and H Γ are presented in Figure 15. It can be seen that there is a difference between H and H Γ at the initial stage of reaction (during the first 10 −19 s from the approximately 10 −9 s of full process). There is one switch of H Γ : from H γ 4 to H γ 2 approximately at the time moment 3.6 × 10 −21 s.

Conclusions and Outlook
For each reaction mechanism, there exists an infinite family of Gorban's conditionally universal Lyapunov functions H Γ (22) indexed by a finite set of n-dimensional vectors Γ, which should include all the stoichiometric vectors of the elementary reactions but may also include arbitrary vectors with at least one positive and at least one negative element. In all the cases, the level sets for H Γ were found significantly different from the level sets of the classical thermodynamic Lyapunov function H(N) (2) (see Figures 2,5,7 and 11).
The comparison of time dependences of H Γ and H along kinetic trajectories gave more tricky results (Figures 3, 6, 8, 10, 13 and 15). Of course, both functions decreased in time. Their values and the rates of descent were different if all elementary reactions were far from their partial equilibria, but if at least one reaction with the stoichiometric vector γ approached closely its partial equilibrium then H Γ (c) ≈ H γ (c) ≈ H(c) and the difference vanished. Nevertheless, if the kinetic trajectory leaved the small vicinity of the partial equilibrium, then the dynamics of H Γ (c) and H(c) became again different.
The new family of the conditionally universal Lyapunov functions gives the answer to an intriguing question about existence of such functions for non-linear reaction mechanisms (for linear reactions, the answer was done by Rényi [9] and elaborated further by several authors [10][11][12][13]). In addition to this theoretical value, we can expect some new fields of applications for these functions.
There may be many applications of the new conditionally universal Lyapunov functions. We can compare this situation to applications of many different divergences in the applied statistical inference problem [46,47]. Moreover, it is possible to use families of different entropies together and find Maximal Entropy (MaxEnt) sets of distributions instead of single distributions. This is the so-called Maximum of All Entropies (MaxAllEnt) approach that takes into account uncertainty in selection of the measure of uncertainty in the inference problem [48].
Another application is the evaluation of the attainability regions. Each Lyapunov function can serve as a tool for evaluation (from above) the region attainable for kinetic curves because the value of this function should decrease in time [49][50][51].
There exists an obvious necessary condition of attainability of a state y from the state x, H(x) ≥ H(y), but it is not sufficient for attainability by a continuous path, along which H decreases monotonically. For example, a 1D system (with n components and n − 1 conservation laws) cannot come from a state x to a state y if they are on the opposite sides of the equilibrium even if H(x) > H(y). Detailed analysis of attainability in several dimensions led to a beautiful chapter of computational convex combinatorial geometry (for more detailed review we refer to [51]). These results and their generalisations are proved to be useful in optimisation of chemical reactors and related problems [33,52,53].
There remain also some problems. It was mentioned that all the H Γ should have an equivalent f -divergence form (4) (possibly, after a monotonic transformation) and this form is still unknown [8]. From the application perspectives, the following question seems to be even more important: are there other families of the conditionally universal Lyapunov functions for non-linear reaction mechanisms? For linear mechanisms, such a question is fully resolved: any conditionally universal Lyapunov function for linear kinetics has the form of f -divergence (or can be produced from an f -divergence by a monotonic transformation) [12,13,54].