Reduction of the Stoichiometric Submodel

Conservation relations and stoichiometric cycles are related to the left and right null space of the stoichiometric matrices, respectively. Consider the mole balances of a thermodynamic model

˙

c=N J+N_eJ_e.

Letb ∈null([N, N_e]^T), then for all flux vectorsJ and J_e the conservation relation b^T c= const holds. Conservation relations define an invariant manifold in the concentration space, and thus the system can be reduced to this manifold. The reduced system contains less components.

Letb ∈null(N), then the effect of a certain flux distributionJ on the change of concentrations

˙

c is indistinguishable from the effect of a flux distribution J +b. Thus, with respect to the concentrations c, the flux vector is redundant and the dynamics of the original system can be described by a reduced system with less fluxes.

4.2.1.1. Reduction of Conservation Relations

We transform the concentrations c by a partitioned transformation such that cˆ₂ is a constant quantity (c˙ˆ₂ = 0) and thuscˆ₂(ˆµ₁,µˆ₂,µˆ_e) = ˆc_2,0. We solve this algebraic equation in order to gain a function µˆ₂(ˆµ₁,µˆ_e). Finally, after applying this function to all occurrences of µˆ₂ in the first subsystem, we consider the first subsystem independently from the second subsystem. Thus, we gain a reduced description of the original system. The following corollary discusses the details of this procedure.

Corollary 4.21 (Reduction of conservation relations). Let M be a thermodynamic model (see Definition 3.1, p. 40). LetT_c = [T_c,1^T , T_c,2^T ]^T be a square and invertible matrix withT_c,2N = 0and T_c,2N_e = 0. Assume that the matrix T_c,2∂c/∂µ T_c,2^T is invertible for all µ∈Rⁱ⁰ and µ_e ∈R^i0,e. Let M˜ be a thermodynamic model with

N˜ =T_c,1N, ˜c(˜µ, µ_e) =T_c,1c(T_c,1^T µ˜+T_c,2^T µˆ₂(˜µ, µ_e), µ_e), S˜=S, R(˜˜ µ, µ_e) =R(T_c,1^T µ˜+T_c,2^T µˆ₂(˜µ, µ_e), µ_e), N˜_e =T_c,1N_e, µ˜₀ = inv^T_Λ−1

c (T_c,1)µ₀ where µˆ₂(˜µ, µ_e) is a solution of

T_c,2c(T_c,1^T µ˜+T_c,2^T µˆ₂, µ_e) = T_c,2c(µ₀, µ_e,0)

and Λ_c is a symmetric, invertible matrix with T_c,1Λ_cT_c,2^T = 0. Then, a solution of M can be reconstructed from a solution of the reduced system M˜ by

c=T_c⁻¹

c˜ T_c,2c(µ₀, µ_e,0)

, µ=T_c,1^T µ˜+T_c,2^T µˆ₂(˜µ, µ_e),

J = ˜J , ∆µ= ∆˜µ.

4. Transformation and Reduction

Proof. TransformM withˆc=T_cc(cf. §4.7, p. 53). The first subsystem (subscript 1) corresponds to the reduced systemM˜. Observe that c˙ˆ₂ = 0 and thus cˆ₂ =T_c,2c(µ₀, µ_e,0). From this we can compute a functionµˆ2(ˆµ1, µe)by solving the equationTc,2c(T_c,1^T µˆ1+T_c,2^T µˆ2, µe) =Tc,2c(µ0, µe,0). If the solution exists, it is locally unique because the matrix T_c,2∂c/∂µ T_c,2^T is invertible. This follows from the Inverse Function Theorem. If T_c,2∂c/∂µ T_c,2^T would be singular for all µ and µ_e, the model equations would not determine a locally unique trajectory forµˆ₂. If the function ˆ

µ₂(ˆµ₁, µ_e) does not exist or is not globally unique, the trajectory ofM does not exist or is not unique, respectively. The possible non-uniqueness of the chemical potentials will be discussed in §4.44 (p. 66). With the function µˆ2(ˆµ1, µe) we get µ= T_c,1^T µˆ1+T_c,2^T µˆ2(ˆµ1, µe). We use this relation to derive the functions R(ˆˆ µ₁, µ_e) and ˆc₁(ˆµ₁, µ_e). Observe that the first part of the equations (cˆ₁, µˆ₁) can now be considered independently of the second (ˆc₂,µˆ₂) and that cˆ₂ and ˆ

