###### Abstract

We study chemical reactions with complex mechanisms under two assumptions: (i)
intermediates are present in small amounts (this is the quasi-steady-state
hypothesis or QSS) and (ii) they are in equilibrium relations with substrates
(this is the quasiequilibrium hypothesis or QE). Under these assumptions, we
prove the generalized mass action law together with the basic relations
between kinetic factors, which are sufficient for the positivity of the
entropy production but hold even without microreversibility, when the detailed
balance is not applicable. Even though QE and QSS produce useful
approximations by themselves, only the combination of these assumptions can
render the possibility beyond the “rarefied gas” limit or the “molecular
chaos” hypotheses. We do not use any a priori form of the kinetic law for the
chemical reactions and describe their equilibria by thermodynamic relations.
The transformations of the intermediate compounds can be described by the
Markov kinetics because of their low density (low density of elementary
events). This combination of assumptions was introduced by Michaelis and
Menten in 1913\. In 1952, Stueckelberg used the same assumptions for the gas
kinetics and produced the remarkable semi-detailed balance relations between
collision rates in the Boltzmann equation that are weaker than the detailed
balance conditions but are still sufficient for the Boltzmann $H$-theorem to
be valid. Our results are obtained within the Michaelis-Menten-Stueckelbeg
conceptual framework.

###### keywords:

chemical kinetics; Lyapunov function; entropy; quasiequilibrium; detailed
balance; complex balance

10.3390/e13050966 13 Received: 25 January 2011; in revised form: 28 March 2011
/ Accepted: 12 May 2011 / Published: 20 May 2011 The Michaelis-Menten-
Stueckelberg Theorem Alexander N. Gorban 1,⋆ and Muhammad Shahzad 1,2 E-Mail:
<EMAIL_ADDRESS>05.70.Ln; 82.20.Db

##  1 Introduction

###  1.1 Main Asymptotic Ideas in Chemical Kinetics

There are several essentially different approaches to asymptotic and scale
separation in kinetics, and each of them has its own area of applicability.

In chemical kinetics various fundamental ideas about asymptotical analysis
were developed Klonowski1983 : Quasieqiulibrium asymptotic (QE), quasi steady-
state asymptotic (QSS), lumping, and the idea of limiting step.

Most of the works on nonequilibrium thermodynamics deal with the QE
approximations and corrections to them, or with applications of these
approximations (with or without corrections). There are two basic formulation
of the QE approximation: The thermodynamic approach, based on entropy maximum,
or the kinetic formulation, based on selection of fast reversible reactions.
The very first use of the entropy maximization dates back to the classical
work of Gibbs Gibb , but it was first claimed for a principle of informational
statistical thermodynamics by Jaynes Jaynes1 . A very general discussion of
the maximum entropy principle with applications to dissipative kinetics is
given in the review Bal . Corrections of QE approximation with applications to
physical and chemical kinetics were developed GKIOeNONNEWT2001 ; GorKar .

QSS was proposed by Bodenstein in 1913 Bodenstein1913 , and the important
Michaelis and Menten work MichaelisMenten1913  was published simultaneously.
It appears that no kinetic theory of catalysis is possible without QSS. This
method was elaborated into an important tool for the analysis of chemical
reaction mechanism and kinetics Semenov1939 ; Christiansen1953 ;
Helfferich1989 . The classical QSS is based on the relative smallness of
concentrations of some of the “active” reagents (radicals, substrate-enzyme
complexes or active components on the catalyst surface) BriggsHaldane1925 ;
Aris1965 ; Segel89 .

Lumping analysis aims to combine reagents into “quasicomponents” for dimension
reduction LumpWei1 ; LumpWei2 ; LumpLiRab1 ; LumpLiRab2 . Wei and Prater Wei62
demonstrated that for (pseudo)monomolecular systems there exist linear
combinations of concentrations which evolve in time independently. These
linear combinations (quasicomponents) correspond to the left eigenvectors of
the kinetic matrix: If $lK=\lambda l$ then
${\mathrm{d}}(l,c)/{\mathrm{d}}t=(l,c)\lambda$, where the standard inner
product $(l,c)$ is the concentration of a quasicomponent. They also
demonstrated how to find these quasicomponents in a properly organized
experiment.

This observation gave rise to a question: How to lump components into proper
quasicomponents to guarantee the autonomous dynamics of the quasicomponents
with appropriate accuracy? Wei and Kuo studied conditions for exact LumpWei1
and approximate LumpWei2  lumping in monomolecular and pseudomonomolecular
systems. They demonstrated that under certain conditions a large monomolecular
system could be well-modelled by a lower-order system.

More recently, sensitivity analysis and Lie group approach were applied to
lumping analysis LumpLiRab1 ; LumpLiRab2 , and more general nonlinear forms of
lumped concentrations were used (for example, concentration of quasicomponents
could be a rational function of $c$).

Lumping analysis was placed in the linear systems theory and the relationships
between lumpability and the concepts of observability, controllability and
minimal realization were demonstrated LumpingOBservability . The lumping
procedures were considered also as efficient techniques leading to nonstiff
systems and demonstrated the efficiency of the developed algorithm on kinetic
models of atmospheric chemistry NonstiffAtmospheric2002 . An optimal lumping
problem can be formulated in the framework of a mixed integer nonlinear
programming (MINLP) and can be efficiently solved with a stochastic
optimization method OptimalLumping2008 .

The concept of limiting step gives the limit simplification: The whole network
behaves as a single step. This is the most popular approach for model
simplification in chemical kinetics and in many areas beyond kinetics. In the
form of a bottleneck approach this approximation is very popular from traffic
management to computer programming and communication networks. Recently, the
concept of the limiting step has been extended to the asymptotology of
multiscale reaction networks GorbaRadul2008 ; GorRadZin2010 .

In this paper, we focus on the combination of the QE approximation with the
QSS approach.

###  1.2 The Structure of the Paper

Almost thirty years ago one of us published a book G1  with Chapter 3 entitled
“Quasiequilibrium and Entropy Maximum”. A research program was formulated
there, and now we are in the position to analyze the achievements of these
three decades and formulate the main results, both theoretical and applied,
and the unsolved problems. In this paper, we start this work and combine a
presentation of theory and application of the QE approximation in physical and
chemical kinetics with exposition of some new results.

We start from the formal description of the general idea of QE and its
possible extensions. In Section 2, we briefly introduce main notations and
some general formulas for exclusion of fast variables by the QE approximation.

In Section 3, we present the history of the QE and the classical confusion
between the QE and the quasi steady state (QSS) approximation. Another
surprising confusion is that the famous Michaelis-Menten kinetics was not
proposed by Michaelis and Menten in 1913 MichaelisMenten1913  but by Briggs
and Haldane BriggsHaldane1925  in 1925. It is more important that Michaelis
and Menten proposed another approximation that is very useful in general
theoretical constructions. We described this approximation for general kinetic
systems. Roughly speaking, this approximation states that any reaction goes
through transformation of fast intermediate complexes (compounds), which (i)
are in equilibrium with the input reagents and (ii) exist in a very small
amount.

One of the most important benefits from this approach is the exclusion of
nonlinear kinetic laws and reaction rate constants for nonlinear reactions.
The nonlinear reactions transform into the reactions of the compounds
production. They are in a fast equilibrium and the equilibrium is ruled by
thermodynamics. For example, when Michaelis and Menten discuss the production
of the enzyme-substrate complex ES from enzyme E and substrate S, they do not
discuss reaction rates. These rates may be unknown. They just assume that the
reaction $E+S\rightleftharpoons ES$ is in equilibrium. Briggs and Haldane
involved this reaction into the kinetic model. Their approach is more general
than the Michaelis–Menten approximation but for the Briggs and Haldane model
we need more information, not only the equilibrium of the reaction
$E+S\rightleftharpoons ES$ but also its rates and constants.

When compounds undergo transformations in a linear first order kinetics, there
is no need to include interactions between them because they are present in
very small amounts in the same volume, and their concentrations are also
small. (By the way, this argument is not applicable to the heterogeneous
catalytic reactions. Although the intermediates are in both small amounts and
in a small volume, i.e., in the surface layer, the concentration of the
intermediates is not small, and their interaction does not vanish when their
amount decreases Yablonskii1991 . Therefore, kinetics of intermediates in
heterogeneous catalysis may be nonlinear and demonstrate bifurcations,
oscillations and other complex behavior.)

In 1952, Stueckelberg Stueckelberg1952  used similar approach in his seminal
paper “$H$-theorem and unitarity of the $S$-matrix”. He studied elastic
collisions of particles as the quasi-chemical reactions

$\mathbf{v+w\to v^{\prime}+w^{\prime}}$

($\mathbf{v,w,v^{\prime},w^{\prime}}$ are velocities of particles) and
demonstrated that for the Boltzmann equation the linear Markov kinetics of the
intermediate compounds results in the special relations for the kinetic
coefficients. These relations are sufficient for the $H$-theorem, which was
originally proved by Boltzmann under the stronger assumption of reversibility
of collisions Boltzmann .

First, the idea of such relations was proposed by Boltzmann as an answer to
the Lorentz objections against Boltzmann’s proof of the $H$-theorem. Lorentz
stated the nonexistence of inverse collisions for polyatomic molecules.
Boltzmann did not object to this argument but proposed the “cyclic balance”
condition, which means balancing in cycles of transitions between states
$S_{1}\to S_{2}\to\ldots\to S_{n}\to S_{1}$. Almost 100 years later,
Cercignani and Lampis CercignaniLamp1981  demonstrated that the Lorenz
arguments are wrong and the new Boltzmann relations are not needed for the
polyatomic molecules under the microreversibility conditions. The detailed
balance conditions should hold.

Nevertheless, Boltzmann’s idea is very seminal. It was studied further by
Heitler Heitler1944  and Coester Coester1951  and the results are sometimes
cited as the “Heitler-Coestler theorem of semi-detailed balance”. In 1952,
Stueckelberg Stueckelberg1952  proved these conditions for the Boltzmann
equation. For the micro-description he used the $S$-matrix representation,
which is in this case equivalent for the Markov microkinetics (see also
Watanabe1955 ).

Later, these relations for the chemical mass action kinetics were rediscovered
and called the complex balance conditions HornJackson1972 ; Feinberg1972 . We
generalize the Michaelis-Menten-Stueckelberg approach and study in Section 5
the general kinetics with fast intermediates present in small amount. In
Subsection 5.7 the big Michaelis-Menten-Stueckelberg theorem is formulated as
the overall result of the previous analysis.

Before this general theory, we introduce the formalism of the QE approximation
with all the necessary notations and examples for chemical kinetics in Section
4.

The result of the general kinetics of systems with intermediate compounds can
be used wider than this specific model of an elementary reaction: The
intermediate complexes with fast equilibria and the Markov kinetics can be
considered as the “construction staging” for general kinetics. In Section 6,
we delete the construction staging and start from the general forms of the
obtained kinetic equations as from the basic laws. We study the relations
between the general kinetic law and the thermodynamic condition of the
positivity of the entropy production.

Sometimes the kinetics equations may not respect thermodynamics from the
beginning. To repair this discrepancy, deformation of the entropy may help. In
Section 7, we show when is it possible to deform the entropy by adding a
linear function to provide agreement between given kinetic equations and the
deformed thermodynamics. As a particular case, we got the “deficiency zero
theorem”.

The classical formulation of the principle of detailed balance deals not with
the thermodynamic and global forms we use but just with equilibria: In
equilibrium each process must be equilibrated with its reverse process. In
Section 7, we demonstrate also that for the general kinetic law the existence
of a point of detailed balance is equivalent to the existence of such a linear
deformation of the entropy that the global detailed balance conditions
(Equation (87) below) hold. Analogously, the existence of a point of complex
balance is equivalent to the global condition of complex balance after some
linear deformation of the entropy.

###  1.3 Main Results: One Asymptotic and Two Theorems

Let us follow the ideas of Michaelis-Menten and Stueckelberg and introduce the
asymptotic theory of reaction rates. Let the list of the components $A_{i}$ be
given. The mechanism of reaction is the list of the elementary reactions
represented by their stoichiometric equations:

$\sum_{i}\alpha_{\rho i}A_{i}\to\sum_{i}\beta_{\rho i}A_{i}\,$ (1)

The linear combinations $\sum_{i}\alpha_{\rho i}A_{i}$ and
$\sum_{i}\beta_{\rho i}A_{i}$ are the complexes. For each complex
$\sum_{i}y_{ji}A_{i}$ from the reaction mechanism we introduce an intermediate
auxiliary state, a compound $B_{j}$. Each elementary reaction is represented
in the form of the “$2n$-tail scheme” (Figure 1) with two intermediate
compounds:

$\sum_{i}\alpha_{\rho i}A_{i}\rightleftharpoons B_{\rho}^{-}\to
B_{\rho}^{+}\rightleftharpoons\sum_{i}\beta_{\rho i}A_{i}\,\vspace{-6pt}$ (2)

Figure 1: A $2n$-tail scheme of an extended elementary reaction.

There are two main assumptions in the Michaelis-Menten-Stueckelberg
asymptotic:

  * •

The compounds are in fast equilibrium with the corresponding input reagents
(QE);

  * •

They exist in very small concentrations compared to other components (QSS).

The smallness of the concentration of the compounds implies that they (i) have
the perfect thermodynamic functions (entropy, internal energy and free energy)
and (ii) the rates of the reactions $B_{i}\to B_{j}$ are linear functions of
their concentrations.

One of the most important benefits from this approach is the exclusion of the
nonlinear reaction kinetics: They are in fast equilibrium and equilibrium is
ruled by thermodynamics.

Under the given smallness assumptions, the reaction rates $r_{\rho}$ for the
elementary reactions have a special form of the generalized mass action law
(see Equation (74) below):

$r_{\rho}=\varphi_{\rho}\exp(\alpha_{\rho},\check{\mu})\,$

where $\varphi_{\rho}>0$ is the kinetic factor and
$\exp(\alpha_{\rho},\check{\mu})$ is the Boltzmann factor. Here and further in
the text $(\alpha_{\rho},\check{\mu})=\sum_{i}\alpha_{\rho i}\check{\mu}_{i}$
is the standard inner product, $\exp(\ ,\ )$ is the exponential of the value
of the inner product and $\check{\mu}_{i}$ are chemical potentials $\mu$
divided on $RT$.

For the prefect chemical systems, $\check{\mu}_{i}=\ln(c_{i}/c_{i}^{*})$,
where $c_{i}$ is the concentration of $A_{i}$ and $c_{i}^{*}>0$ are the
positive equilibrium concentrations. For different values of the conservation
laws there are different positive equilibria. The positive equilibrium
$c_{i}^{*}$ is one of them and it is not important which one is it. At this
point, $\check{\mu}_{i}=0$, hence, the kinetic factor for the perfect systems
is just the equilibrium value of the rate of the elementary reaction at the
equilibrium point $c^{*}$: $\varphi_{\rho}=r_{\rho}(c^{*})$.

The linear kinetics of the compound reactions $B_{i}\to B_{j}$ implies the
remarkable identity for the reaction rates, the complex balance condition
(Equation (89) below)

$\sum_{\rho}\varphi_{\rho}\exp(\check{\mu},\alpha_{\rho})=\sum_{\rho}\varphi_{\rho}\exp(\check{\mu},\beta_{\rho})\,$

for all admissible values of $\check{\mu}$ and given $\varphi$ which may vary
independently. For other and more convenient forms of this condition see
Equation (91) in Section 6. The complex balance condition is sufficient for
the positivity of the entropy production (for decrease of the free energy
under isothermal isochoric conditions). The general formula for the reaction
rate together with the complex balance conditions and the positivity of the
entropy production form the Michaelis-Menten-Stueckelberg theorem (Section
5.7).

The detailed balance conditions (Equation (87) below),

$\varphi_{\rho}^{+}=\varphi_{\rho}^{-}$

for all ${\rho}$, are more restrictive than the complex balance conditions.
For the perfect systems, the detailed balance condition takes the standard
form: $r_{\rho}^{+}(c^{*})=r_{\rho}^{-}(c^{*})$.

We study also some other, less restrictive sufficient conditions for
accordance between thermodynamics and kinetics. For example, we demonstrate
that the $G$-inequality (Equation (92) below)

$\sum_{\rho}\varphi_{\rho}\exp(\check{\mu},\alpha_{\rho})\geq\sum_{\rho}\varphi_{\rho}\exp(\check{\mu},\beta_{\rho})\,$

is sufficient for the entropy growth and, at the same time, weaker than the
condition of complex balance.

If the reaction rates have the form of the generalized mass action law but do
not satisfy the sufficient condition of the positivity of the entropy
production, the situation may be improved by the deformation of the entropy
via addition of a linear function. Such a deformation is always possible for
the zero deficiency systems. Let $q$ be the number of different complexes in
the reaction mechanism, $d$ be the number of the connected components in the
digraph of the transitions between compounds (vertices are compounds and edges
are reactions). To exclude some degenerated cases a hypothesis of weak
reversibility is accepted: For any two vertices $B_{i}$ and $B_{j}$, the
existence of an oriented path from $B_{i}$ to $B_{j}$ implies the existence of
an oriented path from $B_{j}$ to $B_{i}$.

Deficiency of the system is Feinberg1972

$q-d-{\rm rank}\Gamma\geq 0$

where $\Gamma=(\gamma_{ij})=(\beta_{ij}-\alpha_{ij})$ is the stoichiometric
matrix. If the system has zero deficiency then the entropy production becomes
positive after the deformation of the entropy via addition of a linear
function. The deficiency zero theorem in this form is proved in Section 7.3.

Interrelations between the Michaelis-Menten-Stueckelberg asymptotic and the
transition state theory (which is also referred to as the “activated-complex
theory”, “absolute-rate theory”, and “theory of absolute reaction rates”) are
very intriguing. This theory was developed in 1935 by Eyring Eyring1935  and
by Evans and Polanyi EvansPolanyi1935 .

Basic ideas behind the transition state theory are LaidlerTweedale2007 :

  * •

The activated complexes are in a quasi-equilibrium with the reactant
molecules;

  * •

