THE EFFECT OF LUMPING AND EXPANDING ON KINETIC DIFFERENTIAL EQUATIONS

János Tóth Department of Computer Science University of Agricultural Science H-2103 Gödöllő, Hungary

Genyuan Li, Herschel Rabitz†, and Department of Chemistry Princeton University Princeton, New Jersey 08540

Alison S. Tomlin Department of Fuel and Energy The University of Leeds Leeds, LS2 9JT, Great Britain

Abbreviated title. Lumping kinetic differential equations

Keywords. lumping, chemical kinetics, kinetic differential equations, aggregation, global decomposition, factorization, qualitative properties, eigenvalues, Lyapunov functions

AMS(MOS) subject classifications. 80A30, 34Cxx, 58Fxx † Author to whom correspondence should be addressed.

Abstract. Let us consider the differential equation: ẏ(t) = f(y(t)) with an f from R^{N} to R^{N} and suppose that there exists a transformation h from R^{N} to RN̂ (N̂ ≤ N) such that ŷ := h ∘ y obeys a differential equation ˙ ŷ(t) = f̂(ŷ(t)) with some function f̂, then the first equation is said to be lumpable to the second by h. Here mainly the case is investigated when the original differential equation has been induced by a complex chemical reaction. We provided a series of necessary and sufficient conditions for the existence of such functions h and f̂, some of them are formulated in terms of h and f only. Beyond these conditions our main concern here is how lumping changes those properties of the solutions which are either interesting from the point of view of the qualitative theory of differential equations, or from the point of view of formal reaction kinetics.

Results show that each eigenvalue of the Jacobian of the nonlinear lumped system at an equilibrium is an eigenvalue of the original system at the corresponding equilibrium. (Invariant sets, equilibria and periodic solutions are lumped into invariant sets, equilibria and periodic solutions, respectively.) Under certain conditions a Lyapunov function of the lumped system can be used to create a Lyapunov function for the original one (to test relative stability) and vice versa—both around equilibria and far from equilibria.

These general statements do not necessarily imply close qualitative resemblance of lumped and original solutions but provide criteria by which lumping schemes may be selected for this purpose. The precise meaning of the conditions in the general statements is illustrated by examples taken from formal reaction kinetics.

0. Introduction. Mathematical models of chemical reaction kinetics are used in many fields such as in basic research in physical chemistry and biochemistry; in the design and control of chemical reactors, in combustion, in atmospheric chemistry, in pharmacokinetics etc. The use of these models is often restricted because of the large number of dependent variables, which usually are the species’ concentrations. There are a lot of solutions to circumvent this problem. The technique we address here is lumping: the reduction of the number of variables by grouping them via a linear or nonlinear function. This approach was initiated by Wei and Kuo (1969) who gave a lumping analysis of unimolecular reaction systems where the lumped species were certain linear combinations of the original ones. Their work was extended by the authors to exact linear (Li and Rabitz, 1989, 1991a) and exact nonlinear (Li et al., 1994a); approximate linear and nonlinear (Li and Rabitz, 1990, 1991a,b; Li et al., 1993a,b; 1994a,b,c) lumping of arbitrary nonlinear differential equations (see also the review by Li and Rabitz, 1991d), and these methods have also been applied to practical problems arising in petroleum cracking (Li and Rabitz, 1991b), catalytic reforming (Li, Rabitz, 1993a), biochemistry (Li et al., 1993b, 1994b), combustion (Li et al., 1994c), etc.

Here we wish to investigate if lumping (or expanding, i.e. returning to the original system from the lumped one) has any effect on the transient and long range properties of the solutions. The most natural type of question is: if we learn something about a solution of the lumped model, can it also be stated about a certain solution (some solutions) of the original one? However, in order to be able to answer such questions, one has to investigate the corresponding direct problem, namely, how the properties are changed under lumping. Finally, let us remark that it is not a Law of Nature that a lumped (or expanded) system should have anything in common with the original one. We would therefore like to investigate which lumping (or expanding) schemes conserve which properties.

Let us also mention that similar global reduction methods (sometimes under the name aggregation) are needed and have been investigated in many areas such as control theory, econometrics, biology and ecology (Luenberger, 1964; Ijiri, 1971; Łos and Łos, 1974; Iwasa et al., 1987, 1989; Luckyanov et al., 1983).

Standard textbooks on differential equations have nothing to say about lumping or aggregation. Keener’s textbook treated singular perturbation theory of the possible approaches (Keener, 1988). Korolyuk and Turbin (1993) in Sec. 7.1 show that reducing the number of variables by singular perturbation can be investigated within the theory they developed for singularly perturbed semigroups. Wiggins (1990, p.

193) mentioned dimension reduction as a method for simplifying dynamical systems, but he only suggested center manifold theory as a rigorous mathematical theory for this purpose. Whereas that approach admittedly finds its roots in the intrinsic nature of and in the classification of solutions, lumping is a framework which came into being in order to adequately describe some practical systems. There are important applications of lumped species in measurement techniques where the output may be some linear or nonlinear function of the original variables. As it will be seen below, lumping can also be considered as a generalization of C^{2}-conjugacy or C^{2}-equivalence of vector fields from the viewpoint of the theory of ordinary differential equations; a result probably most fundamental from the mathematical point of view. It is also true that lumping is in close relationship with global decomposition of nonlinear control theory (Isidori, 1985), or to factorization (BenArtzi and Helton, 1988). The only mathematical textbook to treat our problem in a similar way, although in less detail, is that by Godbillon (1969, Ch. 5).

The structure of our paper is as follows. Section 1 recapitulates the most important definitions and statements on lumping which will be used in the following sections. In order to make the paper self-contained we included full proofs of our earlier results on necessary and sufficient conditions. We also list some definitions of the qualitative theory of ordinary differential equations and of formal reaction kinetics in the form we shall use them below. Section 2 deals with the effect of lumping on general regularity properties of the right hand sides. Invariant sets, including equilibria, periodic solutions and the first orthant of the phase space are treated in Section 3. Section 4 contains the major results that some of the eigenvalues and Lyapunov functions at equilibria may sometimes be preserved under lumping. Consequences on stability and the selection of possible lumping schemes are treated in detail. Next, in Section 5, time dependent transient solutions including monotone, oligo-oscillatory, periodic and chaotic solutions are treated, or at least mentioned.

Before the Summary and discussion we present a collection of unsolved problems.

One remark on notations might be useful for the reader: y and ŷ normally denote vector valued functions of time, whereas the same letters with any indices, star or tilde denote vectors of the corresponding space.

1. Basic notions.

1.1. Definitions and theorems from calculus and from the qualitative theory of differential equations. Let N be a positive integer (N ∈ N, as the set of positive integers will, as usual, be denoted by N; later we shall also use N_{0} to denote the set of nonnegative integers.)

and consider a differentiable function f from R^{N} to R^{N} with f(0) = 0. Consider the differential equation:

(1.1)  ẏ(t) = f(y(t)),

and suppose, if not said otherwise (as in Section 5), that the domain of all the solutions can be extended to all t ∈ R. (To put it another way: there is no finite escape time to infinity, or the maximal interval of existence of the solutions is R.) This equation is said to induce the dynamical system ψ : R × R^{N} ⟶ R^{N} where ψ(t, y_{0}) is the solution of (1.1) with the initial condition y(0) = y_{0} at time t ∈ R.

We also use the notation ψ_{t} := ψ(t, .) (t ∈ R, R_{0}^{+}, R_{0}^{−}, respectively) to denote the flow (or, the positive and negative semiflow) induced by (1.1).

DEFINITION 1.1. A set S ⊂ R^{N} is said to be positively invariant (negatively invariant, invariant, respectively), if ψ_{t}(S) ⊂ S for all t ∈ R_{0}^{+} (R_{0}^{−}, R, respectively).

DEFINITION 1.2. The point y^{∗} ∈ R^{N} for which

(1.2)  f(y^{∗}) = 0

holds is an equilibrium. It is a sink if all the eigenvalues of f′(y^{∗}) have negative real part; a source if all the eigenvalues of f′(y^{∗}) have positive real part; finally, it is a hyperbolic equilibrium point if none of the eigenvalues have zero real part.

According to the assumption above, 0 is one of the equilibria of (1.1).

The next statement can be found in e.g. Hahn (1963), or in Rouchet et al. (1977, Th. 7.2).

THEOREM 1.1. An equilibrium y^{∗} of (1.1) is stable (asymptotically stable) if and only if there exists a neighborhood U(y^{∗}) of y^{∗} and a function

V ∈ C^{1}(U(y^{∗})\{y^{∗}}, R^{+})

such that

(1.3)  V̇_{(1.1)}(ỹ) := (V′f)(ỹ) ≤ 0 (< 0, respectively)(ỹ ∈ U(y^{∗})\{y^{∗}})

holds.

Such a function V is usually referred to as a Lyapunov function.

A trivial but often used statement follows. Let K, N ∈ N.

STATEMENT 1.1. Suppose that

k ∘ y = l ∘ y

holds for the functions k, l : R^{N} ⟶ R^{K} with some functions y : R ⟶ R^{N} for which ∪_{y}Range(y) = R^{N} is true. Then k = l is also true.

The statement is usually used for all the solutions y of a differential equation with a right hand side defined on the whole R^{N}.

1.2. Fundamentals of lumping. Now let us suppose that f is twice continuously differentiable and consider the differential equation (1.1).

DEFINITION 1.3. Suppose there exists a positive integer N̂ ≤ N and a twice continuously differentiable transformation h from R^{N} to RN̂ with the properties

(i) h(0) = 0,