µ₂ can be reconstructed from ˆc₁ and µˆ₁. Thus, the subsystem Mˆ₁ is a reduced version of Mˆ. Replacing the subscript1 with a tilde (M˜ = ˆM₁) yields the expressions given in the corollary.

The initial value for the reduced chemical potentialsµˆ_1,0 = ˜µ_1,0 follows from the relations given in §4.8 (p. 53). The state variables of the original system can be reconstructed by using the inverse transformation matrixT_c⁻¹, ˆc₂ =T_c,2c(µ₀, µ_e,0)and the function µˆ₂(˜µ, µ_e).

§ 4.22 (Symmetry and positive semi-definiteness of R˜ and ∂˜c/∂µ˜). The matricesR˜ and ∂˜c/∂µ˜ of the reduced model are symmetric and positive semi-definite.

Proof. The resistance matrices of the reduced system and of the original system are equal at corresponding state and input vectors: R(˜˜ µ, µ_e) = R(T_c,1^T µ˜ +T_c,2^T µˆ₂(˜µ, µ_e)). Thus, R˜ is symmetric and positive semi-definite.

Differentiating the expression for˜c(˜µ, µ_e)in Corollary 4.21 yields for the derivative ofc(˜˜µ, µ_e) that ∂˜c/∂µ˜=T_c,1∂c/∂µ T_c,1^T +T_c,1∂c/∂µ T_c,2^T ∂µˆ₂/∂µ˜. The total differential of the conservation relations in Corollary 4.21 isTc,2∂c/∂µ T_c,1^T dµ˜+Tc,2∂c/∂µ T_c,2^T dˆµ2 = 0. This leads to∂µˆ2/∂µ˜ =

−(T_c,2∂c/∂µ T_c,2^T )⁻¹T_c,2∂c/∂µ T_c,1^T . Combining these results yields:

∂˜c/∂µ˜=T_c,1 (∂c/∂µ−∂c/∂µ T_c,2^T (T_c,2∂c/∂µ T_c,2^T )⁻¹T_c,2∂c/∂µ)

| {z }

L

T_c,1^T . Because the matrixL is symmetric, the derivative ∂˜c/∂µ˜ is symmetric.

IfLis positive semi-definite, the derivative∂c/∂˜ µ˜is positive semi-definite. However, from the expression forL given above it is not obvious if L is positive semi-definite. For this reason, an alternative expression forLthat allows to prove its positive semi-definiteness is derived. Assume initially that ∂c/∂µ is positive definite. Then, the matrices W = Tc,1Λc(∂c/∂µ)⁻¹ and Z = [W^T, T_c,2^T ] can be defined. LetL⁰ be the matrix L⁰ = Λ_cT_c,1^T (T_c,1Λ_c(∂c/∂µ)⁻¹Λ_cT_c,1^T )⁻¹T_c,1Λ_c. The matricesL⁰ andL are equal becauseZ is a invertible matrix andL Z =L⁰Z = (Λ_cT_c,1^T ,0). In this representation, it can be seen easily now that the matrix L⁰ = L is symmetric and positive definite. In the limit case where ∂c/∂µ is only positive semi-definite, the limit of L may also be only positive semi-definite.

4. Transformation and Reduction

§ 4.23 (Computation of T_c,1,T_c,2 and Λ_c). Let Y be a matrix with full rank and with span(Y) = \

µ∈Rⁱ⁰,µe∈R^i0,e

span ∂c

∂µ(µ, µ_e)

.

In all cases considered in this thesis, the matrix ∂c/∂µ is diagonal (e. g. in an ideal dilute solution) or the result of a linear transformation of a system with a diagonal matrix ∂c/∂µ. Then, the column space span(∂c/∂µ) is independent of µ and µe, and one may choose Y =