Rates of the reactions are studied by studying the activated complexes at the
saddle point of a potential energy surface.

The similarity is obvious but in the Michaelis-Menten-Stueckelberg asymptotic
an elementary reaction is represented by a couple of compounds with the Markov
kinetics of transitions between them versus one transition state, which moves
along the “reaction coordinate”, in the transition state theory. This is not
exactly the same approach (for example, the theory of absolute reaction rates
uses the detailed balance conditions and does not produce anything similar to
the complex balance).

Important technical tools for the analysis of the Michaelis-Menten-
Stueckelberg asymptotic are the theorem about preservation of the entropy
production in the QE approximation (Section 2 and Appendix 1) and the Morimoto
$H$-theorem for the Markov chains (Appendix 2).

##  2 QE and Preservation of Entropy Production

In this section we introduce informally the QE approximation and the important
theorem about the preservation of entropy production in this approximation. In
Appendix 1, this approximation and the theorem are presented with more formal
details.

Let us consider a system in a domain $U$ of a real vector space $E$ given by
differential equations

$\frac{{\mathrm{d}}x}{{\mathrm{d}}t}=F(x)$ (3)

The QE approximation for (3) uses two basic entities: Entropy and slow
variables.

Entropy $S$ is an increasing concave Lyapunov function for (3) with non-
degenerated Hessian $\partial^{2}S/\partial x_{i}\partial x_{j}$:

$\frac{{\mathrm{d}}S}{{\mathrm{d}}t}\geq 0\,$ (4)

In this approach, the increase of entropy in time is exploited (the Second Law
in the form (4)).

The slow variables $M$ are defined as some differentiable functions of
variables $x$: $M=m(x)$. Here we assume that these functions are linear. More
general nonlinear theory was developed in GorKarQE2006 ; GorKarProjector2004
with applications to the Boltzmann equation and polymer physics. Selection of
the slow variables implies a hypothesis about separation of fast and slow
motion. The slow variables (almost) do not change during the fast motion.
After some initial time, the fast variables with high accuracy are functions
of the slow variables: We can write $x\approx x^{*}_{M}$.

The QE approximation defines the functions $x^{*}_{M}$ as solutions to the
following MaxEnt optimization problem:

$S(x)\to\max\;\;\mbox{subject to}\;\;m(x)=M\,$ (5)

The reasoning behind this approximation is simple: During the fast initial
layer motion, entropy increases and $M$ almost does not change. Therefore, it
is natural to assume that $x^{*}_{M}$ is close to the solution to the MaxEnt
optimization problem (5). Further $x^{*}_{M}$ denotes a solution to the MaxEnt
problem.

A solution to (5), $x^{*}_{M}$, is the QE state, the set of the QE states
$x^{*}_{M}$, parameterized by the values of the slow variables $M$ is the QE
manifold, the corresponding value of entropy

$S^{*}(M)=S(x^{*}_{M})$ (6)

is the QE entropy and the equation for the slow variables

$\frac{{\mathrm{d}}M}{{\mathrm{d}}t}=m(F(x^{*}_{M}))$ (7)

represents the QE dynamics.

The crucial property of the QE dynamics is the preservation of entropy
production.

Theorem about preservation of entropy production. Let us calculate
${\mathrm{d}}S^{*}(M)/{\mathrm{d}}t$ at point $M$ according to the QE dynamics
(7) and find ${\mathrm{d}}S(x)/{\mathrm{d}}t$ at point $x=x^{*}_{M}$ due to
the initial system (3). The results always coincide:

$\frac{{\mathrm{d}}S^{*}(M)}{{\mathrm{d}}t}=\frac{{\mathrm{d}}S(x)}{{\mathrm{d}}t}\,$
(8)

The left hand side in (8) is computed due to the QE approximation (7) and the
right hand side corresponds to the initial system (3). The sketch of the proof
is given in Appendix 1.

The preservation of the entropy production leads to the preservation of the
type of dynamics: If for the initial system (3) entropy production is non-
negative, ${\mathrm{d}}S/{\mathrm{d}}t\geq 0$, then for the QE approximation
(7) the production of the QE entropy is also non-negative,
${\mathrm{d}}S^{*}/{\mathrm{d}}t\geq 0$.

In addition, if for the initial system
$({{\mathrm{d}}S}/{{\mathrm{d}}t})|_{x}=0$ if and only if $F(x)=0$ then the
same property holds in the QE approximation.

##  3 The Classics and the Classical Confusion

###  3.1 The Asymptotic of Fast Reactions

It is difficult to find who introduced the QE approximation. It was impossible
before the works of Boltzmann and Gibbs, and it became very well known after
the works of Jaynes Jaynes1 .

Chemical kinetics has been a source for model reduction ideas for decades. The
ideas of QE appear there very naturally: Fast reactions go to their
equilibrium and, after that, remain almost equilibrium all the time. The
general formalization of this idea looks as follows. The kinetic equation has
the form

$\frac{{\mathrm{d}}N}{{\mathrm{d}}t}=K_{sl}(N)+\frac{1}{\epsilon}K_{fs}(N)$
(9)

Here $N$ is the vector of composition with components $N_{i}>0$, $K_{sl}$
corresponds to the slow reactions, $K_{fs}$ corresponds to fast reaction and
$\epsilon>0$ is a small number. The system of fast reactions has the linear
conservation laws $b_{l}(N)=\sum_{j}b_{lj}N_{j}$: $b_{l}(K_{fs}(N))\equiv 0$.

The fast subsystem

$\frac{{\mathrm{d}}N}{{\mathrm{d}}t}=K_{fs}(N)$

tends to a stable positive equilibrium $N^{*}$ for any positive initial state
$N(0)$ and this equilibrium is a function of the values of the linear
conservation laws $b_{l}(N(0))$. In the plane $b_{l}(N)=b_{l}(N(0))$ the
equilibrium is asymptotically stable and exponentially attractive.

Vector $b(N)=(b_{l}(N))$ is the vector of slow variables and the QE
approximation is