(ii) there exists u : R^{N} ⟶ R^{N−N̂} such that ( \begin{matrix}h′(ỹ) \\ u′(ỹ) \end{matrix} ) is nonsingular for all ỹ ∈ R^{N}, and

(iii) lim_{|_{y}|→∞} |[ \begin{matrix}h(ỹ) \\ u(ỹ) \end{matrix} ]| = ∞.

If for all solutions y to (1.1)

(1.4)  ŷ := h ∘ y

obeys a differential equation

(1.5)  ˙  ŷ(t) = f̂(ŷ(t))

with some continuously differentiable function f̂ : RN̂ ⟶ RN̂, then (1.1) is said to be exactly lumpable to (1.5) by h. The pair of functions h and f̂ is sometimes referred to as a lumping scheme.

Suppose \overline{h} is a mapping from RN̂ to R^{N} that satisfies

(1.6)  h ∘ \overline{h} = IN̂,

where IN̂ is the identity mapping of RN̂. Such a mapping \overline{h} always exists if h fulfills the properties (ii) and (iii) above, and it will be referred to as a generalized inverse of h. (On the existence of a global inverse see e.g. Rădulescu and Rădulescu, 1980, and the references therein.) From now on properties (i)–(iii) will be understood by the name nondegenerate, as well as twice continuous differentiability, and nondegeneracy of h shall usually be assumed from now on. — The case N = N̂ is also allowed so as not to exclude coordinate transformations.

DEFINITION 1.4. Suppose we are given (1.5) and there exist h and f such that (1.1) and (1.4) hold. Then (1.5) is said to be exactly expandable to (1.1).

This definition is only provided to shorten certain statements since all differential equations can be expanded in a trivial way. Let us start from

˙ŷ = f̂ ∘ ŷ

and let y be defined as

y := ( ŷ z )

with an arbitrary scalar-vector function z for which an autonomous equation, say, ż = g ∘ z holds. Then the projection h := pr_{1} or

h(ỹ) := ỹ (ỹ ∈ R^{N})

is nondegenerate and lumps the equation

ẏ = f ∘ y

(where f(˜ŷ, z̃) := ( \begin{matrix}f̂(˜ŷ) \\ g(z̃) \end{matrix} ) , or f := ( \begin{matrix}f̂ ∘ pr_{1} \\ g ∘ pr_{2} \end{matrix} )) into the starting one.

Conceptually, expansion is the same as inverse lumping by Schlijper and Drohm (1988).

For completeness, we briefly restate the conditions given for the existence of a lumping scheme.

THEOREM 1.2. Eqn. (1.1) is exactly lumpable to (1.5) by the nondegenerate function h if and only if

(1.7)  h′f = f̂ ∘ h

holds. The representation

(1.8)  f̂ = (h′f) ∘ \overline{h}

is also valid, where \overline{h} is a generalized inverse of h.

Proof. Taking the (time) derivative of h ∘ y and using (1.1) one obtains

˙ŷ = (h ∘ y)^{·} = h′ ∘ y · ẏ = h′ ∘ y · f ∘ y = (h′f) ∘ y,

and the final expression should equal f̂ ∘ h ∘ y. Therefore Statement 1.1 implies necessity.

To prove sufficiency, let us suppose that (1.7) holds and calculate again the derivative of h ∘ y to obtain

˙ŷ = (h ∘ y)^{·} = h′ ∘ y · ẏ = h′ ∘ y · f ∘ y = (h′f) ∘ y = f̂ ∘ h ∘ y = f̂ ∘ ŷ,

showing that (1.1) is exactly lumpable to (1.5) by h.

Composing both sides of (1.7) with \overline{h} from the right provides the representation (1.8).

Notice that only the continuity of f has been used in the proof of the statement above.

THEOREM 1.3. Let \overline{h} be a generalized inverse of the nondegenerate function h : R^{N} ⟶ RN̂. Then (1.1) is exactly lumpable to (1.5) by h if and only if

(1.9)  h′f = (h′f) ∘ \overline{h} ∘ h

holds.

Proof. Composing both sides of (1.8) with h from the right and comparing the result to (1.7) provides (1.9).

To prove sufficiency, let f̂ := h′f ∘ \overline{h} and use (1.9) to obtain (1.7):

f̂ ∘ h = (h′f ∘ \overline{h}) ∘ h = h′f,

which is known to be a sufficient condition to exact lumpability.

Notice that (1.9) does not contain f̂, and (1.7) does not contain \overline{h}.

THEOREM 1.4. If (1.1) is exactly lumpable to (1.5) by h, then the right hand side f̂ of the lumped equation (1.5) does not depend on the specific choice of the generalized inverse of h.

Proof. Let \overline{h} and h̃ be two generalized inverses of h. Then, by (1.9) and (1.6) we have

(h′f) ∘ \overline{h} = (h′f ∘ h̃ ∘ h) ∘ \overline{h} = (h′f ∘ h̃) ∘ (h ∘ \overline{h}) = h′f ∘ h̃

showing through (1.8) the independence of the lumped right hand side on the choice of the generalized inverse of h.

One can also obtain higher order necessary and sufficient conditions.

THEOREM 1.5. Eqn. (1.1) is exactly lumpable to (1.5) by h if and only if

(1.10)  (h′f)′ = (h′f)′ ∘ \overline{h} ∘ h · \overline{h}′ ∘ h · h′.

Proof. Eqn. (1.10) is obtained from (1.9) by simply taking the derivatives of both sides.

Conversely, integrating (1.10) and taking into consideration the assumptions that f(0) = 0, h(0) = 0 and \overline{h}(0) = 0 gives (1.9).

THEOREM 1.6. Eqn. (1.1) is exactly lumpable to (1.5) by the nondegenerate function h : R^{N} ⟶ RN̂ if and only if there exists a matrix valued function X : R^{N} ⟶ R^{N×(N−N̂)} with rank N − N̂ everywhere for which

(1.11)  h′(ỹ)X(ỹ) = 0 (∈ R^{N̂×(N−N̂)})

and

(1.12)  (h′f)′(ỹ)X(ỹ) = 0 (∈ R^{N̂×(N−N̂)})

simultaneously hold for all ỹ ∈ R^{N}.

Proof. A) The first equation follows from the fact that h is nondegenerate, thus its derivative is of the full rank, therefore the number of independent vectors orthogonal to its rows is exactly N − N̂.

Upon multiplying both sides of (1.10) by X(ỹ) and using (1.11) yields (1.12). B) Suppose h is a function for which there exists an X such that (1.11) and

(1.12) hold. Let us introduce the function u ∈ C^{1}(R^{N}, R^{N−N̂}) in such a way that the transformation T defined by

(1.13)  T := ( \begin{matrix}h \\ u \end{matrix} )

is invertible and the Jacobian of T^{−1} is

(1.14)  T^{−1}′ = ( Y X) ∘ T^{−1}

with a certain function Y and the given function X taking on values in R^{N×N̂} and R^{N×(N−N̂)}, respectively.

Let y be a solution to (1.1) and let us consider the differential equation for z := T ∘ y:

(1.15)

ż = ( \begin{matrix}h′ ∘ y \\ u′ ∘ y \end{matrix} ) f ∘ y

= ( \begin{matrix}h′ ∘ T^{−1} ∘ z \\ u′ ∘ T^{−1} ∘ z \end{matrix} ) f ∘ T^{−1} ∘ z.

Upon denoting the right hand side of the above equation by g—i.e.

(1.16)  g = ( \begin{matrix}h′ ∘ T^{−1} \\ u′ ∘ T^{−1} \end{matrix} ) f ∘ T^{−1} = (( \begin{matrix}h′ \\ u′ \end{matrix} ) f) ∘ T^{−1}

—its Jacobian can be calculated as:

(1.17)

g′ =(( \begin{matrix}h′ \\ u′ \end{matrix} ) f)′ ∘ T^{−1} · T^{−1}′
= ( \begin{matrix}(h′f)′ \\ (u′f)′ \end{matrix} ) ∘ T^{−1} · ( Y X) ∘ T^{−1}
= ( \begin{matrix}(h′f)′Y & (h′f)′X \\ (u′f)′Y & (u′f)′X \end{matrix} ) ∘ T^{−1}.

Condition (1.12) implies

g′ = ( \begin{matrix}(h′f)′Y & 0 \\ (u′f)′Y & (u′f)′X \end{matrix} ) ∘ T^{−1},

i.e. the first N̂ coordinate functions of g do not depend on the last N −N̂ variables. Hence, the first N̂ equations in (1.15) compose a differential equation obtained from (1.1) by exact lumping.

Finally a necessary condition is given which might help to construct the lumping functions explicitly.

THEOREM 1.7. Eqn. (1.1) is exactly lumpable to (1.5) by the nondegenerate function h : R^{N} ⟶ RN̂ only if the manifold (being a closed set)

(1.18)  S := {ỹ ∈ R^{N}; h(ỹ) = 0}

is invariant under the flow induced by (1.1).

Proof. Suppose (1.1) is exactly lumpable to (1.5) by h, then Eq. (1.9) implies for all ỹ ∈ S that

(h′f)(ỹ) = (h′f)(\overline{h}(0)) = const.

However according to our starting assumptions 0 ∈ S, therefore this constant is zero:

(h′f)(\overline{h}(0)) = h′(0)f(0) = 0.

Now let us calculate the derivative of h along a solution y starting from the point y(0) ∈ S:

(h ∘ y)^{·} = (h′f) ∘ y = (h′f ∘ \overline{h} ∘ h) ∘ y.

Putting t = 0 into the above equation yields

(1.19)  (h ∘ y)^{·}(0) = 0.

Consider a trajectory of a solution to (1.1) with the initial condition in the manifold S. Suppose the trajectory leaves the boundary of S. Then, there would be a point in S, where (h ∘ y)^{·} would be nonzero contradicting to (1.19).