∂c/∂µ(µ^∗, µ^∗_e) with an arbitrary state vector µ^∗ and input vector µ^∗_e. If ∂c/∂µ is invertible, then Y can be chosen to be the identity matrix (Y = I). If ∂c/∂µ is diagonal, then Y can be chosen to be a diagonal matrix with diagonal elements Y_ii = 1 if ∂c_i/∂µ_i 6= 0 and Y_ii = 0 if ∂c_i/∂µ_i = 0 for i = 1. . . i₀. With the matrix Y defined above, the matrix Y^T ∂c/∂µ Y has always full rank. Using the matrix Y, the matrices T_c,1, T_c,2 and Λ_c with the properties demanded in Corollary 4.21 can be gained easily by computing kernel matrices: (1) Compute a kernel matrix X with X Y^T N = 0and X Y^TN_e = 0. The choice T_c,2 =X Y^T guarantees that T_c,2N = 0, T_c,2N_e = 0 and that T_c,2∂c/∂µ T_c,2^T =X Y^T ∂c/∂µ Y X^T has full rank. (2) Choose an arbitrary symmetric and invertibleΛ_c, e. g.Λ_c=I. (3) ComputeT_c,1as a right kernel matrix of T_c,2^T Λ_c with T_c,2^T Λ_cT_c,1 = 0.

§ 4.24 (Non-linear equation). The reduction involves the solution of the non-linear equation Tc,2c(µ0, µe,0) = Tc,2c(T_c,1^T µˆ1+T_c,2^T µˆ2, µe). Thus, the computational complexity of the reduction of conservation relations may be large.

§ 4.25 (Invariance of the entropy production). The entropy production of the internal fluxes J is invariant under the above described reduction method: σ[˜s] = ∆˜µ^TJ˜= ∆µ^T J =σ[s].

§ 4.26 (Invariance of the Gibbs energy). For dˆc₂ = T_c,2dc = 0, the differential of the Gibbs energy is invariant under the above described reduction method: dg = µ^T dc = (dµ˜^T Tc,1 + dˆµ₂T_c,2)dc =dµ˜^T d˜c. The information on the dependency of the Gibbs energy on the change of the concentration of the conserved moietydˆc₂ =T_c,2dc is lost during the reduction process.

Example 4.27 (Reduction of a conservation relation). Consider a reactionAB with mass-action kinetics in an ideal dilute solution in a closed system:

N = −1

1

, c(µ) =

c^◦ exp((µ_A−µ^◦_A)/(R^∗T)) c^◦ exp((µ_B−µ^◦_B)/(R^∗T))

, R(µ) =ρ R◦(µ_A/(R^∗T), µ_B/(R^∗T)).

The system contains the conservation relationc_A+c_B = const. We choose the transformation matrices for the reduction as

T_c,1 = −1 1

, T_c,2 = 1 1

withΛc=I. The variablecˆ2 =cA+cBis the concentration of the conserved moiety. The variable ˆ

c₁ = c_B−c_A describes the difference of the concentrations of B and A. Although ˆc₁ may be

4. Transformation and Reduction

negative, it is a concentration in the sense of Definition 3.1 (p. 40). We have µ_A = −ˆµ₁+ ˆµ₂ and µ_B = ˆµ₁ + ˆµ₂. From the conservation relation we get the following dependency of µˆ₁ and ˆ

µ2:

c_A,0+c_B,0 =c^◦ exp((−µˆ₁+ ˆµ₂−µ^◦_A)/(R^∗T)) +c^◦ exp((ˆµ₁+ ˆµ₂−µ^◦_B)/(R^∗T)).

Solving this equation forµˆ₂ and replacing µ˜= ˆµ₁ yields ˆ

µ₂(˜µ) = R^∗T log

c_A,0 +c_B,0

c^◦ · 1

exp((−µ˜−µ^◦_A)/(R^∗T)) + exp((˜µ−µ^◦_B)/(R^∗T))

. Substituting this expression into˜c= ˆc₁ =−c_A+c_B with c_i =c^◦ exp((µ_i−µ^◦_i)/(R^∗T))and into R=ρ R◦(µ_A/(R^∗T), µ_B/(R^∗T)) yields the transformed functions