$\frac{{\mathrm{d}}b}{{\mathrm{d}}t}=b(K_{sl}(N^{*}(b))\,$ (10)

In chemical kinetics, equilibria can be described by conditional entropy
maximum (or conditional extremum of other thermodynamic potentials).
Therefore, for these cases we can apply the formalism of the quasiequilibrium
approximation. The thermodynamic Lyapunov functions serve as tools for
stability analysis and for model reduction Hangos2010 .

The QE approximation, the asymptotic of fast reactions, is well known in
chemical kinetics. Another very important approximation was invented in
chemical kinetics as well. It is the Quasi Steady State (QSS) approximation.
QSS was proposed in Bodenstein1913  and was elaborated into an important tool
for analysis of chemical reaction mechanisms and kinetics Semenov1939 ;
Christiansen1953 ; Helfferich1989 . The classical QSS is based on the relative
smallness of concentrations of some of “active” reagents (radicals, substrate-
enzyme complexes or active components on the catalyst surface) Aris1965 ;
Segel89 . In the enzyme kinetics, its invention was traditionally connected to
the so-called Michaelis-Menten kinetics.

###  3.2 QSS and the Briggs-Haldane Asymptotic

Perhaps the first very clear explanation of the QSS was given by Briggs and
Haldane in 1925 BriggsHaldane1925 . Briggs and Haldane consider the simplest
enzyme reaction $S+E\leftrightharpoons SE\to P+E$ and mention that the total
concentration of enzyme ($[E]+[SE]$) is “negligibly small” compared with the
concentration of substrate $[S]$. After that they conclude that
$\frac{{\mathrm{d}}}{{\mathrm{d}}t}{[SE]}$ is “negligible” compared with
$\frac{{\mathrm{d}}}{{\mathrm{d}}t}{[S]}$ and
$\frac{{\mathrm{d}}}{{\mathrm{d}}t}{[P]}$ and produce the now famous
‘Michaelis-Menten’ formula, which was unknown to Michaelis and Menten:
$k_{1}[E][S]=(k_{-1}+k_{2})[ES]$ or

$[ES]=\frac{[E][S]}{K_{M}+[S]}\;\;{\rm
and}\;\;\frac{{\mathrm{d}}}{{\mathrm{d}}t}{[P]}=k_{2}[ES]=\frac{k_{2}[E][S]}{K_{M}+[S]}\,$
(11)

where the “Michaelis-Menten constant” is

$K_{M}=\frac{k_{-1}+k_{2}}{k_{1}}$

There is plenty of misleading comments in later publications about QSS. Two
most important confusions are:

  * •

Enzymes (or catalysts, or radicals) participate in fast reactions and, hence,
relax faster than substrates or stable components. This is obviously wrong for
many QSS systems: For example, in the Michaelis-Menten system all reactions
include enzyme together with substrate or product. There are no separate fast
reactions for enzyme without substrate or product.

  * •

Concentrations of intermediates are constant because in QSS we equate their
time derivatives to zero. In general case, this is also wrong: We equate the
kinetic expressions for some time derivatives to zero, indeed, but this just
exploits the fact that the time derivatives of intermediates concentrations
are small together with their values, but not obligatory zero. If we accept
QSS then these derivatives are not zero as well: To prove this we can just
differentiate the Michaelis-Menten formula (11) and find that [ES] in QSS is
almost constant when $[S]\gg K_{M}$, this is an additional condition,
different from the Briggs-Haldane condition $[E]+[AE]\ll[S]$ (for more details
see Segel89 ; Klonowski1983 ; Yablonskii1991  and a simple detailed case study
LiShenLi2008 ).

After a simple transformation of variables the QSS smallness of concentration
transforms into a separation of time scales in a standard singular
perturbation form (see, for example Yablonskii1991 ; GorbKarlinChemEngS2003 ).
Let us demonstrate this on the traditional Michaelis-Menten system:

$\begin{split}&\frac{{\mathrm{d}}[S]}{{\mathrm{d}}t}=-k_{1}[S][E]+k_{-1}[SE]\,\\\
&\frac{{\mathrm{d}}[SE]}{{\mathrm{d}}t}=k_{1}[S][E]-(k_{-1}+k_{2})[SE]\,\\\
&[E]+[SE]=e=const,\,[S]+[P]=s=const\,\end{split}$ (12)

This is a homogeneous system with the isochoric (fixed volume) conditions for
which we write the equations. The Briggs-Haldane condition is $e\ll s$. Let us
use dimensionless variables $x=[S]/s$, $\xi=[SE]/e$:

$\begin{split}&\frac{s}{e}\frac{{\mathrm{d}}x}{{\mathrm{d}}t}=-sk_{1}x(1-\xi)+k_{-1}\xi\,\\\
&\frac{{\mathrm{d}}\xi}{{\mathrm{d}}t}=sk_{1}x(1-\xi)-(k_{-1}+k_{2})\xi\,\end{split}$
(13)

To obtain the standard singularly perturbed system with the small parameter at
the derivative, we need to change the time scale. This means that when $e\to
0$ the reaction goes proportionally slower and to study this limit properly we
have to adjust the time scale: ${\mathrm{d}}\tau=\frac{e}{s}{\mathrm{d}}t$:

$\begin{split}&\frac{{\mathrm{d}}x}{{\mathrm{d}}\tau}=-sk_{1}x(1-\xi)+k_{-1}\xi\,\\\
&\frac{e}{s}\frac{{\mathrm{d}}\xi}{{\mathrm{d}}\tau}=sk_{1}x(1-\xi)-(k_{-1}+k_{2})\xi\,\end{split}$
(14)

For small $e/s$, the second equation is a fast subsystem. According to this
fast equation, for a given constant $x$, the variable $\xi$ relaxes to

$\xi_{\rm QSS}=\frac{sx}{K_{M}+sx}$

exponentially, as $\exp(-(sk_{1}x+k_{-1}+k_{2})t)$. Therefore, the classical
singular perturbation theory based on the Tikhonov theorem Tikhonov1952 ;
Vasil'eva1963  can be applied to the system in the form (14) and the QSS
approximation is applicable even on an infinite time interval Hoppensteadt1966
. This transformation of variables and introduction of slow time is a standard
procedure for rigorous proof of QSS validity in catalysis Yablonskii1991 ,
enzyme kinetics Battelu1986  and other areas of kinetics and chemical
engineering Aris1965 .

It is worth to mention that the smallness of parameter $e/s$ can be easily
controlled in experiments, whereas the time derivatives, transformation rates
and many other quantities just appear as a result of kinetics and cannot be
controlled directly.

###  3.3 The Michaelis and Menten Asymptotic

QSS is not QE but the classical work of Michaelis and Menten
MichaelisMenten1913  was done on the intersection of QSS and QE. After the
brilliantly clear work of Briggs and Haldane, the name “Michaelis-Menten” was
attached to the Briggs and Haldane equation and the original work of Michaelis
and Menten was considered as an important particular case of this approach, an
approximation with additional and not necessary assumptions of QE. From our
point of view, the Michaelis-Menten work includes more and may give rise to an
important general class of kinetic models.

Michaelis and Menten studied the “fermentative splitting of cane sugar”. They
introduced three “compounds”: The sucrose-ferment combination, the fructose-
ferment combination and the glucose-ferment combination. The fundamental
assumption of their work was “that the rate of breakdown at any moment is
proportional to the concentration of the sucrose-invertase compound”.

They started from the assumption that at any moment according to the mass
action law

$[S_{i}][E]=K_{i}[S_{i}E]$ (15)

where $[S_{i}]$ is the concentration of the $i$th sugar (here, $i=0$ for
sucrose, 1 for fructose and 2 for glucose), $[E]$ is the concentration of the
free invertase and $K_{i}$ is the $i$th equilibrium constant.

For simplification, they use the assumption that the concentration of any
sugar in question in free state is practically equal to that of the total
sugar in question.

Finally, they obtain

$[S_{0}E]=\frac{e[S_{0}]}{K_{0}(1+q[P])+[S_{0}]}\,$ (16)

where $e=[E]+\sum_{i}[S_{i}E]$, $[P]=[S_{1}]=[S_{2}]$ and
$q=\frac{1}{K_{1}}+\frac{1}{K_{2}}$.

Of course, this formula may be considered as a particular case of the Briggs-
Haldane formula (11) if we take $k_{-1}\gg k_{2}$ in (11) (i.e., the
equilibration $S+E\leftrightharpoons SE$ is much faster than the reaction
$SE\to P+E$) and assume that $q=0$ in (16) (i.e., fructose-ferment combination
and glucose-ferment combination are practically absent).

This is the truth but may be not the complete truth. The Michaelis-Menten
approach with many compounds which are present in small amounts and satisfy
the QE assumption (15) is a seed of the general kinetic theory for perfect and
non-perfect mixtures.

##  4 Chemical Kinetics and QE Approximation

###  4.1 Stoichiometric Algebra and Kinetic Equations

In this section, we introduce the basic notations of the chemical kinetics
formalism. For more details see, for example, Yablonskii1991 .

The list of components is a finite set of symbols $A_{1},\ldots,A_{n}$.

A reaction mechanism is a finite set of the stoichiometric equations of
elementary reactions:

$\sum_{i}\alpha_{\rho i}A_{i}\to\sum_{i}\beta_{\rho i}A_{i}\,$ (17)

where $\rho=1,\ldots,m$ is the reaction number and the stoichiometric
coefficients $\alpha_{\rho i},\beta_{\rho i}$ are nonnegative integers.

A stoichiometric vector $\gamma_{\rho}$ of the reaction (17) is a
$n$-dimensional vector with coordinates

$\gamma_{\rho i}=\beta_{\rho i}-\alpha_{\rho i}\,$ (18)

that is, “gain minus loss” in the $\rho$th elementary reaction.

A nonnegative extensive variable $N_{i}$, the amount of $A_{i}$, corresponds
to each component. We call the vector $N$ with coordinates $N_{i}$ “the
composition vector”. The concentration of $A_{i}$ is an intensive variable
$c_{i}=N_{i}/V$, where $V>0$ is the volume. The vector $c=N/V$ with
coordinates $c_{i}$ is the vector of concentrations.

A non-negative intensive quantity, $r_{\rho}$, the reaction rate, corresponds
to each reaction (17). The kinetic equations in the absence of external fluxes
are

$\frac{{\mathrm{d}}N}{{\mathrm{d}}t}=V\sum_{\rho}r_{\rho}\gamma_{\rho}$ (19)

If the volume is not constant then equations for concentrations include
$\dot{V}$ and have different form (this is typical for the combustion
reactions, for example).

For perfect systems and not so fast reactions, the reaction rates are
functions of concentrations and temperature given by the mass action law for
the dependance on concentrations and by the generalized Arrhenius equation for
the dependance on temperature $T$.

The mass action law states:

$r_{\rho}(c,T)=k_{\rho}(T)\prod_{i}c_{i}^{\alpha_{\rho i}}\,$ (20)

where $k_{\rho}(T)$ is the reaction rate constant.

The generalized Arrhenius equation is:

$k_{\rho}(T)=A_{\rho}\exp\left(\frac{S_{{\rm
a}\rho}}{R}\right)\exp\left(-\frac{E_{{\rm a}\rho}}{RT}\right)\,$ (21)

where $R=8.314\,472~{}\frac{\mathrm{J}}{\mathrm{K~{}mol}}$ is the universal,
or ideal gas constant, $E_{{\rm a}\rho}$ is the activation energy, $S_{{\rm
a}\rho}$ is the activation entropy (i.e., $E_{{\rm a}\rho}-TS_{{\rm a}\rho}$
is the activation free energy), $A_{\rho}$ is the constant pre-exponential
factor. Some authors neglect the $S_{{\rm a}\rho}$ term because it may be less
important than the activation energy, but it is necessary to stress that
without this term it may be impossible to reconcile the kinetic equations with
the classical thermodynamics.

In general, the constants for different reactions are not independent. They
are connected by various conditions that follow from thermodynamics (the
second law, the entropy growth for isolated systems) or microreversibility
assumption (the detailed balance and the Onsager reciprocal relations). In
Section 6.2 we discuss these conditions in more general settings.

For nonideal systems, more general kinetic law is needed. In Section 5 we
produce such a general law following the ideas of the original Michaelis and
Menten paper (this is not the same as the famous “Michaelis-Menten kinetics”).
For this work we need a general formalism of QE approximation for chemical
kinetics.

###  4.2 Formalism of QE Approximation for Chemical Kinetics

####  4.2.1 4.2.1. QE Manifold

In this section, we describe the general formalism of the QE for chemical
kinetics following GorbKarlinChemEngS2003 .

The general construction of the quasi-equilibrium manifold gives the following
procedure. First, let us consider the chemical reactions in a constant volume
under the isothermal conditions. The free energy $F(N,T)=Vf(c,T)$ should
decrease due to reactions. In the space of concentrations, one defines a
subspace of fast motions $L$. It should be spanned by the stoichiometric
vectors of fast reactions.

Slow coordinates are linear functions that annulate $L$. These functions form
a subspace in the space of linear functions on the concentration space.
Dimension of this space is $s=n-\dim L$. It is necessary to choose any basis
in this subspace. We can use for this purpose a basis $b_{j}$ in $L^{\perp}$,
an orthogonal complement to $L$ and define the basic functionals as
$b_{j}(N)=(b_{j},N)$.

The description of the QE manifold is very simple in the Legendre transform.
The chemical potentials are partial derivatives

$\mu_{i}=\frac{\partial F(N,T)}{\partial N_{i}}=\frac{\partial
f(c,T)}{\partial c_{i}}\,$ (22)

Let us use $\mu_{i}$ as new coordinates. In these new coordinates (the
“conjugated coordinates”), the QE manifold is just an orthogonal complement to
$L$. This subspace, $L^{\perp}$, is defined by equations

$\sum_{i}\mu_{i}\gamma_{i}=0\;\;{\rm for\;any}\;\;\gamma\in L$ (23)

It is sufficient to take in (23) not all $\gamma\in L$ but only elements from
a basis in $L$. In this case, we get the system of $n-\dim L$ linear equations
of the form (23) and their solution does not cause any difficulty. For the
actual computations, one requires the inversion from $\mu$ to $c$.

It is worth to mention that the problems of the selection of the slow
variables and of the description of the QE manifold in the conjugated
variables can be considered as the same problem of description of the
orthogonal complement, $L^{\perp}$.

To finalize the construction of the QE approximation, we should find for any
given values of slow variables (and of conservation laws) $b_{i}$ the
corresponding point on the QE manifold. This means that we have to solve the
system of equations for $c$:

$b(N)=b;\;\;(\mu(c,T),\gamma_{\rho})=0\,$ (24)

where $b$ is the vector of slow variables, $\mu$ is the vector of chemical
potentials and vectors $\gamma_{\rho}$ form a basis in $L$. After that, we
have the QE dependence $c_{QE}(b)$ and for any admissible value of $b$ we can
find all the reaction rates and calculate $\dot{b}$.

Unfortunately, the system (24) can be solved analytically only in some special
cases. In general case, we have to solve it numerically. For this purpose, it
may be convenient to keep the optimization statement of the problem:
$F\to\min$ subject to given $b$. There exists plenty of methods of convex
optimization for solution of this problem.

The standard toy example gives us a fast dissociation reaction. Let a
homogeneous reaction mechanism include a fast reaction of the form
$A+B\rightleftharpoons AB$. We can easily find the QE approximation for this
fast reaction. The slow variables are the quantities $b_{1}=N_{A}-N_{B}$ and
$b_{2}=N_{A}+N_{B}+N_{C}$ which do not change in this reaction. Let the
chemical potentials be $\mu_{A}/RT=\ln c_{A}+\mu_{A0}$, $\mu_{B}/RT=\ln
c_{B}+\mu_{B0}$, $\mu_{AB}/RT=\ln c_{AB}+\mu_{AB0}$. This corresponds to the
free energy $F=VRT\sum_{i}c_{i}(\ln c_{i}+\mu_{i0})$, the correspondent free
entropy (the Massieu-Planck potential) is $-F/T$. The stoichiometric vector is
$\gamma=(-1,-1,1)$ and the equations (24) take the form

$c_{A}-c_{B}=\frac{b_{1}}{V}\,,\;\;c_{A}+c_{B}+c_{AB}=\frac{b_{2}}{V}\,,\;\;\frac{c_{AB}}{c_{A}c_{B}}=K\,$

where $K$ is the equilibrium constant $K=\exp(\mu_{A0}+\mu_{B0}-\mu_{AB0})$.

From these equations we get the expressions for the QE concentrations:

$c_{A}(b_{1},b_{2})=\frac{1}{2}\frac{b_{1}}{V}-\frac{1}{K}+\sqrt{\left(\frac{1}{2}\frac{b_{1}}{V}-\frac{1}{K}\right)^{2}+\frac{b_{1}+b_{2}}{KV}}$
$c_{B}(b_{1},b_{2})=c_{A}(b_{1},b_{2})-\frac{b_{1}}{V}\,,\;\;c_{AB}(b_{1},b_{2})=\frac{b_{1}+b_{2}}{V}-2c_{A}(b_{1},b_{2})$

The QE free entropy is the value of the free entropy at this point,
$c(b_{1},b_{2})$.

####  4.2.2 4.2.2. QE in Traditional MM System

Let us return to the simplest homogeneous enzyme reaction
$E+S\rightleftharpoons ES\rightarrow P+S$, the traditional Michaelis-Menten
System (12) (it is simpler than the system studied by Michaelis and Menten
MichaelisMenten1913 ). Let us assume that the reaction $E+S\rightleftharpoons
ES$ is fast. This means that both $k_{1}$ and $k_{-1}$ include large
parameters: $k_{1}=\frac{1}{\epsilon}\kappa_{1}$,
$k_{-1}=\frac{1}{\epsilon}\kappa_{-1}$. For small $\epsilon$, we will apply
the QE approximation. Only three components participate in the fast reaction,
$A_{1}=S$, $A_{2}=E$, $A_{3}=ES$. For analysis of the QE manifold we do not
need to involve other components.

The stoichiometric vector of the fast reaction is $\gamma=(-1,-1,1)$. The
space $L$ is one-dimensional and its basis is this vector $\gamma$. The space
$L^{\perp}$ is two-dimensional and one of the convenient bases is
$b_{1}=(1,0,1)$, $b_{2}=(0,1,1)$. The corresponding slow variables are
$b_{1}(N)=N_{1}+N_{3}$, $b_{2}(N)=N_{2}+N_{3}$. The first slow variable is the
sum of the free substrate and the substrate captured in the enzyme-substrate
complex. The second of them is the conserved quantity, the total amount of
enzyme.

The equation for the QE manifold is (15): $k_{1}c_{1}c_{2}=k_{-1}c_{3}$ or
$\frac{c_{1}}{c_{1}^{*}}\frac{c_{2}}{c_{2}^{*}}=\frac{c_{3}}{c_{3}^{*}}$
because $k_{1}c_{1}^{*}c_{2}^{*}=k_{-1}c_{3}^{*}$, where
$c_{i}^{*}=c_{i}^{*}(T)>0$ are the so-called standard equilibrium values and
for perfect systems $\mu_{i}=RT\ln(c_{i}/c_{i}^{*})$,
$F=RTV\sum_{i}c_{i}(\ln(c_{i}/c_{i}^{*})-1)$.

Let us fix the slow variables and find $c_{1,2,3}$. Equations (24) turn into

$c_{1}+c_{3}=b_{1}\,,\;c_{2}+c_{3}=b_{2}\,,\;k_{1}c_{1}c_{2}=k_{-1}c_{3}$

Here we change dynamic variables from $N$ to $c$ because this is a homogeneous
system with constant volume.

If we use $c_{1}=b_{1}-c_{3}$ and $c_{2}=b_{2}-c_{3}$ then we obtain a
quadratic equation for $c_{3}$:

$k_{1}c_{3}^{2}-(k_{1}b_{1}+k_{1}b_{2}+k_{-1})c_{3}+k_{1}b_{1}b_{2}=0\,$ (25)

Therefore,

$c_{3}(b_{1},b_{2})=\frac{1}{2}\left(b_{1}+b_{2}+\frac{k_{-1}}{k_{1}}\right)-\frac{1}{2}\sqrt{\left(b_{1}+b_{2}+\frac{k_{-1}}{k_{1}}\right)^{2}-4b_{1}b_{2}}$

The sign “$-$” is selected to provide positivity of all $c_{i}$. This choice
provides also the proper asymptotic: $c_{3}\to 0$ if any of $b_{i}\to 0$. For
other $c_{1,2}$ we should use $c_{1}=b_{1}-c_{3}$ and $c_{2}=b_{2}-c_{3}$.

The time derivatives of concentrations are:

$\begin{split}&\dot{c}_{1}=-k_{1}c_{1}c_{2}+k_{-1}c_{3}+v_{\rm in}c_{1}^{\rm
in}-v_{\rm out}c_{1}\,\\\
&\dot{c}_{2}=-k_{1}c_{1}c_{2}+(k_{-1}+k_{2})c_{3}+v_{\rm in}c_{2}^{\rm
in}-v_{\rm out}c_{2}\,\\\
&\dot{c}_{3}=k_{1}c_{1}c_{2}-(k_{-1}+k_{2})c_{3}+v_{\rm in}c_{3}^{\rm
in}-v_{\rm out}c_{3}\,\\\ &\dot{c}_{4}=k_{2}c_{3}+v_{\rm in}c_{4}^{\rm
in}-v_{\rm out}c_{4}\,\end{split}$ (26)

here we added external flux with input and output velocities (per unite
volume) $v_{\rm in}$ and $v_{\rm out}$ and input concentrations $c^{\rm in}$.
This is done to stress that the QE approximation holds also for a system with
fluxes if the fast equilibrium subsystem is fast enough. The input and output
velocities are the same for all components because the system is homogeneous.

The slow system is

$\begin{split}&\dot{b}_{1}=\dot{c}_{1}+\dot{c}_{3}=-k_{2}c_{3}+v_{\rm
in}b_{1}^{\rm in}-v_{\rm out}b_{1}\,\\\
&\dot{b}_{2}=\dot{c}_{2}+\dot{c}_{3}=v_{\rm in}b_{2}^{\rm in}-v_{\rm
out}b_{2}\,\\\ &\dot{c}_{4}=k_{2}c_{3}+v_{\rm in}c_{4}^{\rm in}-v_{\rm
out}c_{4}\,\end{split}$ (27)

where $b_{1}^{\rm in}=c_{1}^{\rm in}+c_{3}^{\rm in}$, $b_{2}^{\rm
in}=c_{2}^{\rm in}+c_{3}^{\rm in}$.

Now, we should use the expression for $c_{3}(b_{1},b_{2})$:

$\begin{split}\dot{b}_{1}=&-k_{2}\frac{1}{2}\left[\left(b_{1}+b_{2}+\frac{k_{-1}}{k_{1}}\right)-\frac{1}{2}\sqrt{\left(b_{1}+b_{2}+\frac{k_{-1}}{k_{1}}\right)^{2}-4b_{1}b_{2}}\,\right]+v_{\rm
in}b_{1}^{\rm in}-v_{\rm out}b_{1}\,\\\
\dot{c}_{4}=&k_{2}\frac{1}{2}\left[\left(b_{1}+b_{2}+\frac{k_{-1}}{k_{1}}\right)-\frac{1}{2}\sqrt{\left(b_{1}+b_{2}+\frac{k_{-1}}{k_{1}}\right)^{2}-4b_{1}b_{2}}\,\right]+v_{\rm
in}c_{4}^{\rm in}-v_{\rm out}c_{4}\,\\\ \dot{b}_{2}=&v_{\rm in}b_{2}^{\rm
in}-v_{\rm out}b_{2}\,\end{split}$ (28)

It is obvious here that in the reduced system (28) there exists one reaction
from the lumped component with concentration $b_{1}$ (the total amount of
substrate in free state and in the substrate-enzyme complex) into the
component (product) with concentration $c_{4}$. The rate of this reaction is
$k_{2}c(b_{1}b_{2})$. The lumped component with concentration $b_{2}$ (the
total amount of the enzyme in free state and in the substrate-enzyme complex)
affects the reaction rate but does not change in the reaction.

Let us use for simplification of this system the assumption of the substrate
excess (we follow the logic of the original Michaelis and Menten paper
MichaelisMenten1913 ):

$[S]\gg[SE]\,,\;\;{\it i.e.},\;\;b_{1}\gg c_{3}\,$ (29)

Under this assumption, the quadratic equation (25) transforms into

$\left(1+\frac{b_{2}}{b_{1}}+\frac{k_{-1}}{k_{1}b_{1}}\right)c_{3}=b_{2}+o\left(\frac{c_{3}}{b_{1}}\right)\,$
(30)

and in this approximation

$c_{3}=\frac{b_{2}b_{1}}{b_{1}+b_{2}+\frac{k_{-1}}{k_{1}}}$ (31)

(compare to (16) and (11): This equation includes an additional term $b_{2}$
in denominator because we did not assume formally anything about the smallness
of $b_{2}$ in (29)).

After this simplification, the QE slow equations (28) take the form

$\begin{split}&\dot{b}_{1}=-\frac{k_{2}b_{2}b_{1}}{b_{1}+b_{2}+\frac{k_{-1}}{k_{1}}}+v_{\rm
in}b_{1}^{\rm in}-v_{\rm out}b_{1}\,\\\ &\dot{b}_{2}=v_{\rm in}b_{2}^{\rm
in}-v_{\rm out}b_{2}\,\\\
&\dot{c}_{4}=\frac{k_{2}b_{2}b_{1}}{b_{1}+b_{2}+\frac{k_{-1}}{k_{1}}}+v_{\rm
in}c_{4}^{\rm in}-v_{\rm out}c_{4}\,\end{split}$ (32)

This is the typical form in the reduced equations for catalytic reactions:
Nominator in the reaction rate corresponds to the “brutto reaction” $S+E\to
P+E$ Yablonskii1991 ; YAbLaz1997 .

####  4.2.3 4.2.3. Heterogeneous Catalytic Reaction

For the second example, let us assume equilibrium with respect to the
adsorption in the CO on Pt oxidation:

CO+Pt$\rightleftharpoons$PtCO; O2+2Pt$\rightleftharpoons$2PtO

(for detailed discussion of the modeling of CO on Pt oxidation, this “Mona
Liza” of catalysis, we address readers to Yablonskii1991 ). The list of
components involved in these 2 reactions is: $A_{1}=$ CO, $A_{2}=$ O2,
$A_{3}=$ Pt, $A_{4}=$ PtO, $A_{5}=$ PtCO (CO2 does not participate in
adsorption and may be excluded at this point).

Subspace $L$ is two-dimensional. It is spanned by the stoichiometric vectors,
$\gamma_{1}=(-1,0,-1,0,1)$, $\gamma_{2}=(0,-1,-2,2,0)$.

The orthogonal complement to $L$ is a three-dimensional subspace spanned by
vectors $(0,2,0,1,0)$, $(1,0,0,0,1)$, $(0,0,1,1,1)$. This basis is not
orthonormal but convenient because of integer coordinates.

The corresponding slow variables are

$\begin{split}&b_{1}=2N_{2}+N_{4}=2N_{\rm O_{2}}+N_{\rm PtO}\,\\\
&b_{2}=N_{1}+N_{5}=N_{\rm CO}+N_{\rm PtCO}\,\\\
&b_{3}=N_{3}+N_{4}+N_{5}=N_{\rm Pt}+N_{\rm PtO}+N_{\rm PtCO}\,\end{split}$
(33)

For heterogeneous systems, caution is needed in transition between $N$ and $c$
variables because there are two “volumes” and we cannot put in (33) $c_{i}$
instead of $N_{i}$: $N_{\rm gas}=V_{\rm gas}c_{\rm gas}$ but $N_{\rm
surf}=V_{\rm surf}c_{\rm surf}$, where where $V_{\rm gas}$ is the volume of
gas, $V_{\rm surf}$ is the area of surface.

There is a law of conservation of the catalyst: $N_{\rm Pt}+N_{\rm PtO}+N_{\rm
PtCO}=b_{3}=const$. Therefore, we have two non-trivial dynamical slow
variables, $b_{1}$ and $b_{2}$. They have a very clear sense: $b_{1}$ is the
amount of atoms of oxygen accumulated in O2 and PtO and $b_{2}$ is the amount
of atoms of carbon accumulated in CO and PtCO.

The free energy for the perfect heterogeneous system has the form

$F=V_{\rm gas}RT\sum_{A_{i}\,{\rm
gas}}c_{i}\left(\ln\left(\frac{c_{i}}{c_{i}^{*}}\right)-1\right)+V_{\rm
surf}RT\sum_{A_{i}\,{\rm
surf}}c_{i}\left(\ln\left(\frac{c_{i}}{c_{i}^{*}}\right)-1\right)\,$ (34)

where $c_{i}$ are the corresponding concentrations and
$c_{i}^{*}=c_{i}^{*}(T)>0$ are the so-called standard equilibrium values. (The
QE free energy is the value of the free energy at the QE point.)

From the expression (34) we get the chemical potentials of the perfect mixture

$\mu_{i}=RT\ln\left(\frac{c_{i}}{c_{i}^{*}}\right)\,$ (35)

The QE manifold in the conjugated variables is given by equations:

$-\mu_{1}-\mu_{3}+\mu_{5}=0\,;\;-\mu_{2}-2\mu_{3}+2\mu_{4}=0$

It is trivial to resolve these equations with respect to $\mu_{3,4}$, for
example:

$\mu_{4}=\frac{1}{2}\mu_{2}+\mu_{3}\,;\;\mu_{5}=\mu_{1}+\mu_{3}$

or with the standard equilibria:

$\frac{c_{4}}{c_{4}^{*}}=\frac{c_{3}}{c_{3}^{*}}\sqrt{\frac{c_{2}}{c_{2}^{*}}}\,,\;\frac{c_{5}}{c_{5}^{*}}=\frac{c_{1}}{c_{1}^{*}}\frac{c_{3}}{c_{3}^{*}}$

or in the kinetic form (we assume that the kinetic constants are in accordance
with thermodynamics and all these forms are equivalent):

$k_{1}c_{1}c_{3}=k_{-1}c_{5}\,,\;k_{2}c_{2}c_{3}^{2}=k_{-2}c_{4}^{2}\,$ (36)

The next task is to solve the system of equations:

$\begin{split}&k_{1}c_{1}c_{3}=k_{-1}c_{5}\,,\;k_{2}c_{2}c_{3}^{2}=k_{-2}c_{4}^{2}\,,2V_{\rm
gas}c_{2}+V_{\rm surf}c_{4}=b_{1}\,,\;\\\ &V_{\rm gas}c_{1}+V_{\rm
surf}c_{5}=b_{2}\,,\;V_{\rm surf}(c_{3}+c_{4}+c_{5})=b_{3}\end{split}$ (37)

This is a system of five equations with respect to five unknown variables,
$c_{1,2,3,4,5}$. We have to solve them and use the solution for calculation of
reaction rates in the QE equations for the slow variables. Let us construct
these equations first, and then return to (37).

We assume the adsorption (the Langmuir-Hinshelwood) mechanism of CO oxidation
(the numbers in parentheses are used below for the numeration of the reaction
rate constants):

$\begin{split}(\pm 1)\;\mbox{CO+Pt$\rightleftharpoons$PtCO}\,\\\ (\pm
2)\;\mbox{O${}_{2}$+2Pt$\rightleftharpoons$2PtO}\,\\\
(3)\;\mbox{PtO+PtCO$\to$CO${}_{2}$+2Pt}\;\end{split}$ (38)

The kinetic equations for this system (including the flux in the gas phase) is

$\displaystyle{\rm CO}\;$ $\displaystyle\dot{N}_{1}=V_{\rm
surf}(-k_{1}c_{1}c_{3}+k_{-1}c_{5})+V_{\rm gas}(v_{\rm in}c_{1}^{\rm
in}-v_{\rm out}c_{1})\,$ $\displaystyle{\rm O_{2}}$
$\displaystyle\dot{N}_{2}=V_{\rm
surf}(-k_{2}c_{2}c_{3}^{2}+k_{-2}c_{4}^{2})+V_{\rm gas}(v_{\rm in}c_{2}^{\rm
in}-v_{\rm out}c_{2})\,$ $\displaystyle{\rm Pt}$
$\displaystyle\dot{N}_{3}=V_{\rm
surf}(-k_{1}c_{1}c_{3}+k_{-1}c_{5}-2k_{2}c_{2}c_{3}^{2}+2k_{-2}c_{4}^{2}+2k_{3}c_{4}c_{5})\,$
$\displaystyle{\rm PtO}$ $\displaystyle\dot{N}_{4}=V_{\rm
surf}(2k_{2}c_{2}c_{3}^{2}-2k_{-2}c_{4}^{2}-k_{3}c_{4}c_{5})\,$
$\displaystyle{\rm PtCO}\;$ $\displaystyle\dot{N}_{5}=V_{\rm
surf}(k_{1}c_{1}c_{3}-k_{-1}c_{5}-k_{3}c_{4}c_{5})\,$ $\displaystyle{\rm
CO_{2}}$ $\displaystyle\dot{N}_{6}=V_{\rm surf}k_{3}c_{4}c_{5}+V_{\rm
gas}(v_{\rm in}c_{6}^{\rm in}-v_{\rm out}c_{6})\,$

Here $v_{\rm in}$ and $v_{\rm out}$ are the flux rates (per unit volume).

For the slow variables this equation gives:

$\begin{split}&\dot{b}_{1}=2\dot{N}_{2}+\dot{N}_{4}=-V_{\rm
surf}k_{3}c_{4}c_{5}+2V_{\rm gas}(v_{\rm in}c_{2}^{\rm in}-v_{\rm
out}c_{2})\\\ &\dot{b}_{2}=\dot{N}_{1}+\dot{N}_{5}=-V_{\rm
surf}k_{3}c_{4}c_{5}+V_{\rm gas}(v_{\rm in}c_{1}^{\rm in}-v_{\rm out}c_{1})\\\
&\dot{b}_{3}=\dot{N}_{3}+\dot{N}_{4}+\dot{N}_{5}=0\\\ &\dot{N}_{6}=V_{\rm
surf}k_{3}c_{4}c_{5}+V_{\rm gas}(v_{\rm in}c_{6}^{\rm in}-v_{\rm
out}c_{6})\end{split}$ (40)

This system looks quite simple. Only one reaction,

PtO+PtCO$\to$CO2+2Pt (41)

is visible. If we know expressions for $c_{3,5}(b)$ then this reaction rate is
also known. In addition, to work with the rates of fluxes, the expressions for
$c_{1,2}(b)$ are needed.

The system of equations (37) is explicitly solvable but the result is quite
cumbersome. Therefore, let us consider its simplification without explicit
analytic solution. We assume the following smallness:

$b_{1}\gg N_{4}\,,\;\;b_{2}\gg N_{5}\,$ (42)

Together with this smallness assumptions equations (37) give:

$\begin{split}&c_{3}=\frac{b_{3}}{V_{\rm
surf}\left(1+\frac{k_{1}}{k_{-1}}\frac{b_{2}}{V_{\rm
gas}}+\sqrt{\frac{1}{2}\frac{k_{2}}{k_{-2}}\frac{b_{1}}{V_{\rm
gas}}}\right)}\\\
&c_{4}=\sqrt{\frac{1}{2}\frac{k_{2}}{k_{-2}}\frac{b_{1}}{V_{\rm
gas}}}\frac{b_{3}}{V_{\rm surf}\left(1+\frac{k_{1}}{k_{-1}}\frac{b_{2}}{V_{\rm
gas}}+\sqrt{\frac{1}{2}\frac{k_{2}}{k_{-2}}\frac{b_{1}}{V_{\rm
gas}}}\right)}\\\ &c_{5}=\frac{k_{1}}{k_{-1}}\frac{b_{2}}{V_{\rm
gas}}\frac{b_{3}}{V_{\rm surf}\left(1+\frac{k_{1}}{k_{-1}}\frac{b_{2}}{V_{\rm
gas}}+\sqrt{\frac{1}{2}\frac{k_{2}}{k_{-2}}\frac{b_{1}}{V_{\rm
gas}}}\right)}\end{split}$ (43)

In this approximation, we have for the reaction (41) rate

$r=k_{3}c_{4}c_{5}=k_{3}\frac{k_{1}}{k_{-1}}\sqrt{\frac{1}{2}\frac{k_{2}}{k_{-2}}}\frac{\sqrt{b_{1}}b_{2}}{V_{\rm
gas}^{3/2}}\frac{b_{3}^{2}}{V_{\rm
surf}^{2}\left(1+\frac{k_{1}}{k_{-1}}\frac{b_{2}}{V_{\rm
gas}}+\sqrt{\frac{1}{2}\frac{k_{2}}{k_{-2}}\frac{b_{1}}{V_{\rm
gas}}}\right)^{2}}$

This expression gives the closure for the slow QE equations (40).

####  4.2.4 4.2.3. Discussion of the QE procedure for Chemical Kinetics

We finalize here the illustration of the general QE procedure for chemical
kinetics. As we can see, the simple analytic description of the QE
approximation is available when the fast reactions have no joint reagents. In
general case, we need either a numerical solver for (24) or some additional
hypotheses about smallness. Michaelis and Menten used, in addition to the QE
approach, the hypothesis about smallness of the amount of intermediate
complexes. This is the typical QSS hypothesis. The QE approximation was
modified and further developed by many authors. In particular, a computational
optimization approach for the numerical approximation of slow attracting
manifolds based on entropy-related and geometric extremum principles for
reaction trajectories was developed Lebiedz2010 .

Of course, validity of all the simplification hypotheses is a crucial
question. For example, for the CO oxidation, if we accept the hypothesis about
the quasiequilibrium adsorption then we get a simple dynamics which
monotonically tends to the steady state. The state of the surface is
unambiguously presented as a continuous function of the gas composition. The
pure QSS hypothesis results for the Langmuir-Hinshelwood reaction mechanism
(38) without quasiequilibrium adsorption in bifurcations and the multiplicity
of steady states Yablonskii1991 . The problem of validity of simplifications
cannot be solved as a purely theoretical question without the knowledge of
kinetic constants or some additional experimental data.

##  5 General Kinetics with Fast Intermediates Present in Small Amount

###  5.1 Stoichiometry of Complexes

In this Section, we return to the very general reaction network.

Let us call all the formal sums that participate in the stoichiometric
equations (17), the complexes. The set of complexes for a given reaction
mechanism (17) is $\Theta_{1},\ldots,\Theta_{q}$. The number of complexes
$q\leq 2m$ (two complexes per elementary reaction, as the maximum) and it is
possible that $q<2m$ because some complexes may coincide for different
reactions.

A complex $\Theta_{i}$ is a formal sum
$\Theta_{i}=\sum_{j=1}^{n}\nu_{ij}A_{j}=(\nu_{i},A)$, where $\nu_{i}$ is a
vector with coordinates $\nu_{ij}$.

Each elementary reaction (17) may be represented in the form
$\Theta^{-}_{\rho}\to\Theta^{+}_{\rho}$, where $\Theta^{\pm}_{\rho}$ are the
complexes which correspond to the right and the left sides (17). The whole
mechanism is naturally represented as a digraph of transformation of
complexes: Vertices are complexes and edges are reactions. This graph gives a
convenient tool for the reaction representation and is often called the
“reaction graph”.

Let us consider a simple example: 18 elementary reactions (9 pairs of mutually
reverse reactions) from the hydrogen combustion mechanism (see, for example,
Conaireatal2004 ).

$\begin{array}[]{ll}{\rm H+O_{2}\rightleftharpoons O+OH;}&{\rm
O+H_{2}\rightleftharpoons H+OH;}\\\ {\rm OH+H_{2}\rightleftharpoons
H+H_{2}O;}&{\rm O+H_{2}O\rightleftharpoons 2OH;}\\\ {\rm
HO_{2}+H\rightleftharpoons H_{2}+O_{2};}&{\rm HO_{2}+H\rightleftharpoons
2OH;}\\\ {\rm H+OH+M\rightleftharpoons H_{2}O+M;}&{\rm
H+O_{2}+M\rightleftharpoons HO_{2}+M;}\\\ {\rm H_{2}O_{2}+H\rightleftharpoons
H_{2}+HO_{2}}&\end{array}$ (44)

There are 16 different complexes here:

$\displaystyle{\rm\Theta_{1}=H+O_{2},\,\Theta_{2}=O+OH,\,\Theta_{3}=O+H_{2},\,\Theta_{4}=H+OH,}$
$\displaystyle{\rm\Theta_{5}=OH+H_{2},\rm\Theta_{6}=H+H_{2}O,\,\Theta_{7}=O+H_{2}O,\,\Theta_{8}=2OH,\,}$
$\displaystyle{\rm\Theta_{9}=HO_{2}+H,\rm\Theta_{10}=H_{2}+O_{2},\,\Theta_{11}=H+OH+M,\,}$
$\displaystyle{\rm\Theta_{12}=H_{2}O+M,\,\rm\Theta_{13}=H+O_{2}+M,\,\Theta_{14}=HO_{2}+M,}$
$\displaystyle{\rm\Theta_{15}=H_{2}O_{2}+H,\,\Theta_{16}=H_{2}+HO_{2}}\,$

The reaction set (44) can be represented as

$\displaystyle\Theta_{1}\rightleftharpoons\Theta_{2},\,\Theta_{3}\rightleftharpoons\Theta_{4},\,\Theta_{5}\rightleftharpoons\Theta_{6},\,\Theta_{7}\rightleftharpoons\Theta_{8}\rightleftharpoons\Theta_{9}\leftrightharpoons\Theta_{10},$
$\displaystyle\Theta_{11}\rightleftharpoons\Theta_{12},\,\Theta_{13}\rightleftharpoons\Theta_{14},\,\Theta_{15}\rightleftharpoons\Theta_{16}\,$

We can see that this digraph of transformation of complexes has a very simple
structure: There are five isolated pairs of complexes and one connected group
of four complexes.

###  5.2 Stoichiometry of Compounds

For each complex $\Theta_{j}$ we introduce an additional component $B_{j}$, an
intermediate compound and $B^{\pm}_{\rho}$ are those compounds $B_{j}$ ($1\leq
j\leq q$), which correspond to the right and left sides of reaction (17).

We call these components “compounds” following the English translation of the
original Michaelis-Menten paper MichaelisMenten1913  and keep “complexes” for
the formal linear combinations $\Theta_{j}$.

An extended reaction mechanism includes two types of reactions: Equilibration
between a complex and its compound ($q$ reactions, one for each complex)

$\Theta_{j}\rightleftharpoons B_{j}$ (45)

and transformation of compounds $B_{\rho}^{-}\to B_{\rho}^{+}$ ($m$ reactions,
one for each elementary reaction from (17). So, instead of the reaction (17)
we can write

$\sum_{i}\alpha_{\rho i}A_{i}\rightleftharpoons B_{\rho}^{-}\to
B_{\rho}^{+}\rightleftharpoons\sum_{i}\beta_{\rho i}A_{i}\,$ (46)

Of course, if the input or output complexes coincide for two reactions then
the corresponding equilibration reactions also coincide.

It is useful to visualize the reaction scheme. In Figure 1 we represent the
$2n$-tail scheme of an elementary reaction sequence (46) which is an extension
of the elementary reaction (17).

The reactions between compounds may have several channels (Figure 2): One
complex may transform to several other complexes.

The reaction mechanism is a set of multichannel transformations (Figure 2) for
all input complexes. In Figure 2 we grouped together the reactions with the
same input complex. Another representation of the reaction mechanism is based
on the grouping of reactions with the same output complex. Below, in the
description of the complex balance condition, we use both representations.

The extended list of components includes $n+q$ components: $n$ initial species
$A_{i}$ and $q$ compounds $B_{j}$. The corresponding composition vector
$N^{\oplus}$ is a direct sum of two vectors, the composition vector for
initial species, $N$, with coordinates $N_{i}$ ($i=1,\ldots,n$) and the
composition vector for compounds, $\Upsilon$, with coordinates $\Upsilon_{j}$
($j=1,\ldots,q$): $N^{\oplus}=N\oplus\Upsilon$.

The space of composition vectors $E$ is a direct sum of $n$-dimensional
$E_{A}$ and $q$-dimensional $E_{B}$: $E=E_{A}\oplus E_{B}$.

For concentrations of $A_{i}$ we use the notation $c_{i}$ and for
concentrations of $B_{j}$ we use $\varsigma_{j}$.

The stoichiometric vectors for reactions $\Theta_{j}\rightleftharpoons B_{j}$
(45) are direct sums: $g^{j}=-\nu_{j}\oplus e_{j}$, where $e_{j}$ is the $j$th
standard basis vector of the space $R^{q}=E_{B}$, the coordinates of $e_{j}$
are $e_{jl}=\delta_{jl}$:

$g^{j}=(-\nu_{j1},-\nu_{j2},\ldots,-\nu_{jn},\underbrace{0,\ldots,0,1}_{l},0,\ldots,0)$
(47)

The stoichiometric vectors of equilibration reactions (45) are linearly
independent because there exists exactly one vector for each $l$.

The stoichiometric vectors $\gamma^{jl}$ of reactions $B_{j}\to B_{l}$ belong
entirely to $E_{B}$. They have $j$th coordinate $-1$, $l$th coordinate $+1$
and other coordinates are zeros.

To exclude some degenerated cases a hypothesis of weak reversibility is
accepted. Let us consider a digraph with vertices $\Theta_{i}$ and edges,
which correspond to reactions from (17). The system is weakly reversible if
for any two vertices $\Theta_{i}$ and $\Theta_{j}$, the existence of an
oriented path from $\Theta_{i}$ to $\Theta_{j}$ implies the existence of an
oriented path from $\Theta_{j}$ to $\Theta_{i}$.

Of course, this weak reversibility property is equivalent to weak
reversibility of the reaction network between compounds $B_{j}$.

Figure 2: A multichannel view on the complex transformation. The hidden
reactions between compounds are included in an oval $\mathbf{S}$.

###  5.3 Energy, Entropy and Equilibria of Compounds

In this section, we define the free energy of the system. The basic hypothesis
is that the compounds are the small admixtures to the system, that is, the
amount of compounds $B_{j}$ is much smaller than amount of initial components
$A_{i}$. Following this hypothesis, we neglect the energy of interaction
between compounds, which is quadratic in their concentrations because in the
low density limit we can neglect the correlations between particles if the
potential of their interactions decay sufficiently fast when the distance
between particles goes to $\infty$ Balescu . We take the energy of their
interaction with $A_{i}$ in the linear approximation, and use the perfect
entropy for $B_{i}$. These standard assumptions for a small admixtures give
for the free energy:

$F=Vf(c,T)+VRT\sum_{j=1}^{q}\varsigma_{j}\left(\frac{u_{j}(c,T)}{RT}+\ln\varsigma_{j}-1\right)$
(48)

Let us introduce the standard equilibrium concentrations for $B_{j}$. Due to
the Boltzmann distribution ($\exp(-u/RT)$) and formula (48)

$\varsigma_{j}^{*}(c,T)=\frac{1}{Z}\exp\left(-\frac{u_{j}(c,T)}{RT}\right)$
(49)

where $1/Z$ is the normalization factor. Let us select here the normalization
$Z=1$ and write:

$F=Vf(c,T)+VRT\sum_{j=1}^{q}\varsigma_{j}\left(\ln\left(\frac{\varsigma_{j}}{\varsigma_{j}^{*}(c,T)}\right)-1\right)\
$ (50)

We assume that the standard equilibrium concentrations
$\varsigma_{j}^{*}(c,T)$ are much smaller than the concentrations of $A_{i}$.
It is always possible because functions $u_{j}$ are defined up to an additive
constant.

The formula for free energy is necessary to define the fast equilibria (45).
Such an equilibrium is the minimizer of the free energy on the straight line
parameterized by $a$: $c_{i}=c_{i}^{0}-a\nu_{ji}$, $\varsigma_{j}=a$.

If we neglect the products
$\varsigma_{j}\partial\varsigma_{j}^{*}(c,T)/\partial c_{i}$ as the second
order small quantities then the minimizers have the very simple form:

$\vartheta_{j}=\sum_{i}\nu_{ji}\frac{\mu_{i}(c,T)}{RT}$ (51)

or

$\varsigma_{j}=\varsigma^{*}_{j}(c,T)\exp\left(\frac{\sum_{i}\nu_{ji}\mu_{i}(c,T)}{RT}\right)$
(52)

where

$\mu_{i}=\frac{\partial f(c,T)}{\partial c_{i}}$

is the chemical potential of $A_{i}$ and

$\vartheta_{j}=\ln\left(\frac{\varsigma_{j}}{\varsigma^{*}_{j}}\right)$

($RT\vartheta_{j}=\frac{1}{V}\frac{\partial F}{\partial\varsigma_{j}}$ is the
chemical potential of $B_{j}$).

The thermodynamic equilibrium of the system of reactions $B_{j}\to B_{l}$ that
corresponds to the reactions (46) is the free energy minimizer under given
values of the conservation laws.

For the systems with fixed volume, the stoichiometric conservation laws of the
monomolecular system of reaction are sums of the concentrations of $B_{j}$
which belong to the connected components of the reaction graph. Under the
hypothesis of weak reversibility there is no other conservation law. Let the
graph of reactions $B_{j}\to B_{l}$ have $d$ connected components $C_{s}$ and
let $V_{s}$ be the set of indexes of those $B_{j}$ which belong to $C_{s}$:
$B_{j}\in C_{s}$ if and only if $j\in V_{s}$. For each $C_{s}$ there exists a
stoichiometric conservation law

$\beta_{s}=\sum_{j\in V_{s}}\varsigma_{j}=const$ (53)

For any set of positive values of $\beta_{s}$ ($s=1,\ldots,d$) and given $c,T$
there exists a unique conditional maximizer $\varsigma^{\rm eq}_{j}$ of the
free energy (50): For the compound $B_{j}$ from the $s$th connected component
($j\in V_{s}$) this equilibrium concentration is

$\varsigma^{\rm eq}_{j}=\beta_{s}\frac{\varsigma_{j}^{*}(c,T)}{\sum_{l\in
V_{s}}\varsigma_{j}^{*}(c,T)}$ (54)

The positive values of concentrations $\varsigma_{j}$ are the equilibrium
concentrations (54) for some values of $\beta_{s}$ if and only if for any
$s=1,\ldots,d$ and all $j,l\in V_{s}$

$\vartheta_{j}=\vartheta_{l}$ (55)

($\vartheta_{j}=\ln(\varsigma_{j}/\varsigma_{j}^{*})$). This means that
compounds are in equilibrium in every connected component $C_{s}$ the chemical
potentials of compounds coincide in each component $C_{s}$. The system of
equations (55) together with the equilibrium conditions (52) constitute the
equilibrium of the systems. All the equilibria form a linear subspace in the
space with coordinates $\mu_{i}/RT$ ($i=1,\ldots,n$) and $\vartheta_{j}$
($j=1,\ldots,q$).

In the expression for the free energy (50) we do not assume anything special
about free energy of the mixture of $A_{i}$. The density of this free energy,
$f(c,T)$, may be an arbitrary smooth function (later, we will add the standard
assumption about convexity of $f(c,T)$ as a function of $c$). For the
compounds $B_{i}$, we assume that they form a very small addition to the
mixture of $A_{i}$, neglect all quadratic terms in concentrations of $B_{i}$
and use the entropy of the perfect systems, $p\ln p$, for this small
admixture.

This approach results in the explicit expressions for the fast equilibria (52)
and expression of the equilibrium compound concentrations through the values
of the stoichiometric conservation laws (54).

###  5.4 Markov Kinetics of Compounds

For the kinetics of compounds transformations $B_{j}\to B_{l}$, the same
hypothesis of the smallness of concentrations leads to the only reasonable
assumption: The linear (monomolecular) kinetics with the rate constant
$\kappa_{lj}>0$. This “constant” is a function of $c,T$: $\kappa_{lj}(c,T)$.
The order of indexes at $\kappa$ is inverse to the order of them in reaction:
$\kappa_{lj}=\kappa_{l\leftarrow j}$.

The master equation for the concentration of $B_{j}$ gives:

$\frac{{\mathrm{d}}\varsigma_{j}}{{\mathrm{d}}t}=\sum_{l,\,l\neq
j}\left(\kappa_{jl}\varsigma_{l}-\kappa_{lj}\varsigma_{j}\right)$ (56)

It is necessary to find when this kinetics respect thermodynamics, i.e., when
the free energy decreases due to the system (56). The necessary and sufficient
condition for matching the kinetics and thermodynamics is: The standard
equilibrium $\varsigma^{*}$ (49) should be an equilibrium for (56), that is,
for every $j=1,\ldots,q$

$\sum_{l,\,l\neq j}\kappa_{jl}\varsigma_{l}^{*}=\sum_{l,\,l\neq
j}\kappa_{lj}\varsigma_{j}^{*}$ (57)

This condition is necessary because the standard equilibrium is the free
energy minimizer for given $c,T$ and
$\sum_{j}\varsigma_{j}=\sum_{j}\varsigma_{j}^{*}$. The sum
$\sum_{j}\varsigma_{j}$ conserves due to (56). Therefore, if we assume that
$F$ decreases monotonically due to (56) then the point of conditional minimum
of $F$ on the plane $\sum_{j}\varsigma_{j}=const$ (under given $c,T$) should
be an equilibrium point for this kinetic system. This condition is sufficient
due to the Morimoto $H$-theorem (see Appendix 2).

For a weakly reversible system, the set of the conditional minimizers of the
free energy (54) coincides with with the set of positive equilibria for the
master equations (56) (see Equation (132) in Appendix 2).

###  5.5 Thermodynamics and Kinetics of the Extended System

In this section, we consider the complete extended system, which consists of
species $A_{i}$ ($i=1,\ldots,n$) and compounds $B_{j}$ ($j=1,\ldots,q$) and
includes reaction of equilibration (45) and transformations of compounds
$B_{j}\to B_{l}$ which correspond to the reactions (46).

Thermodynamic properties of the system are summarized in the free energy
function (50). For kinetics of compounds we accept the Markov model (56) with
the equilibrium condition (57), which guarantees matching between
thermodynamics and kinetics.

For the equilibration reactions (45) we select a very general form of the
kinetic law. The only requirement is: This reaction should go to its
equilibrium, which is described as the conditional minimizer of free energy
$F$ (52). For each reaction $\Theta_{j}\rightleftharpoons B_{j}$ (where the
complex is a formal combination: $\Theta_{j}=\sum_{i}\nu_{ji}A_{i}$) we
introduce the reaction rate $w_{j}$. This rate should be positive if

$\vartheta_{j}<\sum_{i}\nu_{ji}\frac{\mu_{i}(c,T)}{RT}$ (58)

and negative if

$\vartheta_{j}>\sum_{i}\nu_{ji}\frac{\mu_{i}(c,T)}{RT}$ (59)

The general way to satisfy these requirement is to select $q$ continuous
function of real variable $w_{j}(x)$, which are negative if $x>0$ and positive
if $x<0$. For the equilibration rates we take

$w_{j}=w_{j}\left(\vartheta_{j}-\sum_{i}\nu_{ji}\frac{\mu_{i}(c,T)}{RT}\right)$
(60)

If several dynamical systems defined by equations $\dot{x}=J_{1}$, …
$\dot{x}=J_{v}$ on the same space have the same Lyapunov function $F$, then
for any conic combination $J=\sum_{k}a_{k}J_{k}$ ($a_{k}\geq 0$,
$\sum_{k}a_{k}>0$) the dynamical system $\dot{x}=J$ also has the Lyapunov
function $F$.

The free energy (50) decreases monotonically due to any reaction
$\Theta_{j}\rightleftharpoons B_{j}$ with reaction rate $w_{j}$ (60) and also
due to the Markov kinetics (56) with the equilibrium condition (57).
Therefore, the free energy decreases monotonically due to the following
kinetic system:

$\begin{split}&\frac{{\mathrm{d}}c_{i}}{{\mathrm{d}}t}=-\sum_{j=1}^{q}\nu_{ji}w_{j}\,\\\
&\frac{{\mathrm{d}}\varsigma_{j}}{{\mathrm{d}}t}=w_{j}+\sum_{l,\,l\neq
j}\left(\kappa_{jl}\varsigma_{l}-\kappa_{lj}\varsigma_{j}\right)\end{split}$
(61)

where the coefficients $\kappa_{jl}$ satisfy the matching condition (57).

This general system (61) describes kinetics of extended system and satisfies
all the basic conditions (thermodynamics and smallness of compound
concentrations). In the next sections we will study the QE approximations to
this system and exclude the unknown functions $w_{j}$ from it.

###  5.6 QE Elimination of Compounds and the Complex Balance Condition

In this section, we use the QE formalism developed for chemical kinetics in
Section 4 for simplification of the compound kinetics.

First of all, let us describe $L^{\perp}$, where the space $L$ is the subspace
in the extended concentration space spanned by the stoichiometric vectors of
fast equilibration reactions (45). The stoichiometric vector for the
equilibration reactions have a very special structure (47). Dimension of the
space $L$ is equal to the number of complexes: $\dim L=q$. Therefore,
dimension of $L^{\perp}$ is equal to the number of components $A_{i}$: $\dim
L^{\perp}=n$. For each $A_{i}$ we will find a vector $b_{i}\in L^{\perp}$ that
has the following first $n$ coordinates: $b_{ik}=\delta_{ik}$ for
$k=1,\ldots,n$. The condition $(b_{i},g_{j})=0$ gives immediately:
$b_{i,n+j}=\nu_{ji}$. Finally,

$b_{i}=(\overbrace{\underbrace{0,\ldots,0,1}_{i},0,\ldots,0}^{n},\nu_{1i},\nu_{2i},\ldots,\nu_{qi})$
(62)

The corresponding slow variables are

$b_{i}(c,\varsigma)=c_{i}+\sum_{j}\varsigma_{j}\nu_{ji}$ (63)

In the QE approximation all $w_{j}=0$ and the kinetic equations (61) give in
this approximation

$\frac{{\mathrm{d}}b_{i}}{{\mathrm{d}}t}=\sum_{lj,\,l\neq
j}(\kappa_{jl}\varsigma_{l}-\kappa_{lj}\varsigma_{j})\nu_{ji}$ (64)

In these equations, we have to use the dependence $\varsigma(b)$. Here we use
the QSS Michaelis and Menten assumption: The compounds are present in small
amounts

$c_{i}\gg\varsigma_{j}$

In this case, we can take $b_{i}$ instead of $c_{i}$ (i.e., take $\mu(b,T)$
instead of $\mu(c,T)$) in the formulas for equilibria (52):

$\varsigma_{j}=\varsigma^{*}_{j}(b,T)\exp\left(\frac{\sum_{i}\nu_{ji}\mu_{i}(b,T)}{RT}\right)$
(65)

In the final form of the QE kinetic equation there remain two “offprints” of
the compound kinetics: Two sets of functions $\varsigma^{*}_{j}(b,T)\geq 0$
and $\kappa_{jl}(b,T)\geq 0$. These functions are connected by the identity
(57). The final form of the equations is

$\begin{split}\frac{{\mathrm{d}}b_{i}}{{\mathrm{d}}t}=&\sum_{lj,\,l\neq
j}\left(\kappa_{jl}\varsigma^{*}_{l}(b,T)\exp\left(\frac{\sum_{i}\nu_{li}\mu_{i}(b,T)}{RT}\right)-\kappa_{lj}\varsigma^{*}_{j}(b,T)\exp\left(\frac{\sum_{i}\nu_{ji}\mu_{i}(b,T)}{RT}\right)\right)\nu_{ji}\end{split}$
(66)

The identity (57), $\sum_{l,\,l\neq
j}\kappa_{jl}\varsigma_{l}^{*}=\sum_{l,\,l\neq
j}\kappa_{lj}\varsigma_{j}^{*}$, provides a sufficient condition for
decreasing of free energy due to the kinetic equations (66). This is a direct
consequence of two theorem: The theorem about the preservation of entropy
production in the QE approximations (see Section 2 and Appendix 1) and the
Morimoto $H$-theorem (see Appendix 2). Indeed, in the QE state the
equilibrated reactions (45) $\Theta_{j}\rightleftharpoons B_{j}$ do not
produce entropy and all changes in the total free energy are caused by the
Markov kinetics $B_{i}\to B_{j}$. Due to the Morimoto $H$-theorem this change
is negative: The Markov kinetics decrease the perfect free energy of compounds
and do not affect the free energy of $A_{i}$. In the QE approximation, the
concentrations of $A_{i}$ are changing together with concentrations of $B_{j}$
because of the equilibrium conditions for reactions
$\Theta_{j}\rightleftharpoons B_{j}$. Due to the theorem of preservation of
the entropy production, the time derivative of the total free energy in this
QE dynamics coincides with the time derivative of the free energy of $B_{j}$
due to Markov kinetics. In addition to this proof, in Section 6 below we give
the explicit formula for entropy production in (66) and direct proof of its
positivity.

Let us stress that the functions $\varsigma^{*}_{j}(b,T)$ and
$\kappa_{jl}(b,T)$ participate in equations (66) and in identity (57) in the
form of the product. Below we use for this product a special notation:

$\varphi_{jl}(b,T)=\kappa_{jl}(b,T)\varsigma^{*}_{l}(b,T)\;(j\neq l)$ (67)

We call this function $\varphi_{jl}(b,T)$ the kinetic factor. The identity
(57) for the kinetic factor is

$\sum_{l,\,l\neq j}\varphi_{jl}(b,T)=\sum_{l,\,l\neq j}\varphi_{lj}(b,T)\mbox{
for all }j$ (68)

We call the thermodynamic factor (or the Boltzmann factor) the second
multiplier in the reaction rates

$\Omega_{l}(b,T)=\exp\left(\frac{\sum_{i}\nu_{li}\mu_{i}(b,T)}{RT}\right)$
(69)

In this notation, the kinetic equations (66) have a simple form

$\frac{{\mathrm{d}}b_{i}}{{\mathrm{d}}t}=\sum_{lj,\,l\neq
j}(\varphi_{jl}(b,T)\Omega_{l}(b,T)-\varphi_{lj}(b,T)\Omega_{j}(b,T))\nu_{ji}$
(70)

The general equations (70) have the form of “sum over complexes”. Let us
return to the more usual “sum over reactions” form. An elementary reaction
corresponds to the pair of complexes $\Theta_{l},\Theta_{j}$ (46). It has the
form $\Theta_{l}\to\Theta_{j}$ and the reaction rate is
$r=\varphi_{jl}\Omega_{l}$. In the right hand side in (70) this reaction
appears twice: first time with sign “$+$” and the vector coefficient $\nu_{j}$
and the second time with sign “$-$” and the vector coefficient $\nu_{l}$. The
stoichiometric vector of this reaction is $\gamma=\nu_{j}-\nu_{l}$. Let us
enumerate the elementary reactions by index $\rho$, which corresponds to the
pair $(j,l)$. Finally, we transform (46) into the sum over reactions form

$\begin{split}\frac{{\mathrm{d}}b_{i}}{{\mathrm{d}}t}&=\sum_{l,j,\,l\neq
j}\varphi_{jl}(b,T)\Omega_{l}(b,T)(\nu_{ji}-\nu_{li})\\\
&=\sum_{\rho}\varphi_{\rho}(b,T)\Omega_{\rho}(b,T)\gamma_{{\rho}i}\end{split}$
(71)

In the vector form it looks as follows:

$\frac{{\mathrm{d}}b}{{\mathrm{d}}t}=\sum_{\rho}\varphi_{\rho}(b,T)\Omega_{\rho}(b,T)\gamma_{{\rho}}$
(72)

###  5.7 The Big Michaelis-Menten-Stueckelberg Theorem

Let us summarize the results of our analysis in one statement.

Let us consider the reaction mechanism illustrated by Figure 2 (46):

$\sum_{i}\alpha_{\rho i}A_{i}\rightleftharpoons B_{\rho}^{-}\to
B_{\rho}^{+}\rightleftharpoons\sum_{i}\beta_{\rho i}A_{i}$

under the following asymptotic assumptions:

  1. 1.

Concentrations of the compounds $B_{\rho}$ are close to their quasiequilibrium
values (65)

$\varsigma_{j}=(1+\delta)\varsigma_{j}^{\rm
QE}=(1+\delta)\varsigma^{*}_{j}(b,T)\exp\left(\frac{\sum_{i}\nu_{ji}\mu_{i}(b,T)}{RT}\right)\,,\;\;\delta\ll
1$

(this may be due to the fast reversible reactions in (46));

  2. 2.

Concentrations of the compounds $B_{\rho}$ are much smaller than the
concentrations of the components $A_{i}$: There is a small positive parameter
$\varepsilon\ll 1$, $\varsigma^{*}_{j}=\varepsilon\xi^{*}_{j}$ and
$\xi^{*}_{j}$ do not depend on $\varepsilon$;

  3. 3.

Kinetics of transitions between compounds $B_{i}\to B_{j}$ is linear (Markov)
kinetics with reaction rate constants
$k_{ji}=\frac{1}{\varepsilon}\kappa_{ji}$.

Under these assumptions, in the asymptotic $\delta,\varepsilon\to 0$,
$\delta,\varepsilon>0$ kinetics of components $A_{i}$ may be described by the
reaction mechanism

$\sum_{i}\alpha_{\rho i}A_{i}\to\sum_{i}\beta_{\rho i}A_{i}$

with the reaction rates

$r_{\rho}=\varphi_{\rho}\exp\left(\frac{(\alpha_{\rho},{\mu})}{RT}\right)$

where the kinetic factors $\varphi_{\rho}$ satisfy the condition (68):

$\sum_{\rho,\,\alpha_{\rho}=\mathbf{v}}\varphi_{\rho}\equiv\sum_{\rho,\,\beta_{\rho}=\mathbf{v}}\varphi_{\rho}$

for any vector $\mathbf{v}$ from the set of all vectors
$\\{\alpha_{\rho},\beta_{\rho}\\}$. This statement includes the generalized
mass action law for $r_{\rho}$ and the balance identity for kinetic factors
that is sufficient for the entropy growth as it is shown in the next Section
6.

##  6 General Kinetics and Thermodynamics

###  6.1 General Formalism

To produce the general kinetic law and the complex balance conditions, we use
“construction staging”: The intermediate complexes with fast equilibria, the
Markov kinetics and other important and interesting physical and chemical
hypothesis.

In this section, we delete these construction staging and start from the forms
(69), (72) as the basic laws. We use also the complex balance conditions (68)
as a hint for the general conditions which guarantee accordance between
kinetics and thermodynamics.

Let us consider a domain $U$ in $n$-dimensional real vector space $E$ with
coordinates $N_{1},\ldots,N_{n}$. For each $N_{i}$ a symbol (component)
$A_{i}$ is given. A dimensionless entropy (or free entropy, for example,
Massieu, Planck, or Massieu-Planck potential which correspond to the selected
conditions Callen1985 ) $S(N)$ is defined in $U$. “Dimensionless” means that
we use $S/R$ instead of physical $S$. This choice of units corresponds to the
informational entropy ($p\ln p$ instead of $k_{\rm B}p\ln p$).

The dual variables, potentials, are defined as the partial derivatives of $S$:

$\check{\mu}_{i}=-\frac{\partial S}{\partial N_{i}}$ (73)

Warning: This definition differs from the chemical potentials (22) by the
factor ${1}/{RT}$: For constant volume the Massieu-Planck potential is $-F/T$
and we, in addition, divide it on $R$. On the other hand, we keep the same
sign as for the chemical potentials, and this differs from the standard
Legendre transform for $S$. (It is the Legendre transform for function $-S$).

The reaction mechanism is defined by the stoichiometric equations (17)

$\sum_{i}\alpha_{\rho i}A_{i}\to\sum_{i}\beta_{\rho i}A_{i}$

($\rho=1,\ldots,m$). In general, there is no need to assume that the
stoichiometric coefficients $\alpha_{\rho i},\beta_{\rho i}$ are integers.

The assumption that they are nonnegative, $\alpha_{\rho i}\geq 0,\beta_{\rho
i}\geq 0$, may be needed to prove that the kinetic equations preserve
positivity of $N_{i}$. If $N_{i}$ is the number of particles then it is a
natural assumption but we can use other extensive variables instead, for
example, we included energy in the list of variables to describe the non-
isothermal processes BykGOrYab1982 . In this case, the coefficient
$\alpha_{U}$ for the energy component $A_{U}$ in an exothermic reaction is
negative.

So, for variables that are positive (bounded from below) by their physical
sense, we will use the inequalities $\alpha_{\rho i}\geq 0,\beta_{\rho i}\geq
0$, when necessary, but in general, for arbitrary extensive variables, we do
not assume positivity of stoichiometric coefficients. As it is usually, the
stoichiometric vector of reaction is
$\gamma_{\rho}=\beta_{\rho}-\alpha_{\rho}$ (the “gain minus loss” vector).

For each reaction, a nonnegative quantity, reaction rate $r_{\rho}$ is
defined. We assume that this quantity has the following structure:

$r_{\rho}=\varphi_{\rho}\exp(\alpha_{\rho},\check{\mu})$ (74)

where $(\alpha_{\rho},\check{\mu})=\sum_{i}\alpha_{\rho i}\check{\mu}_{i}$.

In the standard formalism of chemical kinetics the reaction rates are
intensive variables and in kinetic equations for $N$ an additional factor—the
volume—appears. For heterogeneous systems, there may be several “volumes”
(including interphase surfaces).

Each reaction has it own “volume”, an extensive variable $V_{\rho}$ (some of
them usually coincide), and we can write

$\frac{{\mathrm{d}}N}{{\mathrm{d}}t}=\sum_{\rho}V_{\rho}\gamma_{\rho}\varphi_{\rho}\exp(\alpha_{\rho},\check{\mu})$
(75)

In these notations, both the kinetic and the Boltzmann factors are intensive
(and local) characteristics of the system.

Let us, for simplicity of notations, consider a system with one volume, $V$
and write

$\frac{{\mathrm{d}}N}{{\mathrm{d}}t}=V\sum_{\rho}\gamma_{\rho}\varphi_{\rho}\exp(\alpha_{\rho},\check{\mu})$
(76)

Below we use the form (76). All our results will hold also for the multi-
volume systems (75) under one important assumption: The elementary reaction

$\sum_{i}\alpha_{\rho i}A_{i}\to\sum_{i}\beta_{\rho i}A_{i}\,$

goes in the same volume as the reverse reaction

$\sum_{i}\beta_{\rho i}A_{i}\to\sum_{i}\alpha_{\rho i}A_{i}\,$

or symbolically

$V_{\rho}^{+}=V_{\rho}^{-}$ (77)

If this condition (77) holds then the detailed balance conditions and the
complex balance conditions will hold separately in all volumes $V_{\rho}$.

An important particular case of (76) gives us the Mass Action Law. Let us take
the perfect free entropy

$S=-\sum_{i}N_{i}\left(\ln\left(\frac{c_{i}}{c_{i}^{*}}\right)-1\right)$ (78)

where $c_{i}=N_{i}/V\geq 0$ are concentrations and $c_{i}^{*}>0$ are the
standard equilibrium concentrations. Under isochoric conditions, $V=const$,
there is no difference between the choice of the main variables, $N$ or $c$.

For the perfect function (78)

$\check{\mu}_{i}=\ln\left(\frac{c_{i}}{c_{i}^{*}}\right)\,,\;\exp(\alpha_{\rho},\check{\mu})=\prod_{i}\left(\frac{c_{i}}{c_{i}^{*}}\right)^{\alpha_{\rho
i}}$ (79)

and for the reaction rate function (74) we get

$r_{\rho}=\varphi_{\rho}\prod_{i}\left(\frac{c_{i}}{c_{i}^{*}}\right)^{\alpha_{\rho
i}}$ (80)

The standard assumption for the Mass Action Law in physics and chemistry is
that $\varphi$ and $c^{*}$ are functions of temperature:
$\varphi_{\rho}=\varphi_{\rho}(T)$ and $c^{*}_{i}=c^{*}_{i}(T)$. To return to
the kinetic constants notation (20) we should write:

$\frac{\varphi_{\rho}}{\prod_{i}{c_{i}^{*}}^{\alpha_{\rho i}}}=k_{\rho}$

Equation (76) is the general form of the kinetic equation we would like to
study. In many senses, this form is too general before we impose restrictions
on the values of the kinetic factors. For physical and chemical systems,
thermodynamics is a source of restrictions:

  1. 1.

The energy of the Universe is constant.

  2. 2.

The entropy of the Universe tends to a maximum.

(R. Clausius, 1865 Clausius .)

The first sentence should be extended: The kinetic equations should respect
several conservation laws: Energy, amount of atoms of each kind (if there is
no nuclear reactions in the system) conservation of total probability and,
sometimes, some other conservation laws. All of them have the form of
conservation of values of some linear functionals: $l(N)=const$. If the input
and output flows are added to the system then

$\frac{{\mathrm{d}}l(N)}{{\mathrm{d}}t}=Vv^{\rm in}l^{\rm in}-v^{\rm out}l(N)$

where $v^{\rm in,out}$ are the input and output fluxes per unit volume,
$l^{\rm in}$ are the input densities (concentration). The standard requirement
is that every reaction respects all these conservation laws. The formal
expression of this requirement is:

$l(\gamma_{\rho})=0\mbox{ for all }\rho$ (81)

There is a special term for this conservation laws: The stoichiometric
conservation laws. All the main conservation laws are assumed to be the
stoichiometric ones.

Analysis of the stoichiometric conservation laws is a simple linear algebra
task: We have to find the linear functionals that annulate all the
stoichiometric vectors $\gamma_{\rho}$. In contrast, entropy is not a linear
function of $N$ and analysis of entropy production is not so simple.

In the next subsection we discuss various conditions which guarantee the
positivity of entropy production in kinetic equations (76).

###  6.2 Accordance Between Kinetics and Thermodynamics

####  6.2.1 6.2.1. General Entropy Production Formula

Let us calculate ${\mathrm{d}}S/{\mathrm{d}}t$ due to equations (76):

$\begin{split}\frac{{\mathrm{d}}S}{{\mathrm{d}}t}&=\sum_{i}\frac{\partial
S}{\partial N_{i}}\frac{{\mathrm{d}}N_{i}}{{\mathrm{d}}t}\\\
&=-\sum_{i}\check{\mu}_{i}V\sum_{\rho}\gamma_{\rho
i}\varphi_{\rho}\exp(\alpha_{\rho},\check{\mu})\\\
&=-V\sum_{\rho}(\gamma_{\rho},\check{\mu})\varphi_{\rho}\exp(\alpha_{\rho},\check{\mu})\end{split}$
(82)

An auxiliary function $\theta(\lambda)$ of one variable $\lambda\in[0,1]$ is
convenient for analysis of ${\mathrm{d}}S/{\mathrm{d}}t$ (it was studied by
Rozonoer and Orlov OrlovRozonoer1984 , see also G1 :

$\theta(\lambda)=\sum_{\rho}\varphi_{\rho}\exp[(\check{\mu},(\lambda\alpha_{\rho}+(1-\lambda)\beta_{\rho}))]$
(83)

With this function, the entropy production (82) has a very simple form:

$\frac{{\mathrm{d}}S}{{\mathrm{d}}t}=V\left.\frac{{\mathrm{d}}\theta(\lambda)}{{\mathrm{d}}\lambda}\right|_{\lambda=1}$
(84)

The auxiliary function $\theta(\lambda)$ allows the following interpretation.
Let us introduce the deformed stoichiometric mechanism with the stoichiometric
vectors,

$\alpha_{\rho}(\lambda)=\lambda\alpha_{\rho}+(1-\lambda)\beta_{\rho}\,,\;\beta_{\rho}(\lambda)=\lambda\beta_{\rho}+(1-\lambda)\alpha_{\rho}$

, which is the initial mechanism when $\lambda=1$, the inverted mechanism with
interchange of $\alpha$ and $\beta$ when $\lambda=0$, the trivial mechanism
(the left and right hand sides of the stoichiometric equations coincide) when
$\lambda=1/2$.

For the deformed mechanism, let us take the same kinetic factors and calculate
the Boltzmann factors with $\alpha_{\rho}(\lambda)$:

$r_{\rho}(\lambda)=\varphi_{\rho}\exp(\alpha_{\rho}(\lambda),\check{\mu})$

In this notation, the auxiliary function $\theta(\lambda)$ is a sum of
reaction rates for the deformed reaction mechanism:

$\theta(\lambda)=\sum_{\rho}r_{\rho}(\lambda)$

In particular, $\theta(1)=\sum_{\rho}r_{\rho}$, this is just the sum of
reaction rates.

Function $\theta(\lambda)$ is convex. Indeed

$\frac{{\mathrm{d}}^{2}\theta(\lambda)}{{\mathrm{d}}\lambda^{2}}=\sum_{\rho}\varphi_{\rho}(\gamma_{\rho},\check{\mu})^{2}\exp[(\check{\mu},(\lambda\alpha_{\rho}+(1-\lambda)\beta_{\rho}))]\geq
0$

This convexity gives the following necessary and sufficient condition for
positivity of entropy production:

$\frac{{\mathrm{d}}S}{{\mathrm{d}}t}>0\mbox{ if and only if
}\theta(\lambda)<\theta(1)\mbox{ for some }\lambda<1$

In several next subsections we study various important particular sufficient
conditions for positivity of entropy production.

####  6.2.2 6.2.2. Detailed Balance

The most celebrated condition which gives the positivity of entropy production
is the principle of detailed balance. Boltzmann used this principle to prove
his famous $H$-theorem Boltzmann .

Let us join elementary reactions in pairs:

$\sum_{i}\alpha_{\rho i}A_{i}\rightleftharpoons\sum_{i}\beta_{\rho i}A_{i}$
(85)

After this joining, the total amount of stoichiometric equations decreases. If
there is no reverse reaction then we can add it formally, with zero kinetic
factor. We will distinguish the reaction rates and kinetic factors for direct
and inverse reactions by the upper plus or minus:

$r_{\rho}^{+}=\varphi_{\rho}^{+}\exp(\alpha_{\rho},\check{\mu})\,,\;r_{\rho}^{-}=\varphi_{\rho}^{-}\exp(\beta_{\rho},\check{\mu})\,,\;r_{\rho}=r_{\rho}^{+}-r_{\rho}^{-}$

$\frac{{\mathrm{d}}N}{{\mathrm{d}}t}=V\sum_{\rho}\gamma_{\rho}r_{\rho}$ (86)

In this notation, the principle of detailed balance is very simple: The
thermodynamic equilibrium in the direction $\gamma_{\rho}$, given by the
standard condition $(\gamma_{\rho},\check{\mu})=0$, is equilibrium for the
corresponding pair of mutually reverse reactions from (85). For kinetic
factors this transforms into the simple and beautiful condition:

$\varphi_{\rho}^{+}\exp(\alpha_{\rho},\check{\mu})=\varphi_{\rho}^{-}\exp(\beta_{\rho},\check{\mu})\Leftrightarrow(\gamma_{\rho},\check{\mu})=0$

therefore

$\varphi_{\rho}^{+}=\varphi_{\rho}^{-}$ (87)

For the systems with detailed balance we can take
$\varphi_{\rho}=\varphi_{\rho}^{+}=\varphi_{\rho}^{-}$ and write for the
reaction rate:

$r_{\rho}=\varphi_{\rho}(\exp(\alpha_{\rho},\check{\mu})-\exp(\beta_{\rho},\check{\mu}))$

M. Feinberg called this kinetic law the “Marselin-De Donder” kinetics
Feinberg1972_a . This representation of the reaction rates gives for the
auxiliary function $\theta(\lambda)$:

$\theta(\lambda)=\sum_{\rho}\varphi_{\rho}(\exp[(\check{\mu},(\lambda\alpha_{\rho}+(1-\lambda)\beta_{\rho}))]+\exp[(\check{\mu},(\lambda\beta_{\rho}+(1-\lambda)\alpha_{\rho}))])$
(88)

Each term in this sum is symmetric with respect to change
$\lambda\mapsto(1-\lambda)$. Therefore, $\theta(1)=\theta(0)$ and, because of
convexity of $\theta(\lambda)$, $\theta^{\prime}(1)\geq 0$. This means
positivity of entropy production.

The principle of detailed balance is a sufficient but not a necessary
condition of the positivity of entropy production. This was clearly explained,
for example, by L. Onsager Onsager1931a ; Onsager1931b . Interrelations
between positivity of entropy production, Onsager reciprocal relations and
detailed balance were analyzed in detail by N.G. van Kampen VKampen1973 .

####  6.2.3 6.2.3. Complex Balance

The principle of detailed balance gives us $\theta(1)=\theta(0)$ and this
equality holds for each pair of mutually reverse reactions.

Let us start now from the equality $\theta(1)=\theta(0)$. We return to the
initial stoichiometric equations (17) without joining the direct and reverse
reactions. The equality reads

$\sum_{\rho}\varphi_{\rho}\exp(\check{\mu},\alpha_{\rho})=\sum_{\rho}\varphi_{\rho}\exp(\check{\mu},\beta_{\rho})$
(89)

Exponential functions $\exp(\check{\mu},y)$ form linearly independent family
in the space of functions of $\check{\mu}$ for any finite set of pairwise
different vectors $y$. Therefore, the following approach is natural: Let us
equalize in (89) the terms with the same Boltzmann-type factor
$\exp(\check{\mu},y)$. Here we have to return to the complex-based
representation of reactions (see Section 5.1).

Let us consider the family of vectors $\\{\alpha_{\rho},\beta_{\rho}\\}$
($\rho=1,\ldots,m$). Usually, some of these 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

$R_{j}^{+}=\\{\rho\,|\,\alpha_{\rho}=y_{j}\\}\,,\;R_{j}^{-}=\\{\rho\,|\,\beta_{\rho}=y_{j}\\}$

We can rewrite the equality (89) in the form

$\sum_{j=1}^{q}\exp(\check{\mu},y_{j})\left[\sum_{\rho\in
R_{j}^{+}}\varphi_{\rho}-\sum_{\rho\in R_{j}^{-}}\varphi_{\rho}\right]=0$ (90)

The Boltzmann factors $\exp(\check{\mu},y_{j})$ form the linearly independent
set. Therefore the natural way to meet these condition is: For any
$j=1,\ldots,q$

$\sum_{\rho\in R_{j}^{+}}\varphi_{\rho}-\sum_{\rho\in
R_{j}^{-}}\varphi_{\rho}=0$ (91)

This is the general complex balance condition. This condition is sufficient
for entropy growth, because it provides the equality $\theta(1)=\theta(0)$.

If we assume that $\varphi_{\rho}$ are constants or, for chemical kinetics,
depend only on temperature, then the conditions (91) give the general solution
to equation (90).

The complex balance condition is more general than the detailed balance.
Indeed, this is obvious: For the master equation (56) the complex balance
condition is trivially valid for all admissible constants. The first order
kinetics always satisfies the complex balance conditions. On the contrary, the
class of the master equations with detailed balance is rather special. The
dimension of the class of all master equations has dimension $n^{2}-n$
(constants for all transitions $A_{i}\to A_{j}$ are independent). For the
time-reversible Markov chains (the master equations with detailed balance)
there is only $n(n+1)/2-1$ independent constants: $n-1$ for equilibrium state
and $n(n-1)/2$ for transitions $A_{i}\to A_{j}$ ($i>j$), because for reverse
transitions the constant can be calculated through the detailed balance.

In general, for nonlinear reaction systems, the complex balance condition is
not necessary for entropy growth. In the next section we will give a more
general condition and demonstrate that there are systems that violate the
complex balance condition but satisfy this more general inequality.

####  6.2.4 6.2.4. $G$-Inequality

Gorban G1  proposed the following inequality for analysis of accordance
between thermodynamics and kinetics: $\theta(1)\geq\theta(0)$. This means that
for any values of $\check{\mu}$

$\sum_{\rho}\varphi_{\rho}\exp(\check{\mu},\alpha_{\rho})\geq\sum_{\rho}\varphi_{\rho}\exp(\check{\mu},\beta_{\rho})$
(92)

In the form of sum over complexes (similarly to (90)) it has the form

$\sum_{j=1}^{q}\exp(\check{\mu},y_{j})\left[\sum_{\rho\in
R_{j}^{+}}\varphi_{\rho}-\sum_{\rho\in R_{j}^{-}}\varphi_{\rho}\right]\geq 0$
(93)

Let us call these inequalities, (92), (93), the $G$-inequalities.

Here, two remarks are needed. First, functions $\exp(\check{\mu},y_{j})$ are
linearly independent but this does not allow us to transform inequalities (93)
similarly to (91) even for constant kinetic factors: Inequality between linear
combinations of independent functions may exist and the “simplified system”

$\sum_{\rho\in R_{j}^{+}}\varphi_{\rho}-\sum_{\rho\in
R_{j}^{-}}\varphi_{\rho}\geq 0\mbox{ for all }j$

is not equivalent to the $G$-inequality.

Second, this simplified inequality is equivalent to the complex balance
condition (with equality instead of $\geq$). Indeed, for any $\rho=1,\ldots,m$
there exist exactly one $j_{1}$ and one $j_{2}\neq j_{1}$ with properties:
$\rho\in R_{j_{1}}^{+}$, $\rho\in R_{j_{2}}^{-}$. Therefore, for any reaction
mechanism with reaction rates (74) the identity holds:

$\sum_{\rho}\left[\sum_{\rho\in R_{j}^{+}}\varphi_{\rho}-\sum_{\rho\in
R_{j}^{-}}\varphi_{\rho}\right]=0$

If all terms in this sum are non-negative then all of them are zeros.

Nevertheless, if at least one of the vectors $y_{j}$ is a convex combination
of others,

$\sum_{k,\,k\neq j}\lambda_{k}y_{k}=y_{j}\mbox{ for some }\lambda_{k}\geq
0,\,\sum_{k,\,k\neq j}\lambda_{k}=1$

then the $G$ inequality has more solutions than the condition of complex
balance. Let us take a very simple example with two components, $A_{1}$ and
$A_{2}$, three reactions and three complexes:

$2A_{1}\rightleftharpoons A_{1}+A_{2},\,2A_{2}\rightleftharpoons
A_{1}+A_{2},\,2A_{1}\rightleftharpoons 2A_{2}$
$y_{1}=(2,0),\,y_{2}=(0,2),\,y_{3}=(1,1)\,,$
$R_{1}^{+}=\\{1,3\\},\,R_{2}^{+}=\\{2,-3\\},\,R_{3}^{+}=\\{-1,-2\\}$
$R_{1}^{-}=\\{-1,-3\\},\,R_{2}^{-}=\\{-2,3\\},\,R_{3}^{-}=\\{1,2\\}$

The complex balance condition for this system is

$\begin{split}&(\varphi_{1}-\varphi_{-1})+(\varphi_{3}-\varphi_{-3})=0\,\\\
&(\varphi_{2}-\varphi_{-2})-(\varphi_{3}-\varphi_{-3})=0\end{split}$ (94)

The $G$-inequality for this system is

$\begin{split}&(\varphi_{1}+\varphi_{3}-\varphi_{-1}-\varphi_{-3})a^{2}+(\varphi_{2}+\varphi_{-3}-\varphi_{-2}-\varphi_{3})b^{2}\\\
&+(\varphi_{-1}+\varphi_{-2}-\varphi_{1}-\varphi_{2})ab\geq 0\mbox{ for all
}a,b>0\end{split}$ (95)

(here, $a,b$ stand for $\exp(\check{\mu}_{1}),\,\exp(\check{\mu}_{2})$). Let
us use for the coefficients at $a^{2}$ and $b^{2}$ notations $\psi_{a}$ and
$\psi_{b}$. Coefficient at $ab$ in (95) is $-(\psi_{a}+\psi_{b})$, linear
combinations $\psi_{a}=\varphi_{1}+\varphi_{3}-\varphi_{-1}-\varphi_{-3}$ and
$\psi_{b}=\varphi_{2}+\varphi_{-3}-\varphi_{-2}-\varphi_{3}$ are linearly
independent functions of variables $\varphi_{i}$ ($i=\pm 1,\pm 2,\pm 3$) and
we get the following task: To find all pairs of numbers
$(\psi_{a},\psi_{b})\in\mathbb{R}^{2}$ which satisfy the inequality

$\psi_{a}a^{2}+\psi_{b}b^{2}\geq(\psi_{a}+\psi_{b})ab\;\mbox{ for all }a,b>0$

Asymptotics $a\to 0$ and $b\to 0$ give $\psi_{a},\psi_{b}\geq 0$.

Let us use homogeneity of functions in (95), exclude one normalization factor
from $a,b$ and one factor from $\psi_{a},\psi_{b}$ and reduce the number of
variables: $b=1-a$, $\psi_{a}=1-\psi_{b}$: We have to find all such
$\psi_{b}\in[0,1]$ that for all $a\in]0,1[$

$a^{2}(1-\psi_{b})+(1-a)^{2}\psi_{b}-a(1-a)\geq 0$

The minimizer of this quadratic function of $a$ is
$a_{\min}=\frac{1}{4}+\frac{1}{2}\psi_{b}$, $a_{\min}\in]0,1[$ for all
$\psi_{b}\in[0,1]$. The minimal value is
$-2(\frac{1}{2}\psi_{b}-\frac{1}{4})^{2}$. It is nonnegative if and only if
$\psi_{b}=\frac{1}{2}$. When we return to the non-normalized variables
$\psi_{a},\psi_{b}$ then we get the general solution of the $G$-inequality for
this example: $\psi_{a}=\psi_{b}\geq 0$. For the kinetic factors this means:

$\begin{split}(\varphi_{1}-\varphi_{-1})+2(\varphi_{3}-\varphi_{-3})-(\varphi_{2}-\varphi_{-2})=0\,\\\
(\varphi_{1}-\varphi_{-1})+(\varphi_{3}-\varphi_{-3})\geq 0\,\\\
(\varphi_{2}-\varphi_{-2})-(\varphi_{3}-\varphi_{-3})\geq 0\end{split}$ (96)

These conditions are wider (weaker) than the complex balance conditions for
this example (94).

In the Stueckelberg language Stueckelberg1952 , the microscopic reasons for
the $G$-inequality instead of the complex balance (68) can be explained as
follows: Some channels of the scattering are unknown (hidden), hence, instead
of unitarity of $S$-matrix (conservation of the microscopic probability) we
have an inequality (the microscopic probability does not increase).

We can use other values of $\lambda_{0}\in[0,1[$ in inequality
$\theta(1)\geq\theta(\lambda_{0})$ and produce constructive sufficient
conditions of accordance between thermodynamics and kinetics. For example,
condition $\theta(1)\geq\theta(1/2)$ is weaker than $\theta(1)\geq\theta(0)$
because of convexity $\theta(\lambda)$.

One can ask a reasonable question: Why we do not use directly positivity of
entropy production ($\theta^{\prime}(1)\geq 0$) instead of this variety of
sufficient conditions. Of course, this is possible, but inequalities like
$\theta(1)\geq\theta(0)$ or equations like $\theta(1)=\theta(0)$ include
linear combinations of exponents of linear functions and often can be
transformed in algebraic equations or inequalities like in the example above.
Inequality $\theta^{\prime}(1)\geq 0$ includes transcendent functions like
$f\exp f$ (where $f$ is a linear function) which makes its study more
difficult.

##  7 Linear Deformation of Entropy

###  7.1 If Kinetics Does not Respect Thermodynamics then Deformation of
Entropy May Help

Kinetic equations in the general form (76) are very general, indeed. They can
be used for the approximation of any continuous time dynamical system on
compact $U$ Ocherki . In previous sections we demonstrated how to construct
the system in the form (76) with positivity of the entropy production when the
entropy function is given.

Let us consider a reverse problem. Assume that a system in the form (76) is
given but the entropy production is not always positive. How to find a new
entropy function for this system to guarantee the positivity of entropy
production?

Existence of such an entropy is very useful for analysis of stability of the
system. For example, let us take an arbitrary Mass Action Law system (80).
This is a rather general system with the polynomial right hand side. Its
stability or instability is not obvious a priori. It is necessary to check
whether bifurcations of steady states, oscillations and other interesting
effects of dynamics are possible for this system.

With the positivity of entropy productions these questions are much simpler
(for application of thermodynamic potentials to stability analysis see, for
example, Yablonskii1991 ; Ocherki ; VolpertKhudyaev1985 ; HangosAtAl1999 ). If
${\mathrm{d}}S/{\mathrm{d}}t\geq 0$ and it is zero only in steady states, then
any motion in compact $U$ converges to a steady state and all the non-
wandering points are steady states. (A non-wandering point is such a point
$x\in U$ that for any $T>0$ and $\varepsilon>0$ there exists such a motion
$c(t)\in U$ that (i) $\|c(0)-x\|<\varepsilon$ and (ii)
$\|c(T^{\prime})-x\|<\varepsilon$ for some $T^{\prime}>T$: A motion returns in
an arbitrarily small vicinity of $x$ after an arbitrarily long time.)
Moreover, the global maximizer of $S$ in $U$ is an asymptotically stable
steady state (at least, locally). It is a globally asymptotically stable point
if there is no other steady state in $U$.

For the global analysis of an arbitrary system of differential equations, it
is desirable either to construct a general Lyapunov function or to prove that
it does not exist. For the Lyapunov functions of the general form this task
may be quite difficult. Therefore, various finite-dimensional spaces of trial
functions are often in use. For example, quadratic polynomials of several
variables provide a very popular class of trial Lyapunov function.

In this section, we discuss the $n$-parametric families of Lyapunov functions
which are produced by the addition of linear function to the entropy:

$S(N)\mapsto S_{\Delta\check{\mu}}(N)=S(N)-\sum_{i}\Delta\check{\mu}_{i}N_{i}$
(97)

The change in potentials $\check{\mu}$ is simply the addition of
$\Delta\check{\mu}$:
$\check{\mu}_{i}\mapsto\check{\mu}_{i}+\Delta\check{\mu}_{i}$.

Let us take a general kinetic equation (76). We are looking for a
transformation that does not change the reaction rates. The Boltzmann factor
$\Omega_{\rho}=\exp(\check{\mu},\alpha_{\rho})$ transforms due to the change
of the entropy:
$\Omega_{\rho}\mapsto\Omega_{\rho}\exp(\Delta\check{\mu},\alpha_{\rho})$.
Therefore, to preserve the reaction rate, the transformation of the kinetic
factors should be
$\varphi_{\rho}\mapsto\varphi_{\rho}\exp(\Delta\check{\mu},\alpha_{\rho})$ in
order to keep the product $r_{\rho}=\Omega_{\rho}\varphi_{\rho}$ constant.

For the new entropy, $S=S_{\Delta\check{\mu}}$, with the new potential and
kinetic factors, the entropy production is given by (98):

$\begin{split}\frac{{\mathrm{d}}S}{{\mathrm{d}}t}=&-\sum_{\rho}(\gamma_{\rho},\check{\mu})\varphi_{\rho}\exp(\alpha_{\rho},\check{\mu})\\\
=&-\sum_{\rho}(\gamma_{\rho},\check{\mu}^{\rm
old}+\Delta\check{\mu})\varphi_{\rho}^{\rm
old}\exp(\alpha_{\rho},\check{\mu}^{\rm old})\\\ =&\frac{{\mathrm{d}}S^{\rm
old}}{{\mathrm{d}}t}-\sum_{\rho}(\gamma_{\rho},\Delta\check{\mu})\varphi_{\rho}^{\rm
old}\exp(\alpha_{\rho},\check{\mu}^{\rm old})\,\end{split}$ (98)

where the superscript “old” corresponds to the non-deformed quantities.

###  7.2 Entropy Deformation for Restoration of Detailed Balance

It may be very useful to find such a vector $\Delta\check{\mu}$ that in new
variables $\varphi_{\rho}^{+}=\varphi_{\rho}^{-}$. For the analysis of the
detailed balance condition, we group reactions in pairs of mutually inverse
reactions (85). Let us consider an equation of the general form (86) with
$r_{\rho}=r_{\rho}^{+}-r_{\rho}^{-}$, $\varphi_{\rho}^{\pm}>0$.

The problem is: To find such a vector $\Delta\check{\mu}$ that

$\varphi_{\rho}^{+}\exp{(\Delta\check{\mu},\alpha_{\rho})}=\varphi_{\rho}^{-}\exp{(\Delta\check{\mu},\beta_{\rho})}$
(99)

or, in the equivalent form of the linear equation

$(\Delta\check{\mu},\gamma_{\rho})=\ln\left(\frac{\varphi_{\rho}^{+}}{\varphi_{\rho}^{-}}\right)$
(100)

The necessary and sufficient conditions for the existence of such
$\Delta\check{\mu}$ are known from linear algebra: For every set of numbers
$a_{\rho}$ (${\rho}=1,\ldots,m$)

$\sum_{\rho}a_{\rho}\gamma_{\rho}=0\Rightarrow\sum_{\rho}a_{\rho}\ln\left(\frac{\varphi_{\rho}^{+}}{\varphi_{\rho}^{-}}\right)=0$
(101)

To check these conditions, it is sufficient to find a basis of solutions of
the uniform systems of linear equations

$\sum_{\rho}a_{\rho}\gamma_{\rho i}=0\;\;(i=1,\ldots,m)$

(that is, to find a basis of the left kernel of the matrix $\Gamma$, ${\rm
coim}\Gamma$, where $\Gamma=(\gamma_{\rho i})$) and then check for these basis
vectors the condition
$\sum_{\rho}a_{\rho}\ln\left(\frac{\varphi_{\rho}^{+}}{\varphi_{\rho}^{-}}\right)=0$
to prove or disprove that the vector with coordinates
$\ln\left(\frac{\varphi_{\rho}^{+}}{\varphi_{\rho}^{-}}\right)$ belongs to the
image of $\Gamma$, ${\rm im}\Gamma$.

For some of the reaction mechanisms it is possible to restore the detailed
balance condition for the general kinetic equation unconditionally. For these
reactions, for any set of positive kinetic factors, there exists such a vector
$\Delta\check{\mu}$ that the detailed balance condition (100) is valid for the
deformed entropy. According to (101) this means that there is no nonzero
solution $a_{\rho}$ for the equation $\sum_{\rho}a_{\rho}\gamma_{\rho}=0$. In
other words, vectors $\gamma_{\rho}$ are independent.

###  7.3 Entropy Deformation for Restoration of Complex Balance

The complex balance conditions (91) are, in general, weaker than the detailed
balance but they are still sufficient for the entropy growth.

Let us consider an equation of the general form (86). We need to find such a
vector $\Delta\check{\mu}$ that in new variables with the new entropy and
kinetic factors the complex balance conditions $\sum_{\rho\in
R_{j}^{+}}\varphi^{\rm new}_{\rho}-\sum_{\rho\in R_{j}^{-}}\varphi^{\rm
new}_{\rho}=0$ hold.

For our purpose, it is convenient to return to the presentation of reactions
as transitions between complexes. The complexes,
$\Theta_{1},\ldots,\Theta_{q}$ are the linear combinations,
$\Theta_{j}=(y_{j},A)$.

Each elementary reaction (17) with the reaction number $\rho$ may be
represented in the form $\Theta_{j}\to\Theta_{l}$, where $\Theta_{j}=\sum
y_{j}A_{j}$, $\rho\in R_{j}^{+}$ ($\alpha_{\rho}=y_{j}$) and $\rho\in
R_{j}^{-}$ ($\beta_{\rho}=y_{l}$). For this reaction, let us use the notation
$\varphi_{\rho}=\varphi_{lj}$. We used this notation in the analysis of
kinetics of compounds (Section 5.6). The complex balance conditions are

$\sum_{j,\,j\neq l}(\varphi_{lj}-\varphi_{jl})=0$ (102)

To obtain these conditions after the entropy deformation, we have to find such
$\Delta\check{\mu}$ that

$\sum_{j,\,j\neq
l}(\varphi_{lj}\exp{(\Delta\check{\mu},y_{j})}-\varphi_{jl}\exp{(\Delta\check{\mu},y_{l})})=0$
(103)

This is exactly the equation for equilibrium of a Markov chain with transition
coefficients $\varphi_{lj}$. Vector $(\Delta\check{\mu},y_{j})$ should be an
equilibrium state for this chain (without normalization to the unit sum of
coordinates).

For this finite Markov chain a graph representation is useful: Vertices are
complexes and oriented edges are reactions. To provide the existence of a
positive equilibrium we assume weak reversibility of the chain: If there
exists an oriented path from $\Theta_{j}$ to $\Theta_{l}$ then there exists an
oriented path from $\Theta_{l}$ to $\Theta_{j}$.

Let us demonstrate how to transform this problem of entropy deformation into a
linear algebra problem. First of all, let us find any positive equilibrium of
the chain, $\varsigma^{*}_{j}>0$:

$\sum_{j,\,j\neq
l}(\varphi_{lj}\varsigma^{*}_{j}-\varphi_{jl}\varsigma^{*}_{l})=0$ (104)

This is a system of linear equations. If we have already an arbitrary
equilibrium of the chain then other equilibria allow a very simple
description. We already found this description for kinetics of compounds (54)

Let us consider the master equation for the Markov chain with coefficients
$\varphi_{lj}$ and apply the formalism from Appendix 2:

$\frac{{\mathrm{d}}\varsigma}{{\mathrm{d}}t}=\sum_{j,\,j\neq
l}(\varphi_{lj}\varsigma_{j}-\varphi_{jl}\varsigma_{l})=0$ (105)

Let the graph of complex transformations $\Theta_{j}\to\Theta_{l}$ have $d$
connected components $C_{s}$ and let $V_{s}$ be the set of indexes of those
$\Theta_{j}$ which belong to $C_{s}$: $\Theta_{j}\in C_{s}$ if and only if
$j\in V_{s}$. For each $C_{s}$ there exists a conservation law
$\beta_{s}(\varsigma)$ for the master equation (53),
$\beta_{s}(\varsigma)=\sum_{j\in V_{s}}\varsigma_{j}$.

For any set of positive values of $\beta_{s}$ ($s=1,\ldots,q$) there exists a
unique equilibrium vector $\varsigma^{\rm eq}$ for (105) with this values
$\beta_{s}$ (54), (132). The set of equilibria is a linear space with the
natural coordinates $\beta_{s}$ ($s=1,\ldots,d$). We are interested in the
positive orthant of this space, $\beta_{s}>0$. For positive $\beta_{s}$,
logarithms of $\varsigma^{\rm eq}$ form a $d$-dimensional linear manifold in
$R^{q}$ (133). The natural coordinates on this manifold are $\ln\beta_{s}$.

Let us notice that the vector $\varsigma^{\circ}$ with coordinates

$\varsigma^{\circ}_{j}=\left(\frac{\varsigma_{j}^{*}(c,T)}{\sum_{l\in
V_{s}}\varsigma_{l}^{*}(c,T)}\right)\mbox{ for }j\in V_{s}$

is also an equilibrium for (105). This equilibrium is normalized to unit
values of all $\beta_{s}(\varsigma^{\circ})$. In the coordinates
$\ln\beta_{s}$ this is the origin. The equations for $\Delta\check{\mu}$ are

$(\Delta\check{\mu},y_{j})-\ln\beta_{s}=\ln\varsigma^{\circ}_{j}\mbox{ for
}j\in V_{s}$ (106)

This is a system of linear equations with respect to $n+d$ variables
$\Delta\check{\mu}_{i}$ ($i=1,\ldots,n$) and $\ln\beta_{s}$ ($s=1,\ldots,d$).
Let the coefficient matrix of this system be denoted by $\mathbf{M}$.

Analysis of solutions and solvability of such equations is one of the standard
linear algebra tasks. If this system has a solution then the complex balance
in the original system can be restored by the linear deformation of the
entropy. If this system is solvable for any right hand side, then for this
reaction mechanism we always can find the entropy, which provides the complex
balance condition.

Unconditional solvability of (106) means that the left hand side matrix of
this system has rank $q$. Let us express this rank through two important
characteristics: It is ${\rm rank}\\{\gamma_{1},\ldots,\gamma_{m}\\}+d$, where
$d$ is the number of connected components in the graph of transformation of
complexes.

To prove this formula, let us write down the matrix $\mathbf{M}$ of the system
(106). First, we change the enumeration of complexes. We group the complexes
from the same connected component together and arrange these groups in the
order of the connected component number. After this change of enumeration,
$\\{1,\ldots,|V_{1}|\\}=V_{1}$,
$\\{|V_{1}|+1,\ldots,|V_{1}|+|V_{2}|\\}=V_{2}$, …,
$\\{|V_{1}|+|V_{2}|+\ldots+|V_{d-1}|+1,\ldots,|V_{1}|+|V_{2}|+\ldots+|V_{d}|\\}=V_{d}$.

Let $y_{j}$ be here the row vector. The matrix is

$\mathbf{M}=\left[\begin{array}[]{lcccc}y_{1}&1&0&\ldots&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots\\\ y_{|V_{1}|}&1&0&\ldots&0\\\
y_{|V_{1}|+1}&0&1&\ldots&0\\\ \vdots&\vdots&\vdots&\vdots&\vdots\\\
y_{|V_{1}|+\ldots+|V_{d}|}&0&0&\ldots&1\end{array}\right]$ (107)

$\mathbf{M}$ consists of $d$ blocks $\mathbf{M}_{s}$, which correspond to
connected components $C_{s}$ of the graph of transformation of complexes:

$\mathbf{M}_{s}\left[\begin{array}[]{lcccc}y_{|V_{1}|+\ldots+|V_{s-1}|+1}&0&\ldots&1&\ldots\\\
\vdots&\vdots&\vdots&\vdots&\vdots\\\
y_{|V_{1}|+\ldots+|V_{s}|}&0&\ldots&1&\ldots\end{array}\right]$ (108)

The first $n$ columns in this matrix are filled by the vectors $y_{j}$ of
complexes, which belong to the component $C_{s}$, then follow $s-1$ columns of
zeros, after that, there is one column of units, and then again zeros. Here,
in (107), (108) we multiplied the last $d$ columns by $-1$. This operation
does not change the rank of the matrix.

Other elementary operations that do not change the rank are: We can add to any
row (column) a linear combination of other rows (columns).

We will use these operations to simplify blocks (108) but first we have to
recall several properties of spanning trees WuChao2004 . Let us consider a
connected, undirected graph $G$ with the set of vertices $\mathcal{V}$ and the
set of edges $\mathcal{E}\subset\mathcal{V}\times\mathcal{V}$. A spanning tree
of $G$ is a selection of edges of $G$ that form a tree spanning every vertex.
For a connected graph with $V$ vertices, any spanning tree has $V-1$ edges.
Let for each vertex $\Theta_{j}$ of $G$ a $n$-dimensional vector $y_{i}$ is
given. Then for every edge $(\Theta_{j},\Theta_{l})\in\mathcal{E}$ a vector
$\gamma_{jl}=y_{j}-y_{l}$ is defined. We identify vectors $\gamma$ and
$-\gamma$ and the order of $j,l$ is not important. Let us use $\Gamma_{G}$ for
this set of $\gamma_{jl}$:

$\Gamma_{G}=\\{y_{j}-y_{l}\,|\,(\Theta_{j},\Theta_{l})\in\mathcal{E}\\}$

For any spanning tree $T$ of graph $G$ we have the following property:

${\rm span}\Gamma_{G}={\rm span}\Gamma_{T}$ (109)

in particular, ${\rm rank}\Gamma_{G}={\rm rank}\Gamma_{T}$.

For the digraphs of reactions between complexes, we create undirected graphs
just by neglecting the directions of edges. We keep for them the same
notations as for original digraphs. Let us select any spanning tree $T_{s}$
for the connected component $C_{s}$ in the graph of transformation of
complexes. In $T_{s}$ we select arbitrarily a root complex. After that, any
other complex $\Theta_{j}$ in $C_{s}$ has a unique parent. This is the vertex
connected to it on the path to the root. For the root complex of $C_{s}$ we
use special notation $\Theta_{s}^{\circ}$.

Now, we transform the block (108) without change of rank: For each non-root
complex we subtract from the corresponding row the row which correspond to its
unique parent. After these transformations (and, may be, some permutations of
rows), the block $\mathbf{M}_{s}$ get the following form:

$\left[\begin{array}[]{lcccc}\gamma^{s}_{1}&0&\ldots&0&\ldots\\\
\gamma^{s}_{2}&0&\ldots&0&\ldots\\\ \vdots&\vdots&\vdots&\vdots&\vdots\\\
\gamma^{s}_{|V_{1}|-1}&0&\ldots&0&\ldots\\\
y_{s}^{\circ}&0&\ldots&1&\ldots\end{array}\right]$ (110)

Here, $\\{\gamma^{s}_{1},\gamma^{s}_{2},\ldots,\gamma^{s}_{|V_{1}|-1}\\}$ is
$\Gamma_{T_{s}}$ for the spanning tree $T_{s}$ and $y_{s}^{\circ}$ is the
coefficient vector for the root complex $\Theta_{s}^{\circ}$.

From the obtained structure of blocks we immediately find that the rank of the
rows with $\gamma$ is ${\rm rank}\\{\gamma_{1},\ldots,\gamma_{m}\\}+d$ due to
(109). Additional $d$ rows with $y_{s}^{\circ}$ are independent due to their
last coordinates and add $d$ to rank. Finally

${\rm rank}\mathbf{M}={\rm rank}\\{\gamma_{1},\ldots,\gamma_{m}\\}+d$ (111)

Obviously, ${\rm rank}\mathbf{M}\leq q$.

In particular, from the formula (111) immediately follows the description of
the reaction mechanisms, for which it is always possible to restore the
thermodynamic properties by the linear deformation of the entropy.

The deficiency zero theorem. If ${\rm rank}\mathbf{M}=q$ then it is always
possible to restore the positivity of the entropy production by the linear
deformation of the entropy.

Feinberg Feinberg1972  called the difference $q-{\rm rank}\mathbf{M}$ the
deficiency of the reaction network. For example, for the “Michaelis-Menten”
reaction mechanism $E+S\rightleftharpoons ES\rightleftharpoons P+S$ ${\rm
rank}\\{\gamma_{1},\gamma_{2}\\}=2$, $d=1$, $q=3$, ${\rm rank}\mathbf{M}=3$
and deficiency is 0.

For the adsorption (the Langmuir-Hinshelwood) mechanism of CO oxidation (38)
${{\rm rank}\\{\gamma_{1},\gamma_{2},\gamma_{3}\\}~{}=~{}3}$, $d=3$, $q=6$,
${\rm rank}\mathbf{M}=6$ and deficiency is 0. To apply the results about the
entropy deformation to this reaction mechanism, it is necessary to introduce
an inverse reaction to the third elementary reaction in (38),
PtO+PtCO$\to$CO2+2Pt with an arbitrarily small but positive constant in order
to make the mechanism weakly reversible.

Let us consider the Langmuir-Hinshelwood mechanism for reduced list of
components. Let us assume that the gas concentrations are constant because of
control or time separation or just as a model “fast” system and just include
them in the reaction rate constants for intermediates. Then the mechanism is
2Pt$\rightleftharpoons$2PtO, Pt$\rightleftharpoons$PtCO, PtO+PtCO$\to$2Pt. For
this system, ${\rm rank}\\{\gamma_{1},\gamma_{2},\gamma_{3}\\}=2$, $d=2$,
$q=5$, ${\rm rank}\mathbf{M}=4$ and deficiency is 1. Bifurcations in this
system are known Yablonskii1991 .

For the fragment of the reaction mechanism of the hydrogen combustion (44),
${\rm rank}\\{\gamma_{1},\ldots,\gamma_{m}\\}=6$, $d=7$, $q=16$, ${\rm
rank}\mathbf{M}=6+7=13$ and deficiency is 3.

###  7.4 Existence of Points of Detailed and Complex Balance

Our formulation of the conditions of detailed and complex balance is not
standard: We formulate them as the identities (87) and (89). These identities
have a global nature and describe the properties of reaction rates for all
states.

The usual approach to the principle of detailed balance is based on
equilibria. The standard formulation is: In all equilibria every process is
balanced with its reverse process. Without special forms of kinetic law this
principle cannot have any consequences for global dynamics. This is a trivial
but not widely known fact. Indeed, let a system $\dot{c}=F(c)$ be given in a
domain $U\subset\mathbb{R}^{n}$ and $\gamma_{1},\ldots,\gamma_{n}$ is an
arbitrary basis in $\mathbb{R}^{n}$. In this basis, we can always write:
$F(c)=\sum_{\rho=1}^{n}r_{\rho}(c)\gamma_{\rho}$. For any equilibrium $c^{*}$,
$r_{\rho}(c^{*})=0$. All the “reaction rates” $r_{\rho}(c)$ vanish
simultaneously. This “detailed balance” means nothing for dynamics because $F$
is an arbitrary vector field. Of course, if the system of vectors
$\\{\gamma_{\rho}\\}$ is not a basis but any complete system of vectors then
such “detailed balance” conditions, $r_{\rho}(c^{*})=0$, also do not imply any
specific features of dynamics without special hypotheses about functions
$r_{\rho}(c)$.

Nevertheless, if we fix the kinetic law then the consequences may be very
important. For example, if kinetics of elementary reactions follow the Mass
Action law then the existence of a positive equilibrium with detailed balance
implies existence of the Lyapunov function in the form of the perfect free
entropy:

$Y=-\sum_{i}c_{i}\left(\ln\left(\frac{c_{i}}{c^{*}_{i}}\right)-1\right)$

where $c_{i}^{*}$ is that positive equilibrium with detailed balance (see, for
example, Yablonskii1991 ).

In this section we demonstrate that for the general kinetic law (74), which
gives the expression of reaction rates through the entropy gradient, if the
kinetic factors are constant (or a function of temperature) then the existence
of the points of detailed (or complex) balance means that the linear
deformation of the entropy exists which restores the global detailed (or
complex) balance conditions (87) (or (89)).

The condition that the kinetic factors are constant means that for a given set
of values $\\{\varphi_{\rho}\\}$ a state with any admissible values of
$\check{\mu}$ is physically possible (admissible). This condition allows us to
vary the potentials $\check{\mu}$ independently of $\\{\varphi_{\rho}\\}$.

Let us assume that for the general kinetic system with the elementary reaction
rates given by (74) a point of detailed balance exists. This means that for
some value of $\check{\mu}=\check{\mu}^{*}$ (the detailed balance point in the
Legendre transform) and for all ${\rho}$ $r_{\rho}^{+}=r_{\rho}^{-}$:

$\varphi_{\rho}^{+}\exp(\alpha_{\rho},\check{\mu}^{*})=\varphi_{\rho}^{-}\exp(\alpha_{\rho},\check{\mu}^{*})$

This formula is exactly the condition (99) of existence of $\Delta\check{\mu}$
which allow us to deform the entropy for restoring the detailed balance in the
global form (87).

If we assume that the point of complex balance exists then there exists such a
value of $\check{\mu}=\check{\mu}^{*}$ (a point of complex balance in the
Legendre transform) that

$\sum_{j,\,j\neq
l}(\varphi_{lj}\exp{(\check{\mu}^{*},y_{j})}-\varphi_{jl}\exp{(\check{\mu}^{*},y_{l})})=0$

This is exactly the deformation condition (103) with
$\Delta\check{\mu}=\check{\mu}^{*}$.

To prove these statements we used an additional condition about possibility to
vary $\check{\mu}$ under given $\\{\varphi_{\rho}\\}$.

So, we demonstrated that for the general kinetic law (74) the existence of a
point of detailed balance is equivalent to the existence of such linear
deformation of the entropy that the global condition (87) holds. Analogously,
the existence of a point of complex balance is equivalent to the global
condition of complex balance after some linear deformation of the entropy.

###  7.5 The Detailed Balance is Needed More Often than the Complex Balance

The complex balance conditions are mathematically nice and more general than
the principle of the detailed balance. They are linked by Stueckelberg to the
Markov models (“$S$-matrix models”) of microscopic kinetics. Many systems
satisfy these conditions (after linear deformation of the entropy) just
because of the algebraic structure of the reaction mechanism (see Section
7.3). Nevertheless, it is used much less than the classical detailed balance.

The reason for the rare use of complex balance is simple: It is less popular
because the stronger condition, the principle of detailed balance, is valid
for most of physical and chemical systems. Onsager revealed the physical
reason for detailed balance Onsager1931a ; Onsager1931b . This is
microreversibility: The microscopic laws of motion are invertible in time: If
we observe the microscopic dynamics of particles in the backward movie then we
cannot find the difference from the real world. This difference occurs in the
macroscopic world.

In microphysics and the $S$-matrix theory this microreversibility property has
the name “$T$-invariance”.

Let us demonstrate how $T$-invariance in micro-world implies detailed balance
in macro-world.

Following Gibbs, we accept the ensemble-based point of view on the macroscopic
states: They are probability distributions in the space of detailed
microscopic states.

First of all, we assume that under given values of conservation laws
equilibrium state exists and is unique.

Second assumption is that the rates of elementary processes are
microscopically observable quantities. This means that somebody (a “demon”),
who observes all the events in the microscopical world can count the rates of
elementary reactions.

Because of $T$-invariance and uniqueness of equilibrium, the equilibrium is
$T$-invariant: If we change all the microscopic time derivatives (velocities)
$v$ to $-v$ then nothing will change.

$T$-transformation changes all reactions to the reverse reactions, just by
reversion of arrows, but the number of the events remains the same: Any
reaction transforms into its reverse reaction but does not change the reaction
rate. This can be formulated also as follows: $T$-transformation maps all
$r_{\rho}^{+}$ into the corresponding $r_{\rho}^{-}$.

Hence, because of the $T$-invariance, the equilibrium rate of each reaction is
equal to the equilibrium rate of the reverse reaction.

The violation of uniqueness of equilibrium for given values of conservation
laws seems improbable. Existence of several equilibria in thermodynamics is
quite unexpected for homogeneous systems but requires more attention for the
systems with phase separation. Nevertheless, if we assume that a multi-phase
system consists of several homogeneous phases, and each of these phases is in
uniform equilibrium, then we return to the previous assumption (with some
white spots for non-uniform interfaces).

$T$-invariance may be violated if the microscopic description is not
reversible in time. Magnetic field and the Coriolis force are the classical
examples for violation of the microscopic reversibility. In a linear
approximation near equilibrium the corresponding modification of the Onsager
relations give the Onsager-Casimir relations Casimir1945 . There are several
attempts for nonlinear formulation of the Onsager-Casimir relations (see
Grmela1993 ).

The principle of detailed balance seems to be still the best nonlinear version
of the Onsager relations for $T$-invariant systems, and the conditions of the
complex balance seem to give the proper relations between kinetic coefficients
in the absence of the microscopic reversibility for nonlinear systems. It is
important to mention here that all these relations are used together with the
general kinetic law (74).

Observability of the rates of elementary reactions deserves a special study.
Two approaches to the reaction rate are possible. If we accept that the
general kinetic law (74) is valid then we can find the kinetic factors by
observation of ${\mathrm{d}}c/{\mathrm{d}}t$ in several points because the
Boltzmann factors are linearly independent. In this sense, they are observable
but one can claim the approximation point of view and state that the general
kinetic law (74) without additional conditions on kinetic factors is very
general and allows to approximate any dynamical system. From this point of
view, kinetic coefficients are just some numbers in the approximation
algorithm and are not observable. This means that there is no such a