Quasi-steady-state and partial-equilibrium approximations in chemical kinetics: one stage beyond the elimination of a fast variable
Classical approximations in chemical kinetics, the quasi-steady-state approximation (QSSA) and the partial-equilibrium approximation (PEA), are used to reduce rate equations for the concentrations and the extents of the reaction steps, respectively. We make precise two conditions on the rate constants necessary and sufficient to eliminate a wellchosen variable in the vicinity of a steady state. The first condition expresses that dynamics admits a small characteristic time associated with a fast variable. The second condition ensures that the fast variable is a concentration for QSSA and an extent for PEA. Both approximations exploit the zeroth order of a singular perturbation method.Eliminating a fast variable does not mean that it has reached a steady state. The fast evolution is considered over and the slow evolution of the eliminated variable is governed by the slow variables. The evolution of the slow variables occurs on a slow manifold in the space of the concentrations or the extents. In some cases the dynamics of the slow variables can be associated with a reduced chemical scheme. QSSA and PEA are applied to three chemical schemes associated with linear and nonlinear dynamics. We find that QSSA cannot be identified with the elimination of a reactive intermediate. The nonlinearities of the rate equations induce a more complex behavior.
Introduction<!>Evolution of the concentrations<!>Evolution of the extents of the reaction steps<!>Linear dynamics<!>Nonlinear dynamics<!>Comparison of QSSA application to linear and nonlinear dynamics<!>Parallel with the Born-Oppenheimer approximation<!>partial-equilibrium approximation (PEA)<!>Conclusion
<p>The laws of chemical kinetics lead to systems of differential equations for the concentrations of the reactive species. Depending on the fluxes of matter and energy, a reactive fluid maintained far from equilibrium may present qualitative changes of behaviors called bifurcations, such as the emergence of regular or chaotic temporal oscillations. 1 Following the dynamical systems theory, all systems of equations exhibiting a typical bifurcation can be reduced to a normal form involving a small number of variables and specific nonlinearities. 2 Knowing these nonlinearities is therefore sufficient to predict the behavior of a complex system. For example, this approach has been successfully applied to the modeling of biological rhythms such as glycolytic oscillations. 3 Even if systems with a large number of variables can be solved numerically, the approximations used to eliminate variables and determine reduced systems of equations remain crucial.</p><p>Two approximations based on timescale analysis are usually applied in chemical kinetics. 4 The quasi-steady-state approximation (QSSA) is used to eliminate a chemical species which evolves faster than the others. The partial-equilibrium approximation (PEA) is employed to remove a reaction step associated with an extent reaching a steady state faster than the other extents. 5 These approximations have been extensively studied over the last century. [6][7][8][9][10][11][12] Both methods are specific examples of the elimination of a fast variable, also called adiabatic elimination. [13][14][15] Nevertheless QSSA and PEA do not simply consist in eliminating a fast variable from a system of rate equations. The eliminated variable must be either a concentration for QSSA 16 or an extent for PEA. 5 Both methods are singular perturbation methods introducing a small parameter in the rate equations. [17][18][19][20][21] However QSSA and PEA are based on the zeroth order of the expansion, so that the small parameter does not explicitly appear in the reduction procedure. Consequently the validity of the expansion is not usually checked. In addition the higher-order expansion does not help, since it reintroduces the fast variable. Both QSSA and PEA exploit the existence of different timescales: 22,23 The trajectories in the space of the concentrations or extents are quickly attracted by a lower-dimensional invariant manifold in which the evolution continues more slowly. [24][25][26] In the two approximations the fast evolution is considered instantaneous and only the slow evolution on the invariant manifold is described.</p><p>In other words the fast variable, concentration or extent, is supposed to instantaneously adapt to the evolution of the slow variables. Quasi-steady state does not mean that the fast variable has reached a steady state but that it keeps evolving at a rate imposed by the slow variables. For this reason, the fast variable is sometimes called a slave variable. 27 Within the framework of the approximation the fast variable is not independent and can be ignored in the reduced dynamics.</p><p>In this article we present necessary and sufficient conditions for eliminating the concentration of a chemical species or the extent of a reaction step in the rate equations in the linear domain around a steady state. These conditions are valid in the entire space of the variables for first-order kinetics. The conditions apply to any reaction scheme and do not require the prior knowledge of the rate constant values but they can be checked when the rate constants are known.</p><p>The manuscript is organized as follows. In section 2 the rate laws associated with a general reaction scheme are written when choosing either concentrations or extents as dynamical variables. Section 3 is devoted to the quasi-steady-state approximation and section 4 to the partial-equilibrium approximation. Section 5 gives conclusion.</p><p>where ν ij and ν ij are the left and right stoichiometric coefficients associated with species C i and reaction j. The forward and backward rate constants of reaction j are denoted k j and k −j for j = 1, .., n r . The system can be described by two sets of variables, either the concentrations or the extents. In the following we present the basics of the resolution of linearized rate equations around a steady state using matrix methods. 28 The eigenvalues and eigenvectors of the Jacobian matrix are involved in the elimination conditions presented in section 3.</p><!><p>The system of ordinary differential equations governing the evolution of the n s concentrations C i is a coupled system of equations</p><p>If n c conservation laws can be written, only n = n s − n c concentrations are independent.</p><p>The system is written in matrix form for the n independent variables</p><p>where</p><p> is the concentration vector and f is a possibly nonlinear function of the variables. The solution of Eq. ( 3) is highly dependent on the nonlinearities of the function f . The steady state obeys the system of equations dC i dt = 0 for i = 1, ..., n. The small deviation from the steady-state vector C 0 is</p><p>Substituting for the concentration vector from Eq. ( 4) into Eq. ( 3) and linearizing around the steady state, we have</p><p>where the Jacobian matrix is</p><p>Equation ( 5) is a linear system of n equations approximately accounting for the dynamics around the steady state. In the case of a linear function f , the Jacobian matrix J is independent of the steady state and the resolution of Eq. ( 5) gives the exact evolution of the concentrations even far from the steady state. We consider the case where the Jacobian matrix J has n eigenvalues denoted λ i for i = 1, .., n. The eigendecomposition of the Jacobian matrix J leads to</p><p>where Λ is a diagonal matrix whose diagonal terms are the eigenvalues and P is the change-of-basis matrix. In the eigenbasis, the rate equations are uncoupled and can be written as</p><p>for i = 1, .., n. The vector of the variables in the eigenbasis</p><p>The solution of Eq. ( 8) is</p><p>for i = 1, .., n, where γ ini i is the initial condition of the variable γ i . At a given time, the state of the system is represented by a point in the n-dimension space of the concentrations. In the same space, the evolution of the system is a curve named trajectory. If the real parts Re(λ i ) of all eigenvalues λ i , i = 1, .., n, are negative the steady state C 0 is stable. A stable steady state is an attractor, i.e. for well-chosen initial conditions the concentrations eventually reach the stable steady concentrations. The trajectory converges towards the steady state. If Re(λ i ) > 0 for at least one eigenvalue λ i , the steady state is unstable.</p><p>The time τ i = 1/ | Re(λ i ) | associated with the variable γ i characterizes the evolution along the i th eigendirection, a large value of τ i being associated with a slow variable and a small value with a fast variable. Equation ( 9) written for t = 0 gives the relation between the initial conditions γ ini i for the variables in the eigenbasis and the initial conditions C ini i for the concentrations, leading to</p><p>The evolution of the concentration C i is then obtained using Eqs. (4,9-11)</p><p>where P ij and P −1 ij are elements of the matrices P and P −1 , respectively. The concentrations are linear combination of the variables in the eigenbasis. Equation (12) shows how each variable γ i contributes to the evolution of the concentrations.</p><!><p>The law governing the evolution of the extent ξ j of step j given in Eq. ( 1) is</p><p>Substituting for the right-hand side from Eq. ( 13) into Eq. ( 2) and integrating on the interval [t, ∞[, we obtain a linear transformation between the two sets of variables, concentrations and extents. We choose that each extent ξ j vanishes at the steady state, leading to</p><p>Substituting C i from Eq. ( 14) into Eq. ( 13), we find the ordinary differential equations governing the evolution of the extents</p><p>The method of resolution used for concentrations can be straightforwardly implemented for the extents. The concentrations and the extents exhibit linear relationships. Consequently the numbers of independent concentrations and independent extents are equal.</p><p>In addition, the concentrations and the extents have the same characteristic times.</p><p>3 Quasi-steady-state approximation (QSSA)</p><p>The quasi-steady-state approximation (QSSA) is a well-known method used to eliminate a fast concentration. [29][30][31][32] For example, when one of the chemical species C i is much more reactive than the other chemical species, the fast variation of C i is over at the timescale at which the other chemical species are still evolving. The approximation consists in writing dC i dt = 0 for a single species C i . According to Eq. ( 2) the concentration C i can then be written as a function of the other concentrations, which reduces by one the number of independent variables in the dynamics. QSSA is used to eliminate the fast concentration C i in the rate equations and we make precise the validity domain of the approximation in the following.</p><p>Fast and slow variables are defined in the eigenbasis of the linearized system, consequently in the vicinity of the steady state. The elimination of C i can be performed if two conditions are fulfilled. First a variable γ k in the eigenbasis must be faster than the other variables, i.e.</p><p>Second the concentration C i must vary as the fast variable γ k during the time interval [0, τ k ]. Using the normalized change-of-basis matrix P and Eq. ( 9) we write the condition as</p><p>According to Eq. ( 17), the slow eigendirections are close to the hyperplan defined by</p><p>According to QSSA, the relaxation with the short characteristic time τ k along the γ k -axis can be considered instantaneous. The conditions given in Eqs. (16,17) involve the rate constants. They are necessary but not sufficient to ensure the validity of the elimination of the variable C i in the entire space of the concentrations since they result from a linearization of dynamics around the steady state.</p><p>Even if these conditions are met, the approximation is a priori not valid at the beginning of the evolution of the concentrations. It is valid only after the fast variable γ k has relaxed. Between 0 and τ k the trajectories in the space of the concentrations are close to the k th fast eigendirection which is close to the C i -axis. After the relaxation of the fast variable γ k the trajectories evolve on an (n − 1)-dimensional space called the slow manifold and defined by the single equation</p><p>At the fixed point corresponding to the steady state, the slow manifold is tangent to the vector space spanned by the n − 1 eigendirections associated with the slow variables. On the slow manifold, the concentration C i slowly evolves and adapts to the evolution of the slow concentrations.</p><p>After the elimination of the concentration C i , the system of differential equations governing the evolution of the other concentrations can be written as a system of dimension n − 1. If this system of differential equations involves polynomials compatible with the rate laws of reaction steps, the evolution of the remaining concentrations can be associated with a reduced set of reaction steps. The existence of a reduced reaction scheme is not always guaranteed as illustrated by Michaelis-Menten kinetics which leads to rational functions of the concentrations in the rate equations. 32 Approximate initial conditions adapted to the elimination of a concentration can be defined as the intersection of the slow manifold and the fast direction passing through the initial conditions before elimination.</p><p>If the remaining characteristic times have different orders of magnitude, an analogous procedure of elimination may be performed to further reduce the number of variables.</p><p>For illustration we consider two reaction schemes differing by one reaction step and associated with either linear or nonlinear rate equations. Using a chemostat, i.e. a reservoir of some species, we limit the number of variables without the need to introduce laws of conservation.</p><!><p>We consider the first-order reaction scheme</p><p>where A is a reservoir of A species also called chemostat, implying that the concentration A is constant. Only the concentrations X and Y associated with chemical species X and Y are independent variables obeying the linear system of ordinary differential equations</p><p>where k 1 , k −1 , k 2 , and k −2 are rate constants. The steady state is</p><p>. Following the method detailed in section 2 for the Jacobian matrix</p><p>we obtain the analytical expressions of the concentrations, which are valid in the entire concentration space (X, Y ) due to the linear nature of the equations given in Eqs. (21,22). The eigenvalues of J are</p><p>and the change-of-basis matrix is</p><p>The eigenvalues are real negative numbers regardless of the values of the rate constants.</p><p>The steady state is thus always stable. Interestingly the rate constant k −2 does not appear in the Jacobian matrix and consequently in the eigenvalues and the change-of-basis matrix. The conditions given in Eqs. (16,17) applied to the elimination of species Y read</p><p>where</p><p>are deduced from Eqs. (24,25) and P 21 and P 22 are given in Eq. (26). Figure 1a displays the domain in the rate constant space in which the conditions given in Eqs. (27,28) are met. The two conditions are numerically found to be equivalent to</p><p>The usual condition for the application of QSSA is recovered: Y can be considered a very reactive intermediate, slowly formed and rapidly consumed. 5 The variable γ 2 associated with the smallest characteristic time τ 2 is the fast variable and evolves like Y . The variable γ 1 associated with the largest characteristic time τ 1 is the slow variable.</p><p>According to Eqs. (18,22) the equation of the slow manifold is dY /dt = 0, i.e.</p><p>QSSA does not mean that the fast variable Y has reached a steady state but that its evolution depends on X and follows the slow manifold. Substituting for the concentration Y from Eq. ( 29) into Eq. ( 21), we obtain the rate equation governing the evolution of the concentration of species X after elimination of Y</p><p>This equation can be associated with the reduced reaction scheme</p><p>Using an expansion in power of k 1 /(k −1 + k 2 ), we check that the characteristic time τ 1 associated with the slow variable converges toward the characteristic time τ = (k</p><p>Approximate initial conditions associated with the reduced reaction scheme can be deduced from the initial conditions X ini and Y ini associated with linear dynamics. Using the second column of the change-of-basis matrix given in Eq. ( 26), we find that the straight line parallel to the fast eigendirection and passing through the point (X ini , Y ini ) is given by</p><p>Solving the system given in Eqs. (29,32), we obtain the initial condition associated with the reduced reaction scheme</p><p>Before studying the behavior of the trajectories in the space (X, Y ) in subsection 3.3, we illustrate QSSA in the case of nonlinear dynamics.</p><!><p>We consider the following reaction scheme</p><p>where A is a chemostat.</p><p>Following the laws of chemical kinetics, the rate equations for the concentrations X and Y are nonlinear and given by</p><p>). Linearizing Eqs. (36,37) around the steady state and following the same method as previously for the Jacobian matrix</p><p>with κ = k 1 k −2 A, we find the eigenvalues</p><p>The change-of-basis matrix is</p><p>The steady state is stable since the eigenvalues are always real and negative regardless of the values of the rate constants.</p><p>The conditions for the elimination of Y are given in Eqs. (27,28) where λ 1 and λ 2 are given in Eqs. (39,40) and P 21 and P 22 , in Eq. ( 41). The variable Y can be eliminated if the parameters meet the conditions given in Eqs. (27,28) as shown in Fig. 1b. The variables γ 2 and γ 1 are the fast and slow variables, respectively.</p><p>According to Eqs. (18,37), the slow manifold obeys</p><p>Substituting for the concentration Y from Eq. (42) into Eq. ( 36), we have</p><p>The associated reduced reaction scheme is</p><p>Linearizing Eq. ( 43) around the steady state X 0 = k −1 k −2 A k 1 k 2 , we find that the variable X obeying the reduced reaction scheme evolves with the characteristic time</p><p>Approximate initial conditions associated with the reduced reaction scheme given in Eq. ( 44) can be deduced from the initial conditions X ini and Y ini associated with the initial reaction scheme. Using the second column of the change-of-basis matrix given in Eq. ( 41), we obtain the equation of a straight line parallel to the fast eigendirection and</p><p>Solving the system given in Eqs. (42,46), we obtain the initial condition X ini of the reduced reaction scheme</p><p>where</p><p>The linear or nonlinear nature of the reaction schemes induces differences in the elimination of the variable Y that are analyzed in the next subsection.</p><!><p>As shown in Fig. 1, the linear and nonlinear dynamics lead to very different parameter domains in which the elimination of Y is valid. Figure 2 displays the eigendirections given in Eq. ( 23) (Eq. (38), resp.) and the slow manifold given in Eq. ( 29) (Eq. ( 42), resp.) for linear (nonlinear, resp.) dynamics.</p><p>Trajectories in the space (X, Y ) deduced from the numerical integration of Eqs. (19,20)   for linear dynamics and Eqs. (20,34) for nonlinear dynamics are also shown for rate constant values in the validity domain of Y elimination. For linear dynamics, the slow manifold is a straight line identical to the slow eigendirection. In this case, the slow manifold is a vector space. For nonlinear dynamics, the slow manifold is curved but, as expected, is tangent to the slow eigendirection at the steady state. For both dynamics, the fast eigendirections are quasi-parallel to the Y -axis. The trajectories are first parallel to the fast eigendirection, then close to the slow manifold, and eventually reach the steady state in the slow eigendirection. The elimination of the variable Y is valid as soon as the trajectories have reached the vicinity of the slow manifold. The rate constants are expressed in arbitrary units.</p><p>The variation of X at large timescales (t τ 2 ) is correctly described by the reduced reaction schemes with the initial conditions given in Eqs. (33,47).</p><p>Figure 3a displays the variation of the absolute characteristic times deduced from Eqs. (39,40) with the rate constant k 2 for nonlinear dynamics. The characteristic time τ given in Eq. ( 45) associated with the reduced reaction scheme is also displayed. In agreement with Fig. 1b, the variation of τ is similar to the variation of τ 1 for sufficiently large k 2 .</p><p>Figure 3b shows the variation of the change-of-basis matrix elements given in Eq. ( 41</p><!><p>As shown in Fig. 3a for the chemical scheme given in Eqs. (34,35), the characteristic times τ 1 and τ 2 display a typical avoided crossing in the sense of quantum chemistry. A parallel between QSSA and the Born-Oppenheimer approximation can be drawn. Both approximations are adiabatic eliminations. 6 The Born-Oppenheimer approximation consists in neglecting the fast movements of the electrons and considering only their response to the slower movement of the nuclei. For degenerate energy states, i.e. identical eigenvalues of the Hamiltonian, the Born-Oppenheimer approximation is not valid. Introducing a perturbation splits the degenerate energy states and leads to an avoided crossing. Similarly, the introduction of a fast variable in a reaction scheme can be considered a perturbation leading to an avoided crossing for the characteristic times. The trajectories are deduced from the numerical integration of Eqs. (21,22) for linear dynamics and Eqs. (36,37) for nonlinear dynamics. The initial conditions (X ini , Y ini ) are (0, 0.</p><!><p>The partial-equilibrium approximation (PEA) offers another way to eliminate a fast variable, systemically leading to a reduced reaction scheme. 25,33 When a reaction step is much faster than the others, the fast variation of its extent ξ i can be eliminated at the timescale of the evolution of the other extents. The approximation consists in writing dξ i dt = 0. This approximation is similar to QSSA, with the only difference that the eliminated variable is an extent and not a concentration.</p><p>The same comments as for the elimination of a concentration can be made for the elimination of an extent. In summary, the extent ξ i can be eliminated in the vicinity of the steady state if the following conditions are met. First a variable γ k in the eigenbasis of the Jacobian matrix associated with the extents must be faster than the other variables, leading to the condition given in Eq. ( 16). Second the variables γ k and ξ i must evolve similarly on the time interval [0, τ k ], leading to the condition given in Eq. ( 17). The elimination of the extent ξ i is valid when the variable γ k associated with the smallest characteristic time τ k has relaxed. In the space of the extents, the trajectories are first parallel to the fast k th eigendirection, i.e. the ξ i -axis, and then follow the slow manifold defined by</p><p>The slow manifold is tangent to the hyperplane defined by the slow eigendirections.</p><p>Contrary to QSSA, PEA always leads to a reduced chemical scheme straightforwardly obtained by removing the eliminated reaction step from the initial scheme.</p><p>We illustrate PEA using the following reaction scheme</p><p>where A is a chemostat. The rate laws for the concentrations are given by</p><p>The steady state is</p><p>The extents ξ 1 and ξ 2 are associated with the steps given in Eq. (49) and Eq. (50), respectively. Using Eq. ( 14) we obtain the relationships between the extents and the concentrations</p><p>Following Eq. ( 13), we find</p><p>Substituting for X and Y from Eqs. (54,55) into Eqs. (56,57) we obtain the rate laws for the extents</p><p>where X 0 and Y 0 are given in Eq. ( 53). We follow the resolution method detailed in section 2. The system given in Eqs. (58,59) is linearized leading to the Jacobian matrix</p><p>The eigenvalues of J are given by</p><p>The change-of-basis matrix is</p><p>The eigenvalues of J are negative for all values of the rate constants. Consequently the steady state (X 0 , Y 0 ) is stable. The conditions given in Eqs. (16,17) applied to the elimination of the extent ξ 1 are written as  According to Eqs. (48,56,58), the slow manifold is defined by</p><p>in the space of the extents and by</p><p>in the space of the concentrations.</p><p>The elimination of the fast extent ξ 1 leads to the reduced reaction scheme</p><p>The mere knowledge of the reduced chemical scheme is not sufficient to correctly describe the evolution of the concentrations X and Y after the elimination of the extent ξ 1 . The concentrations X and Y are not independent but obey Eq. (67). Equations (67,68) are necessary to write the rate laws associated with the reduced reaction scheme.</p><p>Trajectories in the spaces (ξ 1 , ξ 2 ) of the extents and (X, Y ) of the concentrations are shown in Fig. 5 with rate constant values for which the elimination of the extent ξ 1 is valid. The eigendirections are obtained using Eq. ( 63) and the slow manifold using Eq.</p><p>(66) in the space (ξ 1 , ξ 2 ) and Eq. (67) in the space (X, Y ). The trajectories are obtained by numerical integration of Eqs. (58,59) in the space (ξ 1 , ξ 2 ) and Eqs. (51,52) in the space (X, Y ). As expected the fast eigendirection is parallel to the ξ 1 -axis in the space (ξ 1 , ξ 2 ). Both eigendirections are not parallel to any axis in the space (X, Y ) since the evolution of X and Y depends on both characteristic times τ 1 and τ 2 . In both spaces, the slow manifold is tangent to the slow eigendirection. In addition, the trajectories are first parallel to the fast eigendirection and evolve toward the steady state along the slow manifold. The reduced reaction scheme is valid once the trajectory have reached the vicinity of the slow manifold.</p><!><p>Chemical kinetics is usually described by the evolution of two sets of variables, either the concentrations of the chemical species or the extents of the reaction steps. The concentrations and the extents are linked by linear relationships. Laws of conservation are used to reduce dynamics to the evolution of independent variables. The quasi-steady-state and partial-equilibrium approximations require more hypotheses than the mere elimination of a fast variable in rate equations since the eliminated variable must have a physical meaning. It must be a chemical species in the case of QSSA and an extent in the case of PEA.</p><p>We have proposed a general framework to make precise the minimum number of conditions on the rate constants for which either approximation is valid in a domain of the space of the variables. The knowledge of the values of all rate constants is not necessary to apply one or the other approximation. We exploit the eigendecomposition into eigenvalues and eigenvectors of the Jacobian matrix deduced from the linearization of the rate equations around a steady state. The absolute real part of the eigenvalues correspond to the inverse of the characteristic times. Eliminating a fast variable is possible as soon as one characteristic time is one order of magnitude smaller than the others. We write a second condition on the eigenvectors imposing that the eliminated variable is the concentration C i for QSSA or the extent ξ i for PEA. The two conditions are necessary and sufficient in the linear domain around the steady state, i.e. in a sufficiently small neighborhood of the steady state in the space of the concentrations or extents.</p><p>Both approximations are zeroth-order perturbation methods consisting in considering that the relaxation of the fast variable is instantaneous leading to dC i /dt = 0 or dξ i /dt = 0. Either of these equations defines the slow manifold, i.e. a relationship between the fast variable and the slow variables. After the elimination of the fast variable, dynamics occurs on the slow manifold. Substituting for the fast variable from the slow manifold equation into the rate equations of the slow variables leads to a reduced dynamics. For some chemical schemes QSSA and PEA may lead to the same reduced rate equations when the conditions of validity of the two approximations overlap. Far from being stationary, the evolution of the eliminated variable occurs at the timescale on which the reduced rate equations focus. To point out this essential feature, we suggest to use the name quasi-steady-state instead of steady-state approximation and partial-equilibrium instead of equilibrium approximation. It is to be noted that QSSA does not always provide a reduced mechanism since the elimination of the fast variable may lead to non polynomial rate equations for the slow variables as illustrated by the familiar Michaelis-Menten equations. In the case of PEA, the reduced mechanism is straightforwardly obtained but may not be sufficient to account for dynamics. The equation of the slow manifold may be necessary.</p>