˜

c(˜µ) = (c^◦_A+c^◦_B)·−exp((−˜µ−µ^◦_A)/(R^∗T)) + exp((+˜µ−µ^◦_B)/(R^∗T)) + exp((−˜µ−µ^◦_A)/(R^∗T)) + exp((+˜µ−µ^◦_B)/(R^∗T)) and

R(˜ˆ µ) = ρ R◦(−˜µ+ ˆµ₂(˜µ),µ˜+ ˆµ₂(˜µ)).

Thus, the original system with two compounds and simple functions c(µ) and R(µ) can be reduced to one with one compound but more complex functions ˜c(˜µ)and R(˜˜ µ).

4.2.1.2. Reduction of Stoichiometric Cycles

The reduction of stoichiometric cycles reduces the number of fluxes. It proceeds analogously to the reduction of conservation relations. While the reduction of conservation relations is based on the left null space of the stoichiometric matrices, the reduction of stoichiometric cycles uses the right null space.

We transform the fluxesJ by a partitioning transformation such that∆ˆµ₂ = 0 and such that Jˆ₂ vanishes from the mole balances.

Corollary 4.28 (Reduction of cycles). Let M be a thermodynamic model (see Definition 3.1, p. 40) and let TJ = [TJ,1, TJ,2] be a square and invertible matrix with N TJ,2 = 0 and S TJ,2 = 0.

Assume that the matrixT_J,2^T R(µ, µ_e)T_J,2^T is invertible for allµ∈Rⁱ⁰ and µ_e ∈R^i0,e. Let Λ_J be a symmetric, invertible matrix with T_J,2^T Λ_JT_J,1 = 0. Then, a trajectory ofM can be reconstructed from a trajectory of the reduced system M˜ with

N˜ =N T_J,1, S˜=S T_J,1, N˜_e =N_e, c(˜˜µ, µ_e) = c(˜µ, µ_e), µ˜₀ =µ₀, R(˜˜ µ, µ_e) = (T_J,1^T R(˜µ, µ_e)T_J,1)−(T_J,1^T R(˜µ, µ_e)T_J,2) (T_J,2^T R(˜µ, µ_e)T_J,2)⁻¹ (T_J,2^T R(˜µ, µ_e)T_J,1) where

c= ˜c, µ= ˜µ,

J = (T_J,1−T_J,2(T_J,2^T R T_J,2)⁻¹(T_J,2^T R T_J,1)) ˜J , ∆µ= inv_Λ⁻¹

J (T_J,1^T ) ∆˜µ.

4. Transformation and Reduction

Proof. TransformM with J =T_JJˆ(see §4.7, p. 53). Observe that∆ˆµ₂ = 0. BecauseT_J,2^T R T_J,2 is invertible, we can solve the equation T_J,2^T R T_J,1Jˆ₁ + T_J,2^T R T_J,2Jˆ₂ = ∆ˆµ₂ = 0 for Jˆ₂ and get Jˆ2 = −(T_J,2^T R TJ,2)⁻¹(T_J,2^T R TJ,1) ˆJ1. If T_J,2^T R TJ,2 would be singular, the model equation would not uniquely determine the fluxes Jˆ₂. The possible non-uniqueness of the fluxes will be discussed in Corollary 4.71 (p. 76). Entering the Jˆ₂ computed above into the equation T_J,1^T R TJ,1Jˆ1+T_J,1^T R TJ,2Jˆ2 = ∆ˆµ1 yields

(T_J,1^T R(˜µ, µ_e)T_J,1)−(T_J,1^T R(˜µ, µ_e)T_J,2) (T_J,2^T R(˜µ, µ_e)T_J,2)⁻¹ (T_J,2^T R(˜µ, µ_e)T_J,1)

| {z }

R˜

Jˆ₁

|{z}J˜

= ∆ˆµ₁

|{z}∆˜µ

. The formulas for the reconstruction of the original variables from the reduced ones follow from the expressions given in §4.8 (p. 53).