A possible strategy to create lumping functions would be to find all the functions h for which the manifold S defined above is invariant and then investigate each of them to see if they fulfill one of the sufficient conditions given above.

Now we formulate and prove a generalization of a statement on C^{2}-equivalent vector fields, see e.g. Perko (1991; p. 172), which might be considered as the most fundamental result from the mathematical point of view (even if the proof is simple).

THEOREM 1.8. The flow ψ̂_{t} (t ∈ R) induced by (1.5) is the function

(1.20)  h ∘ ψ_{t} ∘ \overline{h}.

Furthermore,

(1.21)  ψ̂_{t} ∘ h = h ∘ ψ_{t}

also holds.

Proof. Let ŷ_{0} := h(y_{0}) and let us show that t ↦ ψ̂_{t}(ŷ_{0}) is the solution of the differential equation (1.5) with the initial condition ŷ(0) = h(y_{0}). Let us start with the initial condition.

ψ̂_{0}(ŷ_{0}) = h(ψ_{0}(\overline{h}(ŷ_{0})) = h(\overline{h}(ŷ_{0})) = ŷ_{0}.

Now let us calculate the time derivative of t ↦ ψ̂_{t}(ŷ_{0}):

(1.22)

(ψ̂_{t})^{·}(ŷ_{0}) = h′(ψ_{t}(\overline{h}(ŷ_{0})))ψ̇_{t}(\overline{h}(ŷ_{0})) =
h′(ψ_{t}(\overline{h}(ŷ_{0})))f(ψ_{t}(\overline{h}(ŷ_{0}))) = (h′f)(ψ_{t}(\overline{h}(ŷ_{0}))) =
(f̂ ∘ h)(ψ_{t}(\overline{h}(ŷ_{0}))) = f̂(h(ψ_{t}(\overline{h}(ŷ_{0})))) = f̂(ψ̂_{t}(ŷ_{0})).

Eqn. (1.21) (which also implies (1.20) through (1.6)) can be proved in a similar way.

We shall need a tool from linear algebra.

THEOREM 1.9. Let N, N̂ ∈ N; N ≥ N̂, and let A ∈ R^{N̂×N̂}, B ∈ R^{N×N}, and let us suppose that the matrix M ∈ R^{N̂×N} has a full rank. If the relation

(1.23)  AM = MB

holds, then each eigenvalue of A is an eigenvalue of B, too.

Proof 1. Suppose the rows of M are orthonormal. If not, consider the decomposition M = QV with Q ∈ R^{N̂×N̂}, V ∈ R^{N̂×N} and VV^{⊤} = IN̂, i.e. with a matrix V having orthonormal rows obtained e.g. from the Gram–Schmidt orthogonalization procedure. The resulting matrix Q will be a triangular, invertible matrix because M is of the full rank. With this decomposition our assumption is AQV = QVB, or, equivalently, Q^{−1}AQ · V = V · B. Introducing the notations Ã := Q^{−1}AQ, M̃ := V, B̃ := B one obtains ÃM̃ = M̃B̃ with M̃ having orthonormal rows. (Alternatively, we could have said that a base transformation changed the rows of M to orthogonal vectors.)

Thus, our assumption can be taken as AM = MB, MM^{⊤} = IN̂. Let’s join N − N̂ more rows to M in such a way that

P := ( \begin{matrix}M \\ U \end{matrix} )

be orthogonal, i.e. PP^{⊤} = P^{⊤}P = (M^{⊤}M + U^{⊤}U =)I_{N}. Then we also have

( \begin{matrix}A 0 \\ UBP^{⊤} \end{matrix} ) P = ( \begin{matrix}A & 0 \\ UBM^{⊤} & UBU^{⊤} \end{matrix} ) ( \begin{matrix}M \\ U \end{matrix} ) = ( \begin{matrix}AM \\ UB(M^{⊤}M + U^{⊤}U) \end{matrix} )
= ( \begin{matrix}MB \\ UB \end{matrix} ) = PB.

This equality however means that the matrices ( \begin{matrix}A & 0 \\ UBM^{⊤} & UBU^{⊤} \end{matrix} ) and B are similar, thus the statement follows from the blocktriangular forms of these matrices, because the eigenvalues of B are the eigenvalues of A and of UBU^{⊤}.

Proof 2. Suppose M is of full rank. Then there exists \overline{M} ∈ R^{N×N̂} such that M\overline{M} = IN̂ holds. We show that no λ ∈ C which is not an eigenvalue of B can be an eigenvalue of A. This statement is equivalent with the original one.

If λ ∈ C is not an eigenvalue of B then there exists (B − λI_{N})^{−1}. We show that in this case there exists (A − λIN̂)^{−1} as well. Furthermore, we explicitely give this inverse as M(B − λI_{N})^{−1}\overline{M}. Let us check this statement using (1.23):

(A − λIN̂)(M(B − λI_{N})^{−1}\overline{M}) = (AM − λM)(B − λI_{N})^{−1}\overline{M}

= (MB − λM)(B − λI_{N})^{−1}\overline{M} = M(B − λI_{N})(B − λI_{N})^{−1}\overline{M}
= M\overline{M} = IN̂.

1.3. Notions of formal reaction kinetics. At the outset we consider only homogeneous isothermal systems. However, the models we introduce are capable of describing some systems outside this area.

Let us suppose we have a complex chemical reaction, i.e. a vessel containing N (a finite number of) chemical species X_{1}, X_{2}, ..., X_{N} and suppose R (a finite number again) reaction steps are taking place as follows:

(1.24)  ∑_{n=1}^{N} α(n, r)X_{n}⟶^{k_{r}} ∑_{n=1}^{N} β(n, r)X_{n} (r = 1, …, R).

Here, as usual, the positive real numbers k_{r} are the reaction rate constants, and the nonnegative integers α(n, r) and β(n, r) are the stoichiometric coefficients.

The order of the rth reaction step in (1.24) is defined to be ∑_{n=1}^{N} α(n, r), whereas the order of the reaction is the maximum of the orders of the steps: max_{r=1,2,...,R} ∑_{n=1}^{N} α(n, r).

A further specialization of first order reactions is obtained if

max_{r=1,2,...,R} ∑_{n=1}^{N} α(n, r) ≤ 1

also holds, i.e. if we only have reaction steps of the type

X_{n} ⟶ X_{p} X_{n} ⟶ O O ⟶ X_{n}.

Such a complex chemical reaction is called a compartmental system.

Sometimes we also include species denoted by the first letters of the alphabet like A and B, or by O as above in the definition of compartmental systems. These species are considered to have constant concentration (usually supposed to be 1) and are only included so as to make the examples more realistic. For example, we write

A + X ⟶ 2X

instead of writing

X ⟶ 2X .

In applications these species are the pools, as e.g. a buffer, or water etc. These species are also called external species as opposed to the internal ones about whose time evolution we are actually interested in.

The rth reaction step in (1.24) is said to be reversible if the reverse reaction also occurs among the reaction steps, otherwise it is irreversible. Reaction (1.24) is said to be reversible if all the reaction steps are reversible. If there is at least one irreversible step, then the reaction is said to be irreversible.

The reaction (1.24) is said to be mass-conserving if there exist positive real numbers ρ_{n} (n = 1, 2, …, N) for which

(1.25)  ∑_{n=1}^{N} α(n, r)ρ_{n} = ∑_{n=1}^{N} β(n, r)ρ_{n} (r = 1, …, R)

holds.

Real chemical reactions without approximation will surely have this property with ρ_{m} as the molecular weight of the species X_{m}. Approximate networks and networks devised to be investigated from theoretical points of view however will not necessarily fulfill this criterion, although we do not want to exclude these ”reactions” from our investigations at the outset.

The usual mass action type deterministic model of this reaction is the polynomial differential equation:

(1.26)  ẏ_{n} = ∑_{r=1}^{R}(β(n, r) − α(n, r))k_{r} ∏_{p=1}^{N} y_{p}^{α(p,r)} (n = 1, …, N)

with the concentrations of species (y_{n} := [X_{n}]) as dependent variables, and with time as independent variable. Equation (1.26) is usually referred to as the (induced) kinetic differential equation of reaction (1.24).

Here we cite a result we shall use below related to the form of kinetic differential equations. Although all kinetic differential equations are polynomial differential equations, it is not true that all polynomial differential equations can be considered as being induced by a reaction. E.g. let us consider the Lorenz-equation

(1.27)  ẋ = −σ(x − y) ẏ = rx − y − xz ż = xy − bz (σ, b, r ∈ R^{+}).

This equation is not a kinetic equation as it contains the term −xz. This term expresses the fact that y decreases in a process in which it does not take part. Such a term is said to express negative cross–effect. With this notion a simple characterization of kinetic differential equations within the class of polynomial ones has been given (Hárs, Tóth 1979) and applied many times (see e.g. Tóth, 1979; Dancsó and Farkas, 1989, Tóth and Hárs 1986a,b; Farkas and Noszticzius, 1992; Poland, 1993).

THEOREM 1.10. A polynomial differential equation is kinetic if and only if it does not contain terms expressing negative cross–effect.

2. General regularity properties. Here we treat regularity properties which are mainly connected with ’setting the stage’ for the solutions of a kinetic differential equation.

We shall show that some regularity properties of the right hand sides are conserved under lumping and might be destroyed under expanding. In the present section differentiability of f is dropped.

THEOREM 2.1. Let h be linear: h(ỹ) := Mỹ with some M ∈ R^{N̂×N} for all ỹ ∈ R^{N}. Suppose f fulfills the Lipschitz property of order α ∈ R^{+} with the constant L ∈ R^{+}. Then f̂ := Mf ∘ \overline{M} also has the Lipschitz property with the same order α ∈ R^{+}.

Proof. Let \overline{M} be any of the generalized inverses of M and let L ∈ R^{+} be the Lipschitz constant of f, i.e. the positive real number for which

|f(y_{1}) − f(y_{2})| ≤ L|y_{1} − y_{2}|^{α}

holds for all y_{1}, y_{2} ∈ R^{N}. Then, |f̂(ŷ_{1}) − f̂(ŷ_{2})| = |Mf(\overline{M}ŷ_{1}) − Mf(\overline{M}ŷ_{2})| ≤

|M||f(\overline{M}ŷ_{1}) − f(\overline{M}ŷ_{2})| ≤ |M|L|\overline{M}ŷ_{1} − \overline{M}ŷ_{2}|^{α} ≤
|M|L|\overline{M}|^{α}|ŷ_{1} − ŷ_{2}|^{α}.

Thus, one Lipschitz constant for f̂ might be L|M||\overline{M}|^{α}.

The best (:=smallest) Lipschitz constant L̂ can be achieved if one takes the infimum of the above expression for all the generalized inverses \overline{M}. An interesting question in linear algebra is when it is possible to attain L̂ = L because in general |M||\overline{M}| ≥ 1.

The statement could have also been formulated and proved for the local Lipschitz property.

EXAMPLE 2.1. Let f(y_{1}, y_{2}) := (sign(y_{1}), y_{2})^{T}, and let h(ỹ_{1}, ỹ_{2}) := ỹ_{2}. Then, f̂(˜ŷ) = ˜ ŷ, thus f̂ has the Lipschitz property, whereas f has not.

The structure of this example will be typical below in the following sense. We start from a system consisting of two qualitatively different parts, and we lump it by a projection to retain the ‘desired’ part only.

Now let us return to the general case of a nonlinear lumping function. For simplicity, the Lipschitz property is considered with α = 1.

THEOREM 2.2. Let us suppose that f fulfills the local Lipschitz property on R^{N} and let h : R^{N} ⟶ RN̂ be a nondegenerate function. Then the function f̂ defined by (1.8) also fulfills the local Lipschitz property.

Proof. Let y_{0} ∈ R^{n} be an arbitrarily fixed point, and U_{ε}(y_{0}) its neighborhood in which

|f(y_{1}) − f(y_{2})| ≤ L|y_{1} − y_{2}| for all y_{1}, y_{2} ∈ U_{ε}(y_{0})

holds with some L ∈ R^{+}.

The function f̂ does not depend on the choice of \overline{h}, therefore let us choose \overline{h} to be a continuously differentiable function for which \overline{h}(ŷ_{0}) = y_{0} holds with some ŷ_{0} ∈ RN̂.

Let ŷ_{1}, ŷ_{2} ∈ U_{δ}(ŷ_{0}). Because of the continuity of \overline{h}, δ can be chosen so small that

y_{1} := \overline{h}(ŷ_{1}) ∈ U_{ε/2}(y_{0}) y_{2} := \overline{h}(ŷ_{2}) ∈ U_{ε/2}(y_{0}).

The set \overline{U}_{ε/2}(y_{0}) being compact, h′, h′′ and f are bounded on it; let, say,

||h′|| < H_{1}, ||h′′|| < H_{2}, ||f|| < F.

On the set U_{δ}(ŷ_{0}) the function \overline{h} fulfills the local Lipschitz property with the Lipschitz constant

L̂ := max_{|ŷ−ŷ_{0}|≤δ} ||\overline{h}′(ŷ)||,

and obviously, h′ also fulfills the local Lipschitz property with the Lipschitz constant H_{2} on U_{ε/2}. Then

|f̂(ŷ_{1}) − f̂(ŷ_{2})| = |h′(\overline{h}(ŷ_{1}))f(\overline{h}(ŷ_{1})) − h′(\overline{h}(ŷ_{2}))f(\overline{h}(ŷ_{2}))|
= |h′(\overline{h}(ŷ_{1}))f(\overline{h}(ŷ_{1})) − h′(\overline{h}(ŷ_{1}))f(\overline{h}(ŷ_{2}))
+ h′(\overline{h}(ŷ_{1}))f(\overline{h}(ŷ_{2})) − h′(\overline{h}(ŷ_{2}))f(\overline{h}(ŷ_{2}))|
≤ |h′(\overline{h}(ŷ_{1}))f(\overline{h}(ŷ_{1})) − h′(\overline{h}(ŷ_{1}))f(\overline{h}(ŷ_{2}))|
+ |h′(\overline{h}(ŷ_{1}))f(\overline{h}(ŷ_{2})) − h′(\overline{h}(ŷ_{2}))f(\overline{h}(ŷ_{2}))|
≤ |h′(\overline{h}(ŷ_{1}))||f(\overline{h}(ŷ_{1})) − f(\overline{h}(ŷ_{2}))|
+ |h′(\overline{h}(ŷ_{1})) − h′(\overline{h}(ŷ_{2}))||f(\overline{h}(ŷ_{2}))|
≤ H_{1}LL̂|ŷ_{1} − ŷ_{2}| + F H_{2}L̂|ŷ_{1} − ŷ_{2}|
= L̂(H_{1}L + F H_{2})|ŷ_{1} − ŷ_{2}|.

Global Lipschitz property can be investigated in a similar way.

THEOREM 2.3. Let k ∈ N_{0}, f ∈ C^{k}(R^{N}, R^{N}) and let h ∈ C^{k+1}(R^{N}, RN̂). Then f̂ ∈ C^{k}(RN̂, RN̂).

Proof. Let us take \overline{h} from C^{k} which can always be done. Then, f̂ being a composition of C^{k}-functions (and being independent from the choice of \overline{h}), is a C^{k}-function itself.

Example 2.1 above also shows that expanding may not conserve the C^{k}-pro-perty.

THEOREM 2.4. Let f be a polynomial and let h be linear: h(ỹ) := Mỹ with some M ∈ R^{N̂×N} for all ỹ ∈ R^{N}. Then f̂ is also a polynomial with a degree not more than the degree of f.

Proof. Let \overline{M} be one of the generalized inverses of M. Then the explicit representation (1.8) specializes into

f̂ = Mf ∘ \overline{M}

in the present case, and this is a polynomial with a degree not more than the degree of f indeed.

EXAMPLE 2.2. The degree can be decreased in lumping as happens when the deterministic model

ẏ_{1} = −y_{1}^{2} + y_{1} ẏ_{2} = y_{1}^{2} + y_{2}

of the reaction

X_{1} ⟶ 2X_{1} ⟶ X_{1} + X_{2} X_{2} ⟶ 2X_{2}

is lumped by h(ỹ_{1}, ỹ_{2}) := ỹ_{1} + ỹ_{2} to obtain

˙ŷ = ŷ,

i.e. to reduce deg(f)=2 to deg(f̂)=1.

3. Invariant sets. INVARIANT SETS

THEOREM 3.1. Positively invariant (negatively invariant, invariant) sets under the flow generated by (1.1) are lumped into positively invariant (negatively invariant, invariant) sets of (1.5).

Proof. Let S be a positively invariant set of (1.1), and let Ŝ be defined as the set obtained from S by lumping: Ŝ := h(S). Let us denote the positive semiflow generated by (1.1) by ψ_{t} (t ≥ 0). Let ˜ ŷ ∈ Ŝ be an arbitrary point. By definition, there exists ỹ ∈ S such that ˜ ŷ = h(ỹ) holds. By Theorem 1.8

ψ̂_{t}(˜ŷ) = ψ̂_{t}(h(ỹ)) = h(ψ_{t}(ỹ)).

However ψ_{t}(ỹ) ∈ S because S is invariant, thus h(ψ_{t}(ỹ)) ∈ Ŝ as it is the result of applying the mapping h to an element of S.

The proof is similar for the negatively invariant and invariant cases.

COROLLARY 3.1. If (1.1) is a kinetic differential equation, then the first orthant (R_{0}^{+})^{N} is lumped into an invariant set of (1.5).

Proof. The statement is a consequence of a theorem by Vol’pert from ’72 (see e.g. Vol’pert and Hudyaev, 1985, or Érdi and Tóth, 1989) according to which the first orthant is a positively invariant set of kinetic differential equations.

The statement does not exclude the possibility of lumping the first orthant of, say, R^{3} into R, i.e. lumping an invariant set important from the point of view of reaction kinetics into an invariant set which is uninteresting or trivial from the point of view of reaction kinetics.

EXAMPLE 3.1 Let us lump the induced kinetic differential equation

ẏ_{1} = −y_{1} ẏ_{2} = y_{1}

of the reaction

X_{1}⟶^{1}X_{2}

by h(ỹ_{1}, ỹ_{2}) := −ỹ_{1} to obtain

˙ŷ = −ŷ.

Then the first orthant (R^{+})^{2} is lumped into R^{−}. Important classes of invariant sets are worth treating separately. EQUILIBRIA

THEOREM 3.2. Equilibria are lumped into equilibria.

Proof. Let y^{∗} be an equilibrium of the original system: f(y^{∗}) = 0, and let ŷ^{∗} := h(y^{∗}). Then,

f̂(ŷ^{∗}) = f̂(h(y^{∗})) = (h′f)(y^{∗}) = h′(y^{∗})f(y^{∗}) = 0.

Notice that h of the example below is not a lumping function in the sense we defined this notion.

EXAMPLE 3.2. The Lotka-Volterra model

A + X_{1}⟶^{1}2X_{1} X_{1} + X_{2}⟶^{1}2X_{2} X_{2}⟶^{1}A

induces the differential equation

ẏ_{1} = y_{1} − y_{1}y_{2} ẏ_{2} = y_{1}y_{2} − y_{2}

with the equilibria (0, 0) and (1, 1). Introducing h(ỹ_{1}, ỹ_{2}) := ỹ_{1}ỹ_{2} exp(−ỹ_{1} − ỹ_{2}) we obtain

˙ŷ = 0

as a reduced model with any real numbers as its equilibrium. The example is in no contradiction with our Theorem above as the Jacobian of h at (ỹ_{1}, ỹ_{2}) being

(ỹ_{2} exp(−ỹ_{1} − ỹ_{2}) − ỹ_{1}ỹ_{2} exp(−ỹ_{1} − ỹ_{2}), ỹ_{1} exp(−ỹ_{1} − ỹ_{2}) − ỹ_{1}ỹ_{2} exp(−ỹ_{1} − ỹ_{2}))

is degenerate (just) at the points (0, 0) and (1, 1). The example is the more special that h′ and f are orthogonal at each points of the state space.

The reason is that singular points of h′ emerge as equilibria of the model as it can be seen from the representation (1.8).

PERIODIC SOLUTIONS

THEOREM 3.3. Periodic solutions of (1.1) are lumped into periodic solutions of (1.5).

Proof. Let y be a periodic solution of (1.1) with the fundamental period T ∈ R^{+}: y(t + T) = y(t) for all t ∈ R. Then ŷ := h ∘ y is a periodic solution to (1.5) as well as

ŷ(t + T) = h(y(t + T)) = h(y(t)) = ŷ(t)

shows, and the fundamental period of ŷ can easily seen to be not larger than T.

Notice that the fact that h is a lumping function has not really been used here, any continuous function h maps a periodic functions into a periodic one.

The next example shows that lumped periodic solutions might be trivially periodic.

EXAMPLE 3.3. The Ivanova reaction (the only three-dimensional second order reaction with a linear first integral to show oscillation (Póta, 1985))

(3.1)  \begin{matrix}X_{1} + X_{2} & ⟶^{1}2X_{2} \\ X_{2} + X_{3} & ⟶^{1}2X_{3} \\ X_{3} + X_{1} & ⟶^{1}2X_{1} \end{matrix}

induces the differential equation:

(3.2)  \begin{aligned}ẏ_{1} & = −y_{1}y_{2} + y_{1}y_{3} \\ ẏ_{2} & = y_{1}y_{2} − y_{2}y_{3} \\ ẏ_{3} & = y_{2}y_{3} − y_{1}y_{3} \end{aligned}

whose solutions can shown to be oscillatory. However, the nondegenerate transformation h(ỹ_{1}, ỹ_{2}, ỹ_{3}) := ỹ_{1} + ỹ_{2} + ỹ_{3} provides the (alas!, kinetic) differential equation

˙ŷ = 0,

which only has constant (trivially periodic, if you like) solutions.

Using a projection as a lumping function it is easy to construct an example in which the fundamental period of the lumped model is strictly smaller than that of the original one.

It is not true that periodic solutions of the lumped model can only be obtained by lumping a periodic solution of the original.

EXAMPLE 3.4. The reaction

(3.3)  \begin{matrix}A + X_{1} & ⟶ 2X_{1} \\ 2A + X_{2} & ⟶ X_{1} + 2X_{2} \\ 2A + X_{1} & ⟶ X_{1} + 2X_{3} \\ A + X_{3} & ⟶ X_{2} + X_{3} \end{matrix}

induces the differential equation:

(3.4)  \begin{aligned}ẏ_{1} & = y_{1} + y_{2} \\ ẏ_{2} & = y_{2} + y_{3} \\ ẏ_{3} & = 2y_{1} \end{aligned}

with no periodic solution in the first orthant. However, with

h(ỹ_{1}, ỹ_{2}, ỹ_{3}) := (ỹ_{1} − ỹ_{2} ỹ_{1} − ỹ_{3})^{T}

one obtains the nonkinetic differential equation

(3.5) ˙ ŷ_{1} = ŷ_{2} ˙ ŷ_{2} = −ŷ_{1}

of the harmonic oscillator.

4. Linearized eigenvalues and Lyapunov functions.

Fundamental to the stability of equilibria are the following two statements. The first one is a generalization of the corresponding statement on C^{2}-conjugate vector fields, see e.g. Perko (1991; p. 173).

THEOREM 4.1. Suppose (1.1) is exactly lumpable to (1.5) by the nondegenerate transformation h : R^{N} ⟶ RN̂ and let y^{∗} be an equilibrium of (1.1), ŷ^{∗} := h(y^{∗}). Then, all the eigenvalues of the Jacobian of the lumped system taken at ŷ^{∗} are eigenvalues of the original system taken at y^{∗}.

Proof. Taking the derivative of (1.7) and substituting the equilibrium y^{∗} of the original equation we obtain:

f̂′(ŷ^{∗})h′(y^{∗}) = h′(y^{∗})f′(y^{∗}).

This however together with the fact that h′(y^{∗}) is of the full rank implies that each eigenvalue of f̂′(ŷ^{∗}) is an eigenvalue of f′(y^{∗}) as well according to Theorem 1.9. Moreover, notice that all y^{∗} being lumped into the same ŷ^{∗} will give the same eigenvalues for f̂′(ŷ^{∗}).

The statement is true for all nonequilibrium points in the case of linear lumping, as the same calculation shows. However, it is not true in general in the case of nonlinear lumping at nonequilibrium points.

COROLLARY 4.1. a. Sinks and sources are lumped into sinks and sources, respectively.

b. Suppose the right hand side f depends on a parameter taking on values in an interval in such a way that the equilibrium point is hyperbolic for all the values of the parameter. Then a Hopf bifurcation (see e.g. Schneider, 1980, Perko, 1991) of the lumped system cannot emerge at any of the parameter values in this interval.

c. Suppose the right hand side f̂ depends on a parameter taking on values in an interval in such a way that the equilibrium point undergoes Hopf bifurcation at a certain parameter value. Then the corresponding equilibrium of the original system undergoes Hopf bifurcation at the corresponding parameter value too.

Eigenvalues with positive real part might be lost under lumping resulting in stabilization of an unstable equilibrium, as can be shown by the projection method.

It might be instructive to collect the eigenvalues of all the examples into a table.

Example Equilibria Eigenvalues Lumped Lumped  Remarks
equilibria eigenvalues
2.1  (0, 0)  —  0  0  f ∉ C
2.2  (0, 0)  1,1  0  1  double eigenvalue
(1, −1)  −1, 1  0  1
3.1  (0, y_{2}^{0})  -1,0  0  −1  y_{2}^{0} ∈ R
3.2  (0, 0)  ±1  ŷ^{0}  0  ŷ^{0} ∈ R, h is degenerate
(1, 1)  ±i  ŷ^{0}  0  ŷ^{0} ∈ R, h is degenerate
3.3  (y_{1}^{0}, 0, 0)  ±y_{1}^{0}, 0  y_{1}^{0}  0  y_{1}^{0} ∈ R
(0, y_{2}^{0}, 0)  ±y_{2}^{0}, 0  y_{2}^{0}  0  y_{2}^{0} ∈ R
(0, 0, y_{3}^{0})  ±y_{3}^{0}, 0  y_{3}^{0}  0  y_{3}^{0} ∈ R
(y_{1}^{0}, y_{1}^{0}, y_{1}^{0}) ±y_{1}^{0}\sqrt{3}i, 0  3y_{1}^{0}  0  y_{1}^{0} ∈ R
3.4  (0, 0, 0)  ±i, 2  (0, 0)  ±i
4.1  (0, 0)  −1, 0  0  0  stability lost
4.2  (0, 0)  −1, 0  0  0
5.1  (0, 0)  0,0  0  0  double eigenvalue
5.2  —  —  —  —
7.1  (0, 0)  ±1  0  −1

The Lyapunov function of the original system can sometimes be used when investigating the lumped one.

THEOREM 4.2. Let us suppose (1.1) is exactly lumpable to (1.5) by the nondegenerate transformation h : R^{N} ⟶ RN̂ and let y^{∗} be an equilibrium of (1.1), U(y^{∗}) its neighborhood, ŷ^{∗} := h(y^{∗}). Let V be a continuously differentiable scalar valued function defined in the neighborhood U(y^{∗}), and suppose it only takes on positive values, except at the given equilibrium where it is zero. If V has the property that h(ỹ_{1}) = h(ỹ_{2}) implies V (ỹ_{1}) = V (ỹ_{2}), then there exists a continuously differentiable function V̂ : h(U(y^{∗})) ⟶ R for which V̂ ∘ h = V holds, such that V̂ will only take on positive values except the equilibrium ŷ^{∗}, where it is zero. The sign of the derivatives of the two functions V̂ and V along the corresponding trajectories of the corresponding differential equations are the same.

Proof. The statement follows from the following series of equalities:

V̂ (ŷ^{∗}) = V̂ (h(y^{∗})) = V (y^{∗}) = 0,

and

V̂ (˜ŷ) = V̂ (h(ỹ)) = V (ỹ) > 0, (˜ŷ ∈ h(U(y^{∗}))\{ŷ^{∗}}),

where ỹ ∈ U(y^{∗}) is any point for which ˜ ŷ = h(ỹ) holds. (ỹ = y^{∗} can only occur for ˜ ŷ = h(y^{∗}) = ŷ^{∗})

Finally, let ˜ ŷ ∈ h(U(y^{∗})) be arbitrary, and suppose that ỹ ∈ U(y^{∗}) is a point for which ˜ ŷ = h(ỹ) holds. Then,

˙V̂_{(1.5)}(˜ŷ) = ˙ V̂_{(1.5)}(h(ỹ)) = (V̂′f̂)(h(ỹ))
= V̂′(h(ỹ))(h′f)(ỹ) = V′(ỹ)f(ỹ) = V̇_{(1.1)}(ỹ).

(Here we used the usual notation V̇_{(k.n)} to denote the derivative of V along the trajectories of the system (k.n), which is sometimes also called the derivative of V with respect to the system (k.n).) Let us remark that a function V̂ fulfilling V̂ ∘ h = V also fulfills V̂ = V ∘ \overline{h} with any of the generalized inverses \overline{h}.

The statement above does not imply stability (asymptotic stability) of the lumped system in cases when the equilibrium point of the original system is stable (asymptotically stable), only under the condition of existence of the function V̂ mentioned above.

EXAMPLE 4.1 Let us lump the deterministic model

ẏ_{1} = −y_{1} ẏ_{2} = −y_{2}^{3}

of the reaction

X_{1} ⟶ A 3X_{2} ⟶ A

by h(ỹ_{1}, ỹ_{2}) := −2y_{2}^{2} to obtain the deterministic model

˙ŷ = ŷ^{2}

of the reaction

A + 2X̂ ⟶ 3X̂.

Although the equilibrium (0, 0) of the original equation was asymptotically stable, 0 is not a stable equilibrium of the lumped model. This is in no contradiction with the statement of Theorem 4.1, because the eigenvalues of the Jacobian of the original system at the equilibrium were −1 and 0, of which 0 remained for the lumped system.

The next example shows a case when a lumped Lyapunov function does exist.

EXAMPLE 4.2 Let us lump the deterministic model

ẏ_{1} = −y_{1} ẏ_{2} = −y_{2}^{3}

of the reaction

X_{1} ⟶ A 3X_{2} ⟶ A

by h(ỹ_{1}, ỹ_{2}) := ỹ_{2} to obtain the deterministic model

˙ŷ = −ŷ^{3}

of the reaction

3X̂ ⟶ A.

Investigation of a lumped system by Lyapunov’s method also yields some information on the original system.

THEOREM 4.3. Let us suppose (1.1) is exactly lumpable to (1.5) by the nondegenerate transformation h : R^{N} ⟶ RN̂ and let y^{∗} be an equilibrium of (1.1), U(y^{∗}) its neighborhood, ŷ^{∗} := h(y^{∗}). Let V̂ be a continuously differentiable scalar valued function defined in the neighborhood Û(ŷ^{∗}) := h(U(y^{∗})) ⊂ RN̂ of ŷ^{∗}, and suppose it only takes on positive values, except at the given equilibrium where it is zero. Then the function defined by V := V̂ ∘ h will only take on positive values on the set

Range(\overline{h}) ∩ U(y^{∗})

with any generalized inverse \overline{h} of h, except the equilibrium y^{∗}, where it is zero. The signs of the derivatives of the two functions V̂ and V along the corresponding trajectories are the same.

Proof. The statements come from the following series of equalities:

V (y^{∗}) = V̂ (h(y^{∗})) = V̂ (ŷ^{∗}) = 0,

and if

ỹ ∈ (Range(\overline{h}) ∩ U(y^{∗}))\\{y^{∗}\},

then there exists ˜ ŷ for which ỹ = \overline{h}(˜ŷ) holds, and thus

V (ỹ) = V̂ (h(ỹ)) = V̂ (h(\overline{h}(˜ŷ))) = V̂ (˜ŷ) > 0,

and finally

V̇_{(1.1)}(ỹ) = ((V̂ ∘ h)′ · f)(\overline{h}(˜ŷ))
= (V̂′ ∘ h)(\overline{h}(˜ŷ)) · h′(\overline{h}(˜ŷ)) · f(\overline{h}(˜ŷ))
= V̂′(h(\overline{h}(˜ŷ))) · f̂(h(\overline{h}(˜ŷ)))
= (V̂′f̂)(˜ŷ) = V̇_{(1.5)}(ỹ).

COROLLARY 4.2. If an equilibrium point of the lumped system is stable (asymptotically stable) then the corresponding equilibrium point of the original system is relatively stable.

The above arguments also show that everywhere defined Lyapunov functions are expanded in such a way that the cosine of the angle between their level surfaces and trajectories does not change sign. This means that Lyapunov functions can also be used for global investigations of the trajectories.

SENSITIVITIES

A similar argument to that used in Theorem 4.1 above shows that sensitivities can also be calculated in the case of equilibrium points (cf. Perko, 1991, p. 83).

THEOREM 4.4. Let us suppose (1.1) is exactly lumpable to (1.5) by h and let y^{∗} be an equilibrium of (1.1), ŷ^{∗} := h(y^{∗}). Let us choose \overline{h} in such a way that \overline{h}(ŷ^{∗}) = y^{∗} (which can always be done). Then the sensitivity of ŷ^{∗} with respect to the initial condition can be expressed as

(4.1)  R ∋ t ↦ exp ((h′f′)(y^{∗})\overline{h}′(ŷ^{∗})t) .

Proof. The sensitivity

\frac{∂ψ̂(t, ŷ)}{∂ŷ}

with respect to the initial condition of the dynamical system ψ̂ induced by (1.5) obeys the sensitivity equations:

(4.2)  \frac{d}{dt} \frac{∂ψ̂(t, ŷ)}{∂ŷ} = f̂′(ψ̂(t, ŷ))\frac{∂ψ̂(t, ŷ)}{∂ŷ}.

If the solution we consider is an equilibrium ŷ^{∗}, then we can easily express the Jacobian of the lumped equation to obtain:

f̂′(ŷ^{∗}) = (h′f′)(y^{∗})\overline{h}′(y^{∗}).

As this is the coefficient matrix of the linear (constant coefficient, in this case) differential equation (4.2), the statement follows.

5. Transient solutions. BLOW UP The definition (1.4) of the lumped solution shows that the next assertion is

true.

THEOREM 5.1. a. If the solutions are defined for all positive times, then the the solutions obtained by lumping are also defined for all positive times.

b. If a lumped solution blows up within a finite time, then the original one also blows up within a finite time.

Again, as in the case of periodic solutions, the fact that h is a lumping function has not really been used here, any continuous function h which is defined everywhere maps a function defined on an interval into another which is defined on the same whole interval.

This result may be especially useful for analyzing the explosive or thermal runaway properties of combustion systems through lumped schemes where asymptotic analysis will be much simpler.

Now we apply the statement above to a simple—artificial—example.

EXAMPLE 5.1. The reaction

(5.1)  \begin{matrix}A + 2X_{1} & ⟶^{1}X_{1} + 2X_{2} \\ 2A + X_{1} + X_{2} & ⟶^{1}4X_{2} \\ A + 2X_{2} & ⟶^{1}X_{2} + 2X_{1} \end{matrix}

induces the differential equation:

(5.2)  \begin{aligned}ẏ_{1} & = −y_{1}^{2} − y_{1}y_{2} + 2y_{2}^{2} \\ ẏ_{2} & = 3y_{1}y_{2} − y_{2}^{2} + 2y_{1}^{2}. \end{aligned}

However, with h(ỹ_{1}, ỹ_{2}) := ỹ_{1} + ỹ_{2} one obtains the kinetic differential equation

˙ŷ = ŷ^{2},

of the reaction

A + 2X̂⟶^{1}3X̂,

i.e. an equation well known to have solutions which blow up within a finite time. This means that the solutions to the original equation (5.2) also blow up, a fact which might have not been easy to show by other methods.

Using again a projection as a lumping function, everywhere defined lumped solutions can be created so that they are expanded into solutions with blowing up coordinates. Or, to put it the other way, examples might be constructed in which solutions with blowing up coordinates are lumped into everywhere defined solutions.

MONOTONICITY

We are trying to show by examples that some properties of transient solutions are not so easy to investigate.

EXAMPLE 5.2. The reaction

(5.3)  A ⟶ 2X_{1} A ⟶ X_{2}

induces the differential equation:

(5.4)  ẏ_{1} = 2 ẏ_{2} = 1,

whose solutions are strictly monotone increasing for all initial values y_{1}(0), y_{2}(0). However, with h(ỹ_{1}, ỹ_{2}) := ỹ_{2} − ỹ_{1} one obtains the nonkinetic differential equation

(5.5)  ˙  ŷ = −1

having only strictly monotone decreasing solutions for all the initial conditions.

The reaction of Example 3.4 induces the differential equation (3.4), whose solutions are monotone increasing (actually, strictly monotone increasing if only at least one of the initial concentrations y_{1}(0), y_{2}(0), y_{3}(0) is positive). However, for h as given one obtains the differential equation (3.5) of the harmonic oscillator which only has oscillatory solutions.

OLIGO-OSCILLATION (THE NUMBER OF LOCAL EXTREMA)

According to the theorem by Póta and Jost (Póta, 1981) all the concentrations in a compartmental system with N chemical species can have not more than N − 2 local extrema.

The consequence for our present purposes is that lumping a compartmental system by any method into a compartmental system with less species reduces the maximum number of local extrema the coordinates are able to produce.

Another consequence is that lumping anything into a compartmental system containing not more than 3 compartments the resulting system cannot show oligooscillation (multiple local extrema, see by Rábai et al., 1979 and Beck and Rábai, 1982).

Conversely, expanding can introduce oligo-oscillation into a system.

6. On the construction of lumping schemes.

Applicability of the method of lumping crucially depends on the possibility of constructing lumping schemes. In the case of linear lumping functions the problem has completely been solved (Li and Rabitz, 1989), whereas in the case of nonlinear lumping functions we have only been able to find some reformulations of the problem not providing the final answer (Li et al., 1994a) which however show the core of the problem. All the statements here have been proved in the two papers mentioned above.

The first results show that a linear lumping function can be found (if it exists at all) for arbitrary nonlinear equations.

THEOREM 6.1. Eqn. (1.1) is exactly lumpable to (1.5) by the linear transformation M ∈ R^{N×N̂} of full rank (with N ≤ N̂) if and only if

(i) there exists a nontrivial subspace M of R^{N} invariant under f′(y)^{⊤} and independent from y such that the rows of M are the basis vectors of M, and

(ii) for any generalized inverse \overline{M} of M and for all y ∈ R^{N} the eigenvectors of f′(y)^{⊤} and f′(\overline{M}My)^{⊤} corresponding to M are the same.

The next two statements provide methods to construct the nontrivial invariant subspace M mentioned above.

THEOREM 6.2. If the decomposition

f′(y)^{⊤} = ∑_{k=1}^{m} a_{k}(y)A_{k} (a_{k} : R^{N} ⟶ R, A_{k} ∈ R^{N×N})

holds and M is a nontrivial subspace of R^{N} invariant under all A_{k} (k = 1, 2, …, m) then M fulfils conditions (i)–(ii) of Theorem 6.1.

THEOREM 6.3. Let the real and complex eigenvalues of f′(y)^{⊤} be λ_{i}(y) (with multiplicity r_{i}, i = 1, 2, …) and σ_{j}(y) ± iτ_{j}(y) (with multiplicity s_{j}, j = 1, 2, …) respectively. Then

Ker {∏_{i}(f′(y)^{⊤} − λ_{i}(y)I_{N})^{r_{i}} ∏_{j}((σ_{j}^{2}(y) + τ_{j}^{2}(y))I_{N} − 2σ_{j}(y)f′(y)^{⊤}(f′(y)^{⊤})^{2})^{s_{j}} }

fulfils conditions (i)–(ii) of Theorem 6.1.

DEFINITION 6.1. The function ψ : R^{N} ⟶ R is a generalized eigenfunction of (1.1) if there exists a function Ω : R ⟶ R such that

Aψ = Ω ∘ ψ

holds, where A := ∑_{n=1}^{N} f_{n}(y)\frac{∂}{∂y_{n}} is the (linear) partial differential operator corresponding to the right hand side of eqn. (1.1).

THEOREM 6.4. Eqn. (1.1) is exactly lumpable to (1.4) by the lumping function h if and only if its coordinate functions are functions of N̂ functionally independent generalized eigenfunctions of (1.1).

How to determine the generalized eigenfunctions? These seem to be a really hard problem because it practically means the solution of a system of (linear) partial differential equations. We have offered two methods. The first one is based on Lie-algebraic methods. The second approach consists in finding the generalized eigenfunctions to the linear part of the original differential equation and then checking if the constructed lumping functions are appropriate for the original nonlinear system.

Finally let us also mention that the eigenfunctions of A (for which Aφ(y) = λ(y)ϕ(y) holds) can be used to transform A into a Jordan canonical form leading to (some) lumped systems—without the need to actually determine the generalized inverse \overline{h} of the lumping function h.

7. Open problems.

Now we present a collection of open problems.

In order to have a ‘good’ lumping function h we need convenient criteria to decide if h is nondegenerate. Such criteria can be found e.g. in Rădulescu and Rădulescu (1980).

The next series of examples and problems relate the right hand side of the equations. Let us consider Examples 3.4, 5.2 and the example by Li and Rabitz (1989, p. 1425) where an 8-dimensional kinetic differential equation has been (linearly) lumped into a 6-dimensional polynomial differential equation. In all the three examples kinetic differential equations have been lumped into nonkinetic ones. On the other side, nonkinetic differential equations might also be lumped into kinetic ones.

EXAMPLE 7.1. Let us consider the nonkinetic differential equation

(7.1)  ẏ_{1} = −y_{2} ẏ_{2} = −y_{1}

and let us lump it by h(ỹ_{1}, ỹ_{2}) := ỹ_{1} + ỹ_{2} to obtain the kinetic differential equation

˙ŷ = −ŷ

of the reaction

X̂⟶^{1}A.

The two sets of examples together show that expanding may also introduce and delete negative cross–effects thereby transforming kinetic differential equations into nonkinetic ones and vice versa.

PROBLEM 7.1. How to characterize those lumping/expanding schemes which transform kinetic differential equations into kinetic differential equations?

This question has been posed and answered for the case of linear lumping for first order reactions by Wei and Kuo (1969), and for second order reactions by Li (1984).

A related result has been obtained by Tóth and Hárs (1986a) who have shown for several examples that nonkinetic equations cannot be transformed into a kinetic equation with orthogonal transformations (:=special linear lumping functions with N = N̂).

Suppose we have a nonkinetic equation, can we expand it to a kinetic one, i.e. can we connect prototype schemes with qualitatively interesting behavior to real chemistry? Based upon earlier results by Korzukhin, Dancsó and Farkas (1989), Farkas and Noszticzius (1992) gave a method to approximately expand the harmonic oscillator into a kinetic differential equation. Poland (1993) did the same for the Lorenz equation (1.27).

(Approximate lumping in the above cases means that the lumped model is a QSSA of the original.)

A quite natural related question follows.

PROBLEM 7.2. Can a kinetic differential equation of a reversible reaction be lumped/expanded into a kinetic differential equation of an irreversible reaction, or conversely?

In order to answer this question a characterization of the right hand sides of differential equations induced by reversible reactions is needed (Simon and Tóth, in preparation).

The general answer to the questions above seems to be very hard, it might be easier to treat first and second order reactions. Another natural restriction— initially proposed and investigated by Wei and Kuo (1969)—might be that we only consider proper exact lumping, that is, the kind of lumping when variables are only collected without coupling between the groups.

In lumping theory it is quite common to use instead of a finite number of components an infinite number of them and reduce this infinite number to a finite one. As a result of lumping these systems “a first order system may be lumped into virtually any kinetic law” (Aris, 1969). Therefore, it would be extremely important to know how reaction orders are transformed.

It would be very important to know the effect of lumping on the existence of linear and nonlinear first integrals, especially on mass conservation.

In order to be able to apply the method used in Example 5.1 one should have general statements. Let us only consider another (more realistic, although not mass action kinetic) example (Scott and Tomlin, 1990) which has shown to be blowing up according to numerical calculations.

PROBLEM 7.3. Find a lumping function to reduce the dimension of the model:

α̇ = μe^{θ} − α(β^{2} + K)
β̇ = α(β^{2} + K) − β
θ̇ = δβ − γθ
(μ, K, δ, γ ∈ R^{+}are parameters.)

The next conjecture, the answer to which would make our remarks on compartmental systems and oligo-oscillation more useful, is: An appropriate choice of the rate constants and the initial conditions may assure that at least one of the concentrations does have as many local extrema as allowed.

It would be important to generally characterize the change in the number of equilibria under lumping, e.g. to provide a slight restriction (such as to exclude the orthogonality of h′ and f) on the lumping scheme so that this number do not increase. A related question is: in which neighborhood of a hyperbolic equilibrium do not change qualitatively the eigenvalues of the Jacobian?

As chaotic behavior can only occur in more than two dimensional systems, it will obviously not survive after lumping with N̂ ≤ 2. Expanding, however, can introduce chaotic behavior into a system. Results on Lyapunov exponents similar to the eigenvalues would provide an insight. It might also be interesting to lump the Lorenz equation into an equation with N̂ = 2.

As to approximate lumping our conjecture is that it also preserves the eigenvalues at least in the hyperbolic case—if the approximation is good enough.

In case one has a negative result (a counterexample) instead of being able to provide sufficient and necessary conditions, one should strive for obtaining positive results (sufficient conditions for the similarity of the lumped and the expanded systems) within restricted classes.

The questions treated in this paper should be (and sometimes has been) investigated for procedures related to lumping such as the application of the quasi-steady state assumption (QSSA), the use of the inertial manifold, of the center manifold, of the slow manifold. A systematic investigation of the effect of these reduction techniques onto kinetic equations and their solutions is missing as yet.

8. Summary and discussion.

In the present paper we have started a systematic investigation of the connections between the qualitative properties of kinetic differential equations and lumping.

Before doing so we provided a series of necessary and sufficient conditions of lumpability. Next we showed how the flow of the lumped system can be calculated from the original one.

We have shown that general regularity properties—like the Lipschitz property, or C^{k}-differentiability—of the right hand side are preserved. Invariant sets, equilibria and periodic solutions are lumped into invariant sets, equilibria and periodic solutions, respectively. The most far reaching result is that all the eigenvalues of the Jacobian of the lumped system at an equilibrium is an eigenvalue of the original system at the corresponding equilibrium. Using this fact it is possible to choose a lumping function so that we retain desired properties of the original system. For example, if our original system contains a number of negative eigenvalues, then we know that by lumping all these eigenvalues the asymptotic stability of the equilibrium point will generally be lost. Similarly, if we would like to retain oscillatory properties, then we must choose a lumping function which preserves the pure imaginary eigenvalues. The appropriate selection of submanifolds depends on the conservation of important qualitative properties.

We also presented examples on the effect of lumping on the transient behavior of solutions. In particular, we showed that if the lumped equation system shows solutions with a finite escape time, then the original system must exhibit the same property. This result has important consequences for lumped equation systems which show thermal runaway or explosion properties.

Acknowledgments. The authors acknowledge the support from the Department of Energy and from the National Scientific Research Found, Hungary (T-014480). Many thanks go to Drs. Gy. Póta, V. Ramakrishna and T. Szilágyi for providing important references. Dr. L. László and Prof. Z. Sebestyén were kind to provide the proofs of Theorem 1.9. Drs. B. Garay, P. L. Simon and L. Szili made many useful remarks when reading a previous version of our manuscript. One of the authors (J. T.) has also been supported by the National Scientific Research Found, Hungary, No. T-014480.

REFERENCES

R. Aris, Mathematical aspects of chemical reactions, Ind. Eng. Chem. Fundamentals, 61 (1969), pp. 17–29.

T. M. Beck and G. Rábai, Kinetics and mechanism of the oxidation of hydroxylamine by iodine, J. C. S. Dalton (1982) pp. 573–576.

A. BenArtzi and J. W. Helton, Ricatti partial differential equations for factoring nonlinear systems, (manuscript, 1988)

A. Dancsó and H. Farkas, On the “simplest” oscillating chemical system, Periodica Polytechnica, 33 (1989), pp. 275–285.

P. Érdi and J. Tóth, Mathematical Models of Chemical reactions: Theory and Applications of Deterministic Models, Princeton University Press, Princeton, NJ, 1989.

H. Farkas and Z. Noszticzius, Explosive, conservative, and dissipative systems and chemical oscillations, in Flow, Diffusion and Rate Processes, S. Sieniutycz and P. Salamon eds., Advances in Thermodynamics, Vol. 6, Taylor and Francis, London, 1992, pp. 299–335.

W. M. Getz and D. H. Jacobson, Sufficient conditions for finite escape times in systems of quadratic differential equations, J. Inst. Math. Appl., 19 (1977), pp. 377-383.

C. Godbillon, Géometrie différentielle et mécanique analitique, Herman, Paris, 1969.

W. Hahn, Theory and Application of Liapunov’s Direct Method, Prentice Hall Inc., Englewood Cliffs, N. J., 1963.

V. Hárs, V. and J. Tóth, On the inverse problem of reaction kinetics, in Qualitative Theory of Differential Equations, M. Farkas, ed., Coll. Math. Soc. J. Bolyai 30, North Holland, Amsterdam, 1981, pp. 363–379.

P. Hutchinson and D. Luss, Lumping of mixtures with many parallel first order reactions, Chem. Engng J., 1 (1970), pp. 129–135.

Y. Ijiri, Fundamental queries in aggregation theory, J. Am. Stat. Assoc., 66 (1971), pp. 153–164.

A. Isidori, Nonlinear Control Systems: An Introduction, M. Thoma, ed. Lecture Notes in Control and Information Sciences 72, Springer Verlag, Berlin, Heidelberg, New York, Tokyo, 1985.

Y. Iwasa, V. Andreasen and S. Levin, Aggregation in model ecosystems, I. Perfect aggregation, Ecol. Mod., 37 (1987), pp. 287–302.

Y. Iwasa, S. Levin and V. Andreasen, Aggregation in model ecosystems, II. Approximate aggregation, IMA J. of Appl. Math. Applied in Medicine and Biology, 6 (1989), pp. 1–23.

J. P. Keener, Principles of Applied Mathematics: Transformation and Application, Advanced book program, Addison Wesley, Redwood City, CA, 1988.

V. S. Korolyuk and A. F. Turbin, Mathematical Foundations of the State Lumping of Large Systems, Mathematics and Its Applications, Vol. 264, Kluwer Academic Publishers, Dordrecht, Boston, London, 1993.

G. Li, A lumping analysis in mono- or/and bimolecular reaction systems, Chem. Engng Sci., 29 (1984), pp. 1261–1270.

G. Li and H. Rabitz, A general analysis of exact lumping in chemical kinetics, Chem. Engng Sci., 44 (1989), pp. 1413–1430.

G. Li and H. Rabitz, A general analysis of approximate lumping in chemical kinetics, Chem. Engng Sci., 44 (1990), pp. 977–1002.

G. Li and H. Rabitz, New approaches to determination of constrained lumping schemes for a reaction system in the whole composition space, Chem. Engng Sci., 46 (1991a), pp. 95–111.

G. Li and H. Rabitz, Determination of constrained lumping schemes for nonisothermal first-order reaction systems, Chem. Engng Sci., 46 (1991b), pp. 583–596.

G. Li and H. Rabitz, A general lumping analysis of a reaction system coupled with diffusion, Chem. Engng Sci., 46 (1991c), pp. 2041–2053.

G. Li and H. Rabitz, A general analysis of lumping in chemical kinetics, in Kinetic and Thermodynamic Lumping of Multicomponent Mixtures, G. Astarita and S. I. Sandler eds., Elsevier Science Publishers B. V., Amsterdam, 1991d, pp. 49–62.

G. Li and H. Rabitz, The direct lumping approach: An application to a catalytic reforming model, Chem. Engng Sci., 48 (1993a), pp. 1903–1909.

G. Li, A. S. Tomlin, H. Rabitz and J. Tóth, Determination of approximate lumping schemes by a singular perturbation method, J. Chem. Phys., 99 (1993b), pp. 3562–3574.

G. Li, H. Rabitz and J. Tóth, A general analysis of exact nonlinear lumping in chemical kinetics, Chem. Engng Sci., 49 (1994a), pp. 343-361.

G. Li, A. S. Tomlin, H. Rabitz and J. Tóth, A general analysis of approximate nonlinear lumping in chemical kinetics I. Unconstrained lumping, J. Chem. Phys., 101 (1994b), pp. 1172–1187.

G. Li, A. S. Tomlin, H. Rabitz and J. Tóth, A general analysis of approximate nonlinear lumping in chemical kinetics II. Constrained lumping, J. Chem. Phys., 101 (1994c), pp. 1188–1201.

J. Łos and M. W. Łos (eds.), Mathematical Models in Economics, North Holland, Amsterdam, 1974.

N. K. Luckyanov, Yu. M. Svirezhev and O. V. Voronkova, Aggregation of variables in simulation models of water ecosystems, Ecol. Mod., 18 (1983), pp. 235–240.

D. G. Luenberger, Observing the state of a linear system, IEEE Trans. Mil. Electron., MIL-8 (1964), pp. 74–80.

L. Perko, Differential Equations and Dynamical Systems, Springer Verlag, New York, Berlin, Heidelberg, London, Paris, Tokyo, Hong Kong, Barcelona, 1991.

N. Peters and B. Rogg (eds.), Reduced Kinetic Mechanisms for Applications in Combustion Systems, Lecture Notes in Physics, M15, Springer Verlag, 1993.

L. M. Pismen, Dynamics of lumped chemically reacting systems near singular bifurcation points, Chem. Engng Sci., 39 (1984), pp. 1063–1077.

L. M. Pismen, Dynamics of lumped chemically reacting systems near singular bifurcation points, II. Almost Hamiltonian dynamics, Chem. Engng Sci., 40 (1985), pp. 905–915.

D. Poland, Cooperative catalysis and chemical chaos: a chemical model for the Lorenz equations, Physica D, 65 (1993), pp. 86–99.

Gy. Póta, On a theorem of overshoot-undershoot kinetics, React. Kinet. Catal. Lett., 17 (1981), pp. 35–39.

Gy. Póta, Irregular behavior of kinetic equations in closed chemical systems. Oscillatory effects, J. Chem. Soc. Faraday Transactions 2, 81 (1985) pp. 115–121.

M. Rădulescu and S. Rădulescu, Global inversion theorems and applications to differential equations, Nonlinear Analysis, Theory, Methods and Applications, 4 (1980), pp. 951–965.

Gy. Rábai, Gy., Bazsa and M. T. Beck, Design of reaction systems exhibiting overshoot-undershoot kinetics, J. Am. Chem. Soc., 101 (1979), pp. 6746–6748.

N. Rouchet, P. Habets and M. Laloy, Stability by Liapunov’s Direct Method, Springer Verlag, New York, 1977.

A. G. Schlijper and J. K. Drohm, Inverse lumping: Estimating compositional data from lumped information, SPE Reservoir Engineering, (Aug. 1988), pp. 1083– 1089.

K. R. Schneider, Andronov-Hopf-Bifurkation im Falle stetiger Parameterabhängigkeit, Z. Angew. Math. Mech., 60 (1980), pp. 19–24.

S. K. Scott and A. S. Tomlin, Period doubling and other complex bifurcations in chemical systems, Philos. Trans. Roy. Soc. London, Ser. A, 332 (1990), pp. 51–68.

J. Tóth, Gradient systems are cross-catalytic, React. Kinet. Catal. Lett. 12 (1979), pp. 253–257.

J. Tóth, and V. Hárs, Orthogonal transforms of the Lorenz- and Rösslerequation, Physica D, 19 (1986a), pp. 135–144.

J. Tóth and V. Hárs, Specification of oscillating chemical reactions starting from a given linearized form, Theor. chim. Acta, 70 (1986b), pp. 143–150.

J. Tóth, G. Li, H. Rabitz and A. Tomlin, Reduction of the number of variables in dynamic models, (in press, 1996).

F. Verhulst, Nonlinear differential equations and dynamical systems, Springer Verlag, Berlin, New York, 1990.

A. I. Vol’pert and S. I. Hudyaev, Analysis in Classes of Discontinuous Functions and the Equations of Mathematical Physics, Martinus Nijhoff Publishers, Dordrecht, Boston, Lancaster, 1985.

J. Wei and J. C. W. Kuo, A lumping analysis in monomolecular reaction systems, Ind. Eng. Chem. Fundamentals, 8 (1969) pp. 114–133.

S. Wiggins, Introduction to Applied Nonlinear Dynamical Systems and Chaos, Springer-Verlag, New York, Berlin, Heidelberg, London, Paris, Tokyo, Hong Kong, 1990.