The Michaelis-Menten-Stueckelberg Theorem

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 ({\em 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.


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 [1]: 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 [2], but it was first claimed for a principle of informational statistical thermodynamics by Jaynes [3]. A very general discussion of the maximum entropy principle with applications to dissipative kinetics is given in the review [4]. Corrections of QE approximation with applications to physical and chemical kinetics were developed [5,6].
QSS was proposed by Bodenstein in 1913 [7], and the important Michaelis and Menten work [8] 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 [9][10][11]. 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) [12][13][14].
Wei and Prater [16] 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 = λl then d(l, c)/dt = (l, c)λ, 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 [15] and approximate [17] 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 [18,19], 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 [20]. 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 [21]. 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 [22].
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 [23,24].
In this paper, we focus on the combination of the QE approximation with the QSS approach.

The Structure of the Paper
Almost thirty years ago one of us published a book [25] 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 [8] but by Briggs and Haldane [12] 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 ⇋ 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 ⇋ 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 [33]. Therefore, kinetics of intermediates in heterogeneous catalysis may be nonlinear and demonstrate bifurcations, oscillations and other complex behavior.) In 1952, Stueckelberg [26] 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 v + w → v ′ + w ′ (v, w, v ′ , w ′ 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 [27].
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 → S 2 → . . . → S n → S 1 . Almost 100 years later, Cercignani and Lampis [28] 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 [29] and Coester [30] and the results are sometimes cited as the "Heitler-Coestler theorem of semi-detailed balance". In 1952, Stueckelberg [26] 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 [31]).
Later, these relations for the chemical mass action kinetics were rediscovered and called the complex balance conditions [51,63]. 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.

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: The linear combinations i α ρi A i and i β ρi A i are the complexes. For each complex 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: 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 → 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 ρ for the elementary reactions have a special form of the generalized mass action law (see Equation (74) below): where ϕ ρ > 0 is the kinetic factor and exp(α ρ ,μ) is the Boltzmann factor. Here and further in the text (α ρ ,μ) = i α ρiμi is the standard inner product, exp( , ) is the exponential of the value of the inner product andμ i are chemical potentials µ divided on RT .
For the prefect chemical systems,μ 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,μ 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 * : ϕ ρ = r ρ (c * ).
The linear kinetics of the compound reactions B i → B j implies the remarkable identity for the reaction rates, the complex balance condition (Equation (89) below) for all admissible values ofμ and given ϕ 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), for all ρ, are more restrictive than the complex balance conditions. For the perfect systems, the detailed balance condition takes the standard form: r + ρ (c * ) = r − ρ (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) 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 [63] q − d − rankΓ ≥ 0 where Γ = (γ ij ) = (β ij − α 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 [35] and by Evans and Polanyi [36].
Basic ideas behind the transition state theory are [37]: • 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).

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 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 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 [38,39] 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 ≈ x * M . The QE approximation defines the functions x * M as solutions to the following MaxEnt optimization problem: 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 is the QE entropy and the equation for the slow variables 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 dS * (M)/dt at point M according to the QE dynamics (7) and find dS(x)/dt at point x = x * M due to the initial system (3). The results always coincide: 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, dS/dt ≥ 0, then for the QE approximation (7) the production of the QE entropy is also non-negative, dS * /dt ≥ 0.
In addition, if for the initial system (dS/dt)| x = 0 if and only if F (x) = 0 then the same property holds in the QE approximation.

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 [3].
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 Here N is the vector of composition with components N i > 0, K sl corresponds to the slow reactions, K f s corresponds to fast reaction and ǫ > 0 is a small number. The system of fast reactions has the linear 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 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 [40].
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 [7] and was elaborated into an important tool for analysis of chemical reaction mechanisms and kinetics [9][10][11]. 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) [13,14]. In the enzyme kinetics, its invention was traditionally connected to the so-called Michaelis-Menten kinetics.

QSS and the Briggs-Haldane Asymptotic
Perhaps the first very clear explanation of the QSS was given by Briggs and Haldane in 1925 [12]. Briggs where the "Michaelis-Menten constant" is 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) [1,14,33] and a simple detailed case study [41]).
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 [33,34]). Let us demonstrate this on the traditional Michaelis-Menten system: 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 → 0 the reaction goes proportionally slower and to study this limit properly we have to adjust the time scale: dτ = e s dt: For small e/s, the second equation is a fast subsystem. According to this fast equation, for a given constant x, the variable ξ relaxes to ξ QSS = 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 [42,43] can be applied to the system in the form (14) and the QSS approximation is applicable even on an infinite time interval [44]. This transformation of variables and introduction of slow time is a standard procedure for rigorous proof of QSS validity in catalysis [33], enzyme kinetics [45] and other areas of kinetics and chemical engineering [13].
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.

The Michaelis and Menten Asymptotic
QSS is not QE but the classical work of Michaelis and Menten [8] 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 where [S i ] is the concentration of the ith 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 ith 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  (11) if we take k −1 ≫ k 2 in (11) (i.e., the equilibration S + E ⇌ SE is much faster than the reaction SE → 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.

Stoichiometric Algebra and Kinetic Equations
In this section, we introduce the basic notations of the chemical kinetics formalism. For more details see, for example, [33].
The list of components is a finite set of symbols A 1 , . . . , A n . A reaction mechanism is a finite set of the stoichiometric equations of elementary reactions: where ρ = 1, . . . , m is the reaction number and the stoichiometric coefficients α ρi , β ρi are nonnegative integers.
A stoichiometric vector γ ρ of the reaction (17) is a n-dimensional vector with coordinates that is, "gain minus loss" in the ρ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 A non-negative intensive quantity, r ρ , the reaction rate, corresponds to each reaction (17). The kinetic equations in the absence of external fluxes are If the volume is not constant then equations for concentrations includeV 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: where k ρ (T ) is the reaction rate constant. The generalized Arrhenius equation is: where R = 8.314 472 J K mol is the universal, or ideal gas constant, E aρ is the activation energy, S aρ is the activation entropy (i.e., E aρ − T S aρ is the activation free energy), A ρ is the constant pre-exponential factor. Some authors neglect the S aρ 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.

QE Manifold
In this section, we describe the general formalism of the QE for chemical kinetics following [34]. 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 ) = V f (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 ⊥ , 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 Let us use µ i as new coordinates. In these new coordinates (the "conjugated coordinates"), the QE manifold is just an orthogonal complement to L. This subspace, L ⊥ , is defined by equations It is sufficient to take in (23) not all γ ∈ 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 µ 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 ⊥ .
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: where b is the vector of slow variables, µ is the vector of chemical potentials and vectors γ ρ 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ḃ. 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 → 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 ⇋ 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 µ A /RT = ln c A + µ A0 , µ B /RT = ln c B +µ B0 , µ AB /RT = ln c AB +µ AB0 . This corresponds to the free energy F = V RT i c i (ln c i +µ i0 ), the correspondent free entropy (the Massieu-Planck potential) is −F/T . The stoichiometric vector is γ = (−1, −1, 1) and the equations (24) take the form From these equations we get the expressions for the QE concentrations: The QE free entropy is the value of the free entropy at this point, c(b 1 , b 2 ).

QE in Traditional MM System
Let us return to the simplest homogeneous enzyme reaction E + S ⇋ ES → P + S, the traditional Michaelis-Menten System (12) (it is simpler than the system studied by Michaelis and Menten [8]). Let us assume that the reaction E + S ⇋ ES is fast. This means that both k 1 and k −1 include large parameters: For small ǫ, 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 γ = (−1, −1, 1). The space L is one-dimensional and its basis is this vector γ. The space L ⊥ is two-dimensional and one of the convenient bases is 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): are the so-called standard equilibrium values and for perfect systems . Let us fix the slow variables and find c 1,2,3 . Equations (24) turn into 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 : Therefore, The sign "−" is selected to provide positivity of all c i . This choice provides also the proper asymptotic: The time derivatives of concentrations are: here we added external flux with input and output velocities (per unite volume) v in and v out and input concentrations c 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ḃ 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 [8]): Under this assumption, the quadratic equation (25) transforms into and in this approximation (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 forṁ This is the typical form in the reduced equations for catalytic reactions: Nominator in the reaction rate corresponds to the "brutto reaction" S + E → P + E [33,49].

Heterogeneous Catalytic Reaction
For the second example, let us assume equilibrium with respect to the adsorption in the CO on Pt oxidation: (for detailed discussion of the modeling of CO on Pt oxidation, this "Mona Liza" of catalysis, we address readers to [33]). The list of components involved in these 2 reactions is: does not participate in adsorption and may be excluded at this point).
The corresponding slow variables are For heterogeneous systems, caution is needed in transition between N and c variables because there are two "volumes" and we cannot put in (33) There is a law of conservation of the catalyst: N Pt + N PtO + N 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 O 2 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 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 The QE manifold in the conjugated variables is given by equations: It is trivial to resolve these equations with respect to µ 3,4 , for example: or with the standard equilibria: or in the kinetic form (we assume that the kinetic constants are in accordance with thermodynamics and all these forms are equivalent): The next task is to solve the system of equations: 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): The kinetic equations for this system (including the flux in the gas phase) is Here v in and v out are the flux rates (per unit volume). For the slow variables this equation gives: This system looks quite simple. Only one reaction, 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: Together with this smallness assumptions equations (37) give: In this approximation, we have for the reaction (41) rate This expression gives the closure for the slow QE equations (40).

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 [47].
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 [33]. 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.

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 Θ 1 , . . . , Θ q . The number of complexes q ≤ 2m (two complexes per elementary reaction, as the maximum) and it is possible that q < 2m because some complexes may coincide for different reactions.
where ν i is a vector with coordinates ν ij . Each elementary reaction (17) may be represented in the form Θ − ρ → Θ + ρ , where Θ ± ρ 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, [48] There are 16 different complexes here: The reaction set (44) can be represented as 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.

Stoichiometry of Compounds
For each complex Θ j we introduce an additional component B j , an intermediate compound and B ± ρ are those compounds B j (1 ≤ j ≤ 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 [8] and keep "complexes" for the formal linear combinations Θ j .
An extended reaction mechanism includes two types of reactions: Equilibration between a complex and its compound (q reactions, one for each complex) and transformation of compounds B − ρ → B + ρ (m reactions, one for each elementary reaction from (17). So, instead of the reaction (17) we can write 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 ⊕ is a direct sum of two vectors, the composition vector for initial species, N, with coordinates N i (i = 1, . . . , n) and the composition vector for compounds, Υ, with coordinates Υ j (j = 1, . . . , q): N ⊕ = N ⊕ Υ.
The space of composition vectors E is a direct sum of n-dimensional E A and q-dimensional E B : For concentrations of A i we use the notation c i and for concentrations of B j we use ς j . The stoichiometric vectors for reactions Θ j ⇋ B j (45) are direct sums: g j = −ν j ⊕ e j , where e j is the jth standard basis vector of the space R q = E B , the coordinates of e j are e jl = δ jl : The stoichiometric vectors of equilibration reactions (45) are linearly independent because there exists exactly one vector for each l.
The stoichiometric vectors γ jl of reactions B j → B l belong entirely to E B . They have jth coordinate −1, lth 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 Θ i and edges, which correspond to reactions from (17). The system is weakly reversible if for any two vertices Θ i and Θ j , the existence of an oriented path from Θ i to Θ j implies the existence of an oriented path from Θ j to Θ i .
Of course, this weak reversibility property is equivalent to weak reversibility of the reaction network between compounds B j .

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 ∞ [50]. 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: Let us introduce the standard equilibrium concentrations for B j . Due to the Boltzmann distribution (exp(−u/RT )) and formula (48) where 1/Z is the normalization factor. Let us select here the normalization Z = 1 and write: We assume that the standard equilibrium concentrations ς * 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 0 i − aν ji , ς j = a. If we neglect the products ς j ∂ς * j (c, T )/∂c i as the second order small quantities then the minimizers have the very simple form: or where is the chemical potential of A i and ∂ς j is the chemical potential of B j ). The thermodynamic equilibrium of the system of reactions B j → B l that corresponds to the reactions (46) The positive values of concentrations ς j are the equilibrium concentrations (54) for some values of β s if and only if for any s = 1, . . . , d and all j, l ∈ V s ϑ j = ϑ l (55) (ϑ j = ln(ς j /ς * 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 µ i /RT (i = 1, . . . , n) and ϑ j (j = 1, . . . , 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).

Markov Kinetics of Compounds
For the kinetics of compounds transformations B j → B l , the same hypothesis of the smallness of concentrations leads to the only reasonable assumption: The linear (monomolecular) kinetics with the rate constant κ lj > 0. This "constant" is a function of c, T : κ lj (c, T ). The order of indexes at κ is inverse to the order of them in reaction: κ lj = κ l←j .
The master equation for the concentration of B j gives: 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 ς * (49) should be an equilibrium for (56), that is, for every j = 1, . . . , q l, l =j This condition is necessary because the standard equilibrium is the free energy minimizer for given c, T and j ς j = j ς * j . The sum j ς 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 j ς 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).

Thermodynamics and Kinetics of the Extended System
In this section, we consider the complete extended system, which consists of species A i (i = 1, . . . , n) and compounds B j (j = 1, . . . , q) and includes reaction of equilibration (45) and transformations of compounds B j → 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 Θ j ⇋ B j (where the complex is a formal combination: Θ j = i ν ji A i ) we introduce the reaction rate w j . This rate should be positive if and negative if 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 If several dynamical systems defined by equationsẋ = J 1 , ...ẋ = J v on the same space have the same Lyapunov function F , then for any conic combination J = k a k J k (a k ≥ 0, k a k > 0) the dynamical systemẋ = J also has the Lyapunov function F . The free energy (50) decreases monotonically due to any reaction Θ j ⇋ 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: where the coefficients κ 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.

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 ⊥ , 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 ⊥ is equal to the number of components A i : dim L ⊥ = n. For each A i we will find a vector b i ∈ L ⊥ that has the following first n coordinates: b ik = δ ik for k = 1, . . . , n. The condition (b i , g j ) = 0 gives immediately: b i,n+j = ν ji . Finally, The corresponding slow variables are In the QE approximation all w j = 0 and the kinetic equations (61) give in this approximation In these equations, we have to use the dependence ς(b). Here we use the QSS Michaelis and Menten assumption: The compounds are present in small amounts In this case, we can take b i instead of c i (i.e., take µ(b, T ) instead of µ(c, T )) in the formulas for equilibria (52): In the final form of the QE kinetic equation there remain two "offprints" of the compound kinetics: Two sets of functions ς * j (b, T ) ≥ 0 and κ jl (b, T ) ≥ 0. These functions are connected by the identity (57). The final form of the equations is The identity (57), l, l =j κ jl ς * l = l, l =j κ lj ς * 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) Θ j ⇋ B j do not produce entropy and all changes in the total free energy are caused by the Markov kinetics B i → 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 Θ j ⇋ 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 ς * j (b, T ) and κ 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: We call this function ϕ jl (b, T ) the kinetic factor. The identity (57) for the kinetic factor is l, l =j We call the thermodynamic factor (or the Boltzmann factor) the second multiplier in the reaction rates In this notation, the kinetic equations (66) have a simple form 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 Θ l , Θ j (46). It has the form Θ l → Θ j and the reaction rate is r = ϕ jl Ω l . In the right hand side in (70) this reaction appears twice: first time with sign "+" and the vector coefficient ν j and the second time with sign "−" and the vector coefficient ν l . The stoichiometric vector of this reaction is γ = ν j − ν l . Let us enumerate the elementary reactions by index ρ, which corresponds to the pair (j, l). Finally, we transform (46) into the sum over reactions form In the vector form it looks as follows:

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): under the following asymptotic assumptions: 1. Concentrations of the compounds B ρ are close to their quasiequilibrium values (65) (this may be due to the fast reversible reactions in (46)); 2. Concentrations of the compounds B ρ are much smaller than the concentrations of the components A i : There is a small positive parameter ε ≪ 1, ς * j = εξ * j and ξ * j do not depend on ε; 3. Kinetics of transitions between compounds B i → B j is linear (Markov) kinetics with reaction rate constants k ji = 1 ε κ ji .
Under these assumptions, in the asymptotic δ, ε → 0, δ, ε > 0 kinetics of components A i may be described by the reaction mechanism where the kinetic factors ϕ ρ satisfy the condition (68): for any vector v from the set of all vectors {α ρ , β ρ }. This statement includes the generalized mass action law for r ρ and the balance identity for kinetic factors that is sufficient for the entropy growth as it is shown in the next Section 6.

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 , . . . , 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 [52]) 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 B p ln p).
The dual variables, potentials, are defined as the partial derivatives of S: 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) 1, . . . , m). In general, there is no need to assume that the stoichiometric coefficients α ρi , β ρi are integers.
The assumption that they are nonnegative, α ρi ≥ 0, β ρi ≥ 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 [53]. In this case, the coefficient α 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 α ρi ≥ 0, β ρi ≥ 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 γ ρ = β ρ − α ρ (the "gain minus loss" vector).
For each reaction, a nonnegative quantity, reaction rate r ρ is defined. We assume that this quantity has the following structure: where (α ρ ,μ) = i α ρiμ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 ρ (some of them usually coincide), and we can write 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 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 goes in the same volume as the reverse reaction If this condition (77) holds then the detailed balance conditions and the complex balance conditions will hold separately in all volumes V ρ . An important particular case of (76) gives us the Mass Action Law. Let us take the perfect free entropy For the perfect function (78) and for the reaction rate function (74) we get The standard assumption for the Mass Action Law in physics and chemistry is that ϕ and c * are functions of temperature: ϕ ρ = ϕ ρ (T ) and c * i = c * i (T ). To return to the kinetic constants notation (20) we should write: (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. The energy of the Universe is constant. 2. The entropy of the Universe tends to a maximum. (R. Clausius, 1865 [54].) 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 dl(N) where v in,out are the input and output fluxes per unit volume, l 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: 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 γ ρ . 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).

General Entropy Production Formula
Let us calculate dS/dt due to equations (76): An auxiliary function θ(λ) of one variable λ ∈ [0, 1] is convenient for analysis of dS/dt (it was studied by Rozonoer and Orlov [55], see also [25]: With this function, the entropy production (82) has a very simple form: The auxiliary function θ(λ) allows the following interpretation. Let us introduce the deformed stoichiometric mechanism with the stoichiometric vectors, , which is the initial mechanism when λ = 1, the inverted mechanism with interchange of α and β when λ = 0, the trivial mechanism (the left and right hand sides of the stoichiometric equations coincide) when λ = 1/2.
For the deformed mechanism, let us take the same kinetic factors and calculate the Boltzmann factors with α ρ (λ): r ρ (λ) = ϕ ρ exp(α ρ (λ),μ) In this notation, the auxiliary function θ(λ) is a sum of reaction rates for the deformed reaction mechanism: In particular, θ(1) = ρ r ρ , this is just the sum of reaction rates. Function θ(λ) is convex. Indeed This convexity gives the following necessary and sufficient condition for positivity of entropy production: dS dt > 0 if and only if θ(λ) < θ(1) for some λ < 1 In several next subsections we study various important particular sufficient conditions for positivity of entropy production.

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 [27].
Let us join elementary reactions in pairs: 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: In this notation, the principle of detailed balance is very simple: The thermodynamic equilibrium in the direction γ ρ , given by the standard condition (γ ρ ,μ) = 0, is equilibrium for the corresponding pair of mutually reverse reactions from (85). For kinetic factors this transforms into the simple and beautiful condition: For the systems with detailed balance we can take ϕ ρ = ϕ + ρ = ϕ − ρ and write for the reaction rate: M. Feinberg called this kinetic law the "Marselin-De Donder" kinetics [46]. This representation of the reaction rates gives for the auxiliary function θ(λ): Each term in this sum is symmetric with respect to change λ → (1 − λ). Therefore, θ(1) = θ(0) and, because of convexity of θ(λ), θ ′ (1) ≥ 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 [56,57]. Interrelations between positivity of entropy production, Onsager reciprocal relations and detailed balance were analyzed in detail by N.G. van Kampen [58].

Complex Balance
The principle of detailed balance gives us θ(1) = θ(0) and this equality holds for each pair of mutually reverse reactions.
Let us start now from the equality θ(1) = θ(0). We return to the initial stoichiometric equations (17) without joining the direct and reverse reactions. The equality reads Exponential functions exp(μ, y) form linearly independent family in the space of functions ofμ 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(μ, y). Here we have to return to the complex-based representation of reactions (see Section 5.1). Let us consider the family of vectors {α ρ , β ρ } (ρ = 1, . . . , m). Usually, some of these vectors coincide. Assume that there are q different vectors among them. Let y 1 , . . . , y q be these vectors. For each j = 1, . . . , q we take We can rewrite the equality (89) in the form q j=1 exp(μ, y j ) The Boltzmann factors exp(μ, y j ) form the linearly independent set. Therefore the natural way to meet these condition is: For any j = 1, . . . , q This is the general complex balance condition. This condition is sufficient for entropy growth, because it provides the equality θ(1) = θ(0).
If we assume that ϕ ρ 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 → 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 → 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.

G-Inequality
Gorban [25] proposed the following inequality for analysis of accordance between thermodynamics and kinetics: θ(1) ≥ θ(0). This means that for any values ofμ In the form of sum over complexes (similarly to (90)) it has the form q j=1 exp(μ, y j ) Let us call these inequalities, (92), (93), the G-inequalities.
Here, two remarks are needed. First, functions exp(μ, 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" is not equivalent to the G-inequality.
Second, this simplified inequality is equivalent to the complex balance condition (with equality instead of ≥). Indeed, for any ρ = 1, . . . , m there exist exactly one j 1 and one j 2 = j 1 with properties: ρ ∈ R + j 1 , ρ ∈ R − j 2 . Therefore, for any reaction mechanism with reaction rates (74) the identity holds: 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, k, k =j λ k y k = y j for some λ k ≥ 0, k, k =j λ 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: The complex balance condition for this system is The G-inequality for this system is (here, a, b stand for exp(μ 1 ), exp(μ 2 )). Let us use for the coefficients at a 2 and b 2 notations ψ a and ψ b . Coefficient at ab in (95) is −(ψ a + ψ b ), linear combinations ψ a = ϕ 1 + ϕ 3 − ϕ −1 − ϕ −3 and ψ b = ϕ 2 + ϕ −3 − ϕ −2 − ϕ 3 are linearly independent functions of variables ϕ i (i = ±1, ±2, ±3) and we get the following task: To find all pairs of numbers (ψ a , ψ b ) ∈ R 2 which satisfy the inequality Asymptotics a → 0 and b → 0 give ψ a , ψ b ≥ 0. Let us use homogeneity of functions in (95), exclude one normalization factor from a, b and one factor from ψ a , ψ b and reduce the number of variables: b = 1 − a, ψ a = 1 − ψ b : We have to find all such ψ b ∈ [0, 1] that for all a ∈]0, 1[ The minimizer of this quadratic function of a is a min = 1 It is nonnegative if and only if ψ b = 1 2 . When we return to the non-normalized variables ψ a , ψ b then we get the general solution of the G-inequality for this example: ψ a = ψ b ≥ 0. For the kinetic factors this means: These conditions are wider (weaker) than the complex balance conditions for this example (94).
In the Stueckelberg language [26], 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).
One can ask a reasonable question: Why we do not use directly positivity of entropy production (θ ′ (1) ≥ 0) instead of this variety of sufficient conditions. Of course, this is possible, but inequalities like θ(1) ≥ θ(0) or equations like θ(1) = θ(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 θ ′ (1) ≥ 0 includes transcendent functions like f exp f (where f is a linear function) which makes its study more difficult.

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 [59]. 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, [33,[59][60][61]). If dS/dt ≥ 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 ∈ U that for any T > 0 and ε > 0 there exists such a motion c(t) ∈ U that (i) c(0) − x < ε and (ii) c(T ′ ) − x < ε for some T ′ > 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: The change in potentialsμ is simply the addition of ∆μ:μ i →μ i + ∆μ 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 Ω ρ = exp(μ, α ρ ) transforms due to the change of the entropy: Ω ρ → Ω ρ exp(∆μ, α ρ ). Therefore, to preserve the reaction rate, the transformation of the kinetic factors should be ϕ ρ → ϕ ρ exp(∆μ, α ρ ) in order to keep the product r ρ = Ω ρ ϕ ρ constant.
For the new entropy, S = S ∆μ , with the new potential and kinetic factors, the entropy production is given by (98): where the superscript "old" corresponds to the non-deformed quantities.

Entropy Deformation for Restoration of Detailed Balance
It may be very useful to find such a vector ∆μ that in new variables ϕ + ρ = ϕ − ρ . 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 ρ = r + ρ − r − ρ , ϕ ± ρ > 0. The problem is: To find such a vector ∆μ that or, in the equivalent form of the linear equation The necessary and sufficient conditions for the existence of such ∆μ are known from linear algebra: For every set of numbers a ρ (ρ = 1, . . . , m) To check these conditions, it is sufficient to find a basis of solutions of the uniform systems of linear equations ρ a ρ γ ρi = 0 (i = 1, . . . , m) (that is, to find a basis of the left kernel of the matrix Γ, coimΓ, where Γ = (γ ρi )) and then check for these basis vectors the condition ρ a ρ ln ϕ + ρ ϕ − ρ = 0 to prove or disprove that the vector with coordinates belongs to the image of Γ, imΓ. 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 ∆μ that the detailed balance condition (100) is valid for the deformed entropy. According to (101) this means that there is no nonzero solution a ρ for the equation ρ a ρ γ ρ = 0. In other words, vectors γ ρ are independent.

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 ∆μ that in new variables with the new entropy and kinetic factors the complex balance conditions For our purpose, it is convenient to return to the presentation of reactions as transitions between complexes. The complexes, Θ 1 , . . . , Θ q are the linear combinations, Θ j = (y j , A).
Each elementary reaction (17) with the reaction number ρ may be represented in the form Θ j → Θ l , where Θ j = y j A j , ρ ∈ R + j (α ρ = y j ) and ρ ∈ R − j (β ρ = y l ). For this reaction, let us use the notation ϕ ρ = ϕ lj . We used this notation in the analysis of kinetics of compounds (Section 5.6). The complex balance conditions are To obtain these conditions after the entropy deformation, we have to find such ∆μ that j, j =l (ϕ lj exp (∆μ, y j ) − ϕ jl exp (∆μ, y l )) = 0 This is exactly the equation for equilibrium of a Markov chain with transition coefficients ϕ lj . Vector (∆μ, 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 Θ j to Θ l then there exists an oriented path from Θ l to Θ 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, ς * j > 0: 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 ϕ lj and apply the formalism from Appendix 2: Let the graph of complex transformations Θ j → Θ l have d connected components C s and let V s be the set of indexes of those Θ j which belong to C s : Θ j ∈ C s if and only if j ∈ V s . For each C s there exists a conservation law β s (ς) for the master equation (53) Let us notice that the vector ς • with coordinates is also an equilibrium for (105). This equilibrium is normalized to unit values of all β s (ς • ). In the coordinates ln β s this is the origin. The equations for ∆μ are (∆μ, y j ) − ln β s = ln ς • j for j ∈ V s (106) This is a system of linear equations with respect to n + d variables ∆μ i (i = 1, . . . , n) and ln β s (s = 1, . . . , d). Let the coefficient matrix of this system be denoted by 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 rank{γ 1 , . . . , γ 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 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, . . . , Let y j be here the row vector. The matrix is 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 [62]. Let us consider a connected, undirected graph G with the set of vertices V and the set of edges E ⊂ V × 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 Θ j of G a n-dimensional vector y i is given. Then for every edge (Θ j , Θ l ) ∈ E a vector γ jl = y j − y l is defined. We identify vectors γ and −γ and the order of j, l is not important. Let us use Γ G for this set of γ jl : For any spanning tree T of graph G we have the following property: in particular, rankΓ G = rankΓ 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 Θ 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 Θ • s . 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 M s get the following form: Here, {γ s 1 , γ s 2 , . . . , γ s |V 1 |−1 } is Γ Ts for the spanning tree T s and y • s is the coefficient vector for the root complex Θ • s . From the obtained structure of blocks we immediately find that the rank of the rows with γ is rank{γ 1 , . . . , γ m } + d due to (109). Additional d rows with y • s are independent due to their last coordinates and add d to rank. Finally Obviously, rankM ≤ 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 rankM = q then it is always possible to restore the positivity of the entropy production by the linear deformation of the entropy.
For the adsorption (the Langmuir-Hinshelwood) mechanism of CO oxidation (38) rank{γ 1 , γ 2 , γ 3 } = 3, d = 3, q = 6, rankM = 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→CO 2 +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⇋2PtO, Pt⇋PtCO, PtO+PtCO→2Pt. For this system, rank{γ 1 , γ 2 , γ 3 } = 2, d = 2, q = 5, rankM = 4 and deficiency is 1. Bifurcations in this system are known [33].

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ċ = F (c) be given in a domain U ⊂ R n and γ 1 , . . . , γ n is an arbitrary basis in R n . In this basis, we can always write: F (c) = n ρ=1 r ρ (c)γ ρ . For any equilibrium c * , r ρ (c * ) = 0. All the "reaction rates" r ρ (c) vanish simultaneously. This "detailed balance" means nothing for dynamics because F is an arbitrary vector field. Of course, if the system of vectors {γ ρ } is not a basis but any complete system of vectors then such "detailed balance" conditions, r ρ (c * ) = 0, also do not imply any specific features of dynamics without special hypotheses about functions r ρ (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: where c * i is that positive equilibrium with detailed balance (see, for example, [33]). 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 {ϕ ρ } a state with any admissible values ofμ is physically possible (admissible). This condition allows us to vary the potentialsμ independently of {ϕ ρ }.
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μ =μ * (the detailed balance point in the Legendre transform) and for all ρ r + ρ = r − ρ : This formula is exactly the condition (99) of existence of ∆μ 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μ =μ * (a point of complex balance in the Legendre transform) that j, j =l (ϕ lj exp (μ * , y j ) − ϕ jl exp (μ * , y l )) = 0 This is exactly the deformation condition (103) with ∆μ =μ * .
To prove these statements we used an additional condition about possibility to varyμ under given {ϕ ρ }.
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.

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 [56,57]. 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 + ρ into the corresponding r − ρ . 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 [64]. There are several attempts for nonlinear formulation of the Onsager-Casimir relations (see [65]). 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 dc/dt 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 microscopic thing as the rate of elementary reaction, and the set of reactions serves just for the approximation of the right hand side of the kinetic equation. We cannot fully disprove this point of view but can just say that in some cases the collision-based approach with physically distinguished elementary reactions is based on the solid experimental and theoretical background. If the elementary reactions physically exist then the detailed balance for T -invariant systems is proved.

Conclusions
We present the general formalism of the Quasiequilibrium approximation (QE) with the proof of the persistence of entropy production in the QE approximation (Section 2).
We demonstrate how to apply this formalism to chemical kinetics and give several examples for the Mass Action law kinetic equation. We discuss the difference between QE and Quasi-Steady-State (QSS) approximations and analyze the classical Michaelis-Menten and Briggs-Haldane model reduction approaches (Section 3). After that, we use ideas of Michaelis, Menten and Stueckelberg to create a general approach to kinetics.
Let us summarize the main results of our discussion. First of all, we believe that this is the finish of the Michaelis-Menten-Stueckelberg program. The approach to modeling of the reaction kinetics proposed by Michaelis and Menten in 1913 [8] for enzyme reactions was independently in 1952 applied by Stueckelberg [26] to the Boltzmann equation.
The idea of the complex balance (cyclic balance) relations was proposed by Boltzmann as an answer to the Lorentz objections against Boltzmann's proof of the H-theorem. Lorentz stated that the collisions of the polyatomic molecules may have no inverse collisions. Cercignani and Lampis [28] demonstrated that the Boltzmann H-theorem based on the detailed balance conditions is valid for the polyatomic molecules under the microreversibility conditions and this new Boltzmann's idea was not needed.
Nevertheless, this seminal idea was studied further by many authors [29][30][31] mostly for linear systems. Stueckelberg [26] proved these conditions for the Boltzmann equation. He used in his proof the S-matrix representation of the micro-kinetics.
Some consequences of the Stueckelberg approach were rediscovered for the Mass Action law kinetics by Horn and Jackson in 1972 [51] and supplemented by the "zero deficiency theorem" [63]. This is the history.
In our work, we develop the Michaelis-Menten-Stueckelberg approach to general kinetics. This is a combination of the QE (fast equilibria) and the QSS (small amounts) approaches to the real or hypothetical intermediate states. These intermediate states (compounds) are included in all elementary reactions (46) as it is illustrated in Figure 1. Because of the small amount, the free energy for these compounds B i is perfect (48), the kinetics of compounds is the first order Markov kinetics and satisfies the master equation.
After that, we use the combination of QE and QSS approximations and exclude the concentrations of compounds. For the general kinetics the main result of this approach is the general kinetic law (74). Earlier, we just postulated this law because of its convenient and natural form [25,53], now we have the physical framework where this law can be proved.
We do not assume anything about reaction rates of the main reactions (17). We use only thermodynamic equilibrium, the hypothesis about fast equilibrium with compounds and the smallness of concentration of compounds. This smallness implies the perfect entropy and the first order kinetics for compounds. After that, we get the reaction rate functions from the qualitative assumptions about compounds and the equilibrium thermodynamic data.
For example, if we relax the assumption about fast equilibrium and use just smallness of compound concentrations (the Briggs-Haldane QSS approach [12][13][14]) then we immediately need the formulas for reaction rates of compound production. Equilibrium data become insufficient. If we relax the assumption about smallness of concentrations then we lose the perfect entropy and the first order Markov kinetics. So, only the combination of QE and QSS gives the desired result.
For the kinetics of rarefied gases the mass action law for elastic collisions (the Boltzmann equation) or for inelastic processes like chemical reactions follows from the "molecular chaos" hypothesis and the low density limits. The Michaelis-Menten-Stueckelberg approach substitutes low density of all components by low density of the elementary events (or of the correspondent compounds) together with the QE assumption.
The general kinetic law has a simple form: For an elementary reaction the reaction rate is r = ϕΩ, where Ω > 0 is the Boltzmann factor, Ω = exp ( i α iμi ),μ i = −∂S/∂N i is the chemical potential µ divided by RT , and ϕ ≥ 0 is the kinetic factor. Kinetic factors for different reactions should satisfy some conditions. Two of them are connected to the basic physics: • The detailed balance: The kinetic factors for mutually reverse reaction should coincide, ϕ + = ϕ − .
This identity is proven for systems with microreversibility (Section 7.5). For the general kinetic law we studied several sufficient conditions of accordance between thermodynamics and kinetics: Detailed balance, complex balance and G-inequality.
In the practice of modeling, a kinetic model may, initially, do not respect thermodynamic conditions. For these cases, we solved the problem of whether it is possible to add a linear function to entropy in order to provide agreement with the given kinetic model and deformed thermodynamics. The answer is constructive (Section 7) and allows us to prove the general algebraic conditions for the detailed and complex balance.
Finally, we have to mention that Michaelis, Menten and Stueckelberg did not prove their "big theorem". Michaelis and Menten did not recognize that their beautiful result of mass action law produced from the equilibrium relations between substrates and compounds, the assumption about smallness of compound concentrations and the natural hypothesis about linearity of compound kinetics is a general theorem. Stueckelberg had much more and fully recognized that his approach decouples the Boltzmann H-theorem and the microreversibility (detailed balance). This is important because for every professional in theoretical physics it is obvious that the microreversibility cannot be important necessary condition for the H-theorem. Entropy production should be positive without any relation to detailed balance (the proof of the H-theorem for systems with detailed balance is much simpler but it does not matter: Just the Markov microkinetics is sufficient for it). Nevertheless, Stueckelberg did not produce the generalized mass action law and did not analyze the general kinetic equation. Later, Horn, Jackson and Feinberg approached the complex balance conditions again and studied the generalized mass action law but had no significant interest in the microscopic assumptions behind these properties. Therefore, this paper is the first publication of the Michaelis-Menten-Stueckelberg theorem.