§ 4.29 (Symmetry and positive semi-definiteness of R˜ and ∂˜c/∂µ˜). The matricesR˜ and ∂˜c/∂µ˜ of the reduced model are symmetric and positive semi-definite.

Proof. The derivative∂˜c/∂µ˜is symmetric and positive semi-definite because ˜c(˜µ, µ_e) =c(˜µ, µ_e) and ∂c/∂µ is symmetric and positive semi-definite.

The symmetry of R˜ is obvious from the expression given in Corollary 4.28.

The resistance matrix R˜ resulting from the reduction of stoichiometric cycles and the deriva-tive∂˜c/∂µ˜resulting from the reduction of conservation relations (see §4.22, p. 59) have the same structure. For this reason, the proof of the positive semi-definiteness of R˜ is completely analog to the proof of the positive semi-definiteness of∂˜c/∂µ˜ in §4.22 (p. 59).

§ 4.30 (Computation of T_J,1, T_J,2 and Λ_J). Let Y be a matrix with full rank and with span(Y) = \

µ∈Rⁱ⁰,µe∈R^i0,e

span (R(µ, µ_e)).

In all cases considered in this thesis, the matrixR is diagonal or the result of a linear transfor-mation of a system with a diagonal matrixR. Then, the column space span(R)is independent of µ and µ_e, and one may choose Y = R(µ^∗, µ^∗_e) with an arbitrary state vector µ^∗ and input vector µ^∗_e. If R is invertible, then Y can be chosen to be the identity matrix (Y = I). If R is diagonal, then Y can be chosen to be a diagonal matrix with diagonal elements Yii = 1 if R_ii 6= 0 and Y_ii = 0 if R_ii = 0 for i = 1. . . i₀. With the matrix Y defined above, the matrix Y^T R Y has always full rank. Using the matrix Y, the matrices T_J,1, T_J,2 and Λ_J with the properties demanded in Corollary 4.28 can be gained easily by computing kernel matrices: (1) Compute a kernel matrixXwithN Y X = 0andS Y X = 0. The choiceT_J,2 =Y X guarantees that N T_J,2 = 0, S T_J,2 = 0 and that T_J,2^T R T_J,2 =X^T Y^T R Y X has full rank. (2) Choose an arbitrary invertible, symmetric ΛJ, e. g. ΛJ =I. (3) Compute TJ,1 as a right kernel matrix of T_J,2^T Λ_J with T_J,2^T Λ_JT_J,1 = 0.

§ 4.31(Invariance of entropy production). The entropy production is invariant under the above described reduction method: σ[˜s] =σ[s].

4. Transformation and Reduction

Proof. Using the equations for the reconstruction of the original states given in Corollary 4.28 we get for σ[s] = ∆µ^T J^T:

σ[s] = ∆˜µ^T (T_J,1^T Λ_JT_J,1)⁻¹T_J,1^T Λ_J(T_J,1−T_J,2(T_J,2^T R T_J,2)⁻¹(T_J,2^T R T_J,1))

| {z }

I

J˜=σ[˜s].

§ 4.32(Invariance of the Gibbs energy). The differential of the Gibbs energy is invariant under the above described reduction method: dg=µ^T dc= ˜µ^T d˜c=d˜g.

Example 4.33. The network in Example 3.11 (Equation 3.2, p. 42) contains a stoichiometric cycle. The matrices

TJ,1 =



 1 0 0 1 0 0



, TJ,2 =



 1 1 1



, ΛJ =





+1 −1 0

−1 0 1

0 1 1





fulfill the conditions T_J,2^T Λ_JT_J,1 = 0 and N T_J,2 = 0. Here,T_J,1 is chosen such that J˜₁ =J₁ and J˜₂ =J₂. This leads to the given matrix Λ_J. A reduction with these matrices yields a system with

N˜ =





−1 0 1 −1

0 1



, R˜ = 1 R₁ +R₂+R₃

R₁(R₂+R₃) −R₁R₂

−R₁R₂ R₂(R₁+R₃)

.

This system contains a minimal number of fluxes and corresponds to the system given in Equa-tion 3.3 (p. 42).

Im Dokument Thermokinetic modeling and model reduction of reaction networks (Seite 58-63)