# Simulink-Based Analysis for Coupled Metabolic Systems

Authors: Shinq-Jen Wu; Cheng-Tao Wu
Journal: Applied Computational Intelligence and Soft Computing
 (2018)

---

## Abstract

Stability analysis and dynamic simulation are important for researchers to capture the performance and the properties of underling systems. S-systems have good potential for characterizing dynamic interactive behaviour of large scale metabolic and genetic systems. It is important to develop a platform to achieve timely dynamic behaviour of S-systems to various situations. In this study, we first set up the respective block diagrams of S-systems for module-based simulation. We then derive reasonable theorems to examine the stability of S-systems and find out what kinds of environmental situations will make systems stable. Three canonical systems are used to examine the results which are carried out in the Matlab/Simulink environments.


---

## Body

## 1. Introduction

A model in state space representation described as nonlinearly coupled ordinary differential equations (ODEs) is able to extract biologically information of underlying systems. S-systems [1, 2] which are composed of a set of coupled power-law-based ODEs are popular and have potential for large scale systems because S-systems are able (a) to capture dynamic behaviour of various biological systems such as gene regulation networks, metabolic pathways, and signal transduction cascades, (b) to identifynet interactive strength of constitutes, and (c) to possess good generation characteristics.System identification is reduced to infer the structure and estimate the parameters in the case that a given family of ODEs is chosen. One-stage identification carries out these two things at the same time. The identification becomes a multiobjective optimization problem. Extremely good performances in globally and locally searching abilities for computational modelling challenge researchers. We previously proposed various methods to smarten up the existing intelligent technologies for S-system modelling. A self-interactive learning was proposed to integrate an error performance, a skeleton structure index and a smooth evolution index [3]. Cockroach genetic algorithms (CGA) remained the globally searching abilities of GA, and incorporated the competition behaviour of cockroach swarms for food shortage into GA to improve the locally searching ability [4]. Additionally, we proposed cockroach swarm evolution (CSE) which was based on the cooperatively foraging behaviour of cockroach swarms and the completive behaviour and swarm migration were treated as event-induced operations [5]. Seeding-inspired chemotaxis genetic algorithms (SCGA) were designed to be attracted to an attracter purposely and then jump from it successfully through seeding-inspired strategies and winner-chemotaxis population migration [6]. Kimura and collaborators introduced SVM-based linear programming classification and extracted gene interactive information from the classifiers through a genetic local search with a distance-independent diversity control [7–9].Some researchers divided S-system modelling into two stages (structure identification and parameter estimation). In this way a multiobjective optimization problem becomes two single-objective optimization problems (two-stage system identification). S-system modelling also reduces to parameter estimation when the interactive relation of constitutes is partially known, or the static biological pathway is known. Therefore, parameter estimation for S-systems is required in the case that (a) system structures are inferred, (b) the relationship between genes and/or proteins is known, or (c) the qualitative pathways of underlying systems are known. Two review articles for metaheuristic developments in systems biology were published recently [10, 11].Chowdhury and coworkers currently introduced time delay into S-systems, wherein kinetic constants become real numbers instead of integer values [12]. They further proposed stochastic models for noise contaminated biological systems [13]. When the structure and parameters of S-systems are available the research point is to theoretically analyse S-systems in order to realize the performance, the properties, and the limitation of underlying systems. In this study, we propose a dynamic graph model. This model is able to timely predict the response of S-systems to various environmental situations. We then theoretically derive the stability conditions of S-systems and examine them with the proposed graph models.


## 2. Simulink-Based Dynamic Behaviour Prediction

The steady state values and the sensitivity at equilibrium points of S-systems were theoretically derived through algebraic equations because of the special power-law structures [14, 15]. Voit and coworkers discussed the operating principles of S-systems from a nominal equilibrium point to a newgiven equilibrium state through algebraic equations and linear programming [16]. They developed COSMOS (computation of sensitivity in model ODE systems) to estimate the steady state values and the sensitivity at steady state [15]. Sriyudthsak and coworkers used logarithm gain variables and approximation methods to get the approximate values of the dynamic sensitivity of S-systems through SoftCADS (software for sensitivity calculation) [17]. In this section we shall develop block diagrams of S-systems and use the diagrams to simulate dynamic behaviour of underlying systems. Dynamic stability analysis will be done in the next section and tested in Matlab/Simulink environments.S-systems derived from biochemical system theory (BST) are used to describe interactive behaviour of metabolites. The rate change of constitutes (metabolites, proteins or genes)xi, i=1⋯n, is the net influx minus the net efflux (vi+-vi- ). For a fully connected network the prototype of S-systems is described as(1)x˙i=vi+-vi-=αi∏j=1n+mxjgij-βi∏j=1n+mxjhij,wherem is the number of independent constitutes, αi and βi are the rate constants, and gij and hij denote the interactive strength from xj on xi (positive values for excitatory effect and negative values for inhibitory effect).Block diagrams are a kind of systems models, wherein the principal parts or functions are shown as blocks which are connected to each other through directed lines (showing signal flow). Simulink is one of the toolboxes in Matlab (MATrix LABoratory, a software developed by The MathWorks Company). Researchers are able to draw block diagrams of underlying systems and do simulation through the Simulink toolbox. In the following we draw the respective block diagrams of three biological systems in the Simulink environment and test the system responses to various experimental environments.Case 1 (cascade pathways (3 genes) [18]).
There are one constant sourcex4 and three dependent constituents (x1,x2,x3) in Figure 1. The constitute x1 is generated from x4 and the production is inhibited by both x2 and x3. The generated x1 induces x2’s generation and further the generation of x3. Based on this pathway, we write the respective S-system as(2)x˙1=10x2-0.1x3-0.05x4-5x10.5,x˙2=2x10.5-1.44x20.5,x˙3=3x20.5-7.2x0.5.We then draw the respective block diagram in the Simulink environment, as shown in Figure 2. Figure 3 is the simulation results for the system in an initial condition (x1,x2,x3)=(0.2,0.5,0.1) and an experimental condition x4=0.75 (Case 8 in Table 1).Table 1
Steady state values of the branch pathway.
Case 1 2 3 4 5 6 7 8 initial values x 1 (0) 1.6 1.6 2.7 2.7 1.6 0.4 2.2 2 x 2 (0) 1.6 1.6 1.6 0.4 1.6 2.7 2 1.2 x 3 (0) 0.4 0.4 2 0.4 0.4 0.4 1.6 2.3 x 4 (0) 1.6 1.6 0.4 1.6 1.6 2.7 2.3 0.8  Independent variables x 5 0.6 0.3 0.3 0.9 0.9 0.9 0.6 0.6  steady-state values x 1 0.3995605 0.2127742 0.2127742 0.5776509 0.5776509 0.5776509 0.3995605 0.3995605 x 2 2.006073 1.317966 1.317966 2.564889 2.564889 2.564889 2.006073 2.006073 x 3 2.228371 1.38912 1.38912 2.937986 2.937986 2.937986 2.228371 2.228371 x 4 0.1427537 0.0962826 0.0962826 0.1797378 0.1797378 0.1797378 0.1427537 0.1427537Figure 1
Cascade pathways [18].Figure 2
The block diagrams for the cascade pathway.Figure 3
(3 genes) simulation results of Case 8 in Table1.Case 2 (genetic branch pathways (4 genes) [19]).
The branch pathway in Figure4 has one constant source x0 and four dependent constituents (x1,x2,x3,x4). The constitute x1 is generated from x0 and the production is inhibited by x3. x1 further generates both x2 and x4. The x3 comes from x2 and the degradation is excited byx4. We have the respective S-systems in (3)x˙1=20x3-0.8x0-10x10.5,x˙2=8x10.5-3x20.75,x˙3=3x20.75-5x30.5x40.2,x˙4=2x10.5-6x0.8.Figure 5 is the block diagram of the genetic branch pathway shown in Simulink environment. Figure 6 is simulation results for the case that the initial condition (x1,x2,x3,x4)=(2.7,1.6,2,0.4) and the experimental condition x0=0.3 (Case 3 in Table 2).Table 2
Steady state values of the cascade pathway.
Case 1 2 3 4 5 6 7 8 initial values x 1 (0) 0.1 0.4 0.6 0.6 0.6 0.6 0.4 0.2 x 2 (0) 0.1 0.6 0.1 0.6 0.1 0.6 0.1 0.5 x 3 (0) 0.1 0.6 0.4 0.4 0.6 0.1 0.6 0.1  Independent variables x 4 0.75 0.9 0.6 0.6 0.9 0.75 0.75 0.75  steady-state values x 1 1.834666 2.428703 1.301562 1.301562 2.428703 1.834666 1.834666 1.834666 x 2 3.539094 4.684998 2.510729 2.510729 4.684998 3.539094 3.539094 3.539094 x 3 0.6144261 0.8133678 0.4358904 0.4358904 0.8133678 0.6144261 0.6144261 0.6144261Figure 4
Genetic branch pathways [19].Figure 5
The block diagrams for the genetic branch pathway.Figure 6
(4 genes) simulation results of Case 3 in Table 2.Case 3 (a small scale genetic network [20]).
The genetic system in Figure7 have three constant sources (x6,x7,x8) and five dependent constitutes (x1,x2,x3,x4,x5). The solid lines show the generation and degradation of constitutes. The inhibitory and excitatory behaviour are described as dash lines with (-) and (+) notations, respectively. Equation (4) describes the interaction of constitutes. We then draw a block diagram in Figure 8 and execute simulation for Case 8 in Table 3, as shown in Figure 9. The experiment was conducted at (x6,x7,x8)=(0.75,1,1.25) with an initial condition set at (x1,x2,x3,x4,x5)=(0.1,0.7,0.4,0.1,0.7).(4)x˙1=5x3x5-1x6-10x12,x˙2=10x12x7-10x22,x˙3=10x2-1x8-10x2-1x32,x˙4=8x32x5-1x6-10x42,x˙5=10x42x7-10x52,Table 3
Steady state values of the small scale network.
Case 1 2 3 4 5 6 7 8 initial values x 1 (0) 0.4 0.7 0.1 0.4 0.4 0.7 0.1 0.1 x 2 (0) 0.4 0.4 0.4 0.4 0.1 0.1 0.1 0.7 x 3 (0) 0.7 0.1 0.1 0.4 0.7 0.1 0.4 0.4 x 4 (0) 0.1 0.4 0.4 0.7 0.4 0.1 0.4 0.1 x 5 (0) 0.7 0.1 0.4 0.4 0.7 0.4 0.7 0.7  Independent variables x 6 1 1 1.25 1.25 1.25 0.75 0.75 0.75 x 7 0.75 1 1.25 1 1.25 1 0.75 1 x 8 0.75 1.25 1.25 0.75 1 0.75 1 1.25  steady-state values x 1 0.7517064 0.7476744 0.7760044 0.7718421 0.7617077 0.6509968 0.6995416 0.679307 x 2 0.6509968 0.7476744 0.8675993 0.7718421 0.8516151 0.6509968 0.6058208 0.679307 x 3 0.8660254 1.118034 1.118034 0.8660254 1 0.8660254 1 1.118034 x 4 0.8848579 1 1.037891 0.9085603 0.9634925 0.7663094 0.8848579 0.90856 x 5 0.7663094 1 1.160397 0.9085603 1.077217 0.7663094 0.7663094 0.90856Figure 7
Small scale genetic networks [20].Figure 8
The block diagrams for the small scale genetic network.Figure 9
(5 genes) simulation results of Case 8 in Table 3.


## 3. Dynamic Behaviour Analysis

The S-system in (1) is a highly nonlinear system. Nontrivial steady state values (equilibrium points) are obtained by setting the equilibrium flux v-i+=v-i-. Let the steady state of the dependent variables y-j=ln⁡x-j (j=1,…,n), the independent variables yj=ln⁡xj (j=n+1,…,n+m), the kinetic orders ghij=gij-hij, and the rate constants αβi=ln⁡(βi/αi). We have(5)∑j=1nghij·y-j+∑j=n+1n+mghij·yj=αβi,i=1,2,…,n,(6)Ad·Y¯d+AI·YI=b,where Y¯d denotes Yd in equilibrium, Ad=[ghij] with i=1,…,n and j=1,…,n, AI=[ghij] with i=1,…,n and j=n+1,…,n+m, Yd=[yi] with i=1,…,n, YI=[yi] with i=n+1,…,n+m, and b=[αβi] with i=1,…,n. The solution uniquely exists when Ad is nonsingular:(7)Y-d=Ad-1·b-Ad-1AI·YI.In the case of Ad is singular the system has multiple solutions. We can use Moore-Penrose psuedoinverse Ad† to get the best fit (least square) solution:(8)Y-d=Ad†·b-Ad†AI·YI.Ad and AI are kinetic-order-related parameters, andb is a rate-constant-related parameter. In other words, the steady state (X-d) depends on the reaction parameters (αi,βi,gij,hij) and the independent variable (XI). Tables 1–3 [14] show the steady state values for these three genetic networks. The branch pathway system achieves the same steady state values in Cases 1, 7, and 8 of Table 1 because their independent variables are the same, even that the initial conditions of the system are totally different. Cases 2 and 3 or Cases 4, 5, and 6 of Table 1 also show the same results. The system in Cases 1, 2, and 5 of Table 1 is operated at the same initial conditions but in different environments (independent variables). Therefore, the system approaches different steady state values. Table 2 shows the same phenomena for the cascade pathway system. As for the small scale network, totally different steady state values, as shown in Table 3, are obtained because the independent variables are totally different.The change of the values of independent variables denotes that a cell faces persistent changes in survival environments. Systems will show different dynamic evolution in response to such a change, which may be induced by such a stress environment as heat shock. Lee and coworkers [16] discuss that how independent variables (XI) should be changed such that the steady states of dependent variables (X-d)jump from a nominal value to a target value. Three cases (n=m, n<m, n>m) with respect to a unique solution, no solution and infinite solutions are considered. They proposed an algebraic approach to get approximate solutions for the no solution case through choosing some key dependent variables. They also adopted a numerical approach to get an optimal solution for the infinite solution case. In other words, they focus on the steady state algebraic relation (static behaviour) between dependent and independent variables. In this section we shall look intodynamic behaviour directly. We concern withwhat kinds of independent variables are preferred by biological systems from the viewpoints of system dynamics.As we know a biological system always operates at a steady state and will be temporarily deviated from the state. For example, blood serum undergoes short term changes for the intake and absorption of water and food and for kidney and liver’s operations. In a living system there exist regulatory mechanisms to effectively regulate various concentrations back to their nominal steady state levels. So what we concern is whether a biological system is locally stable. What kinds of independent variables can generate an equilibrium state such that a system can asymptotically stabilized to the state.Theorem 1.
If all of the eigenvalues of the Hadamard productAd∘E are located at the left-hand plane, then the S-system in (1) is stable, wherein the steady state parameterE is the outer product of the equilibrium flux vector V-+ and equilibrium state inverse vector X-^d; i.e., E=V-+⊗X-^d where the vector X-^d=[1/x-i] with i=1,…,n, and x-i denotes the dependent variable xi in equilibrium.Proof.
By the converse theorem [21], we know the stability of the nonlinear S-system at the neighbourhood of X-d concurs with that of the respective linearized system. So we linearize the system in (1) around the equilibrium state vector X-d as follows. The perturbed independent variable δxi=xi-x-i,i=1,…,n, varies with time by(9)ddtδxi=∂vi+∂Xd-∂vi-∂XdXd=X¯d=αi∑j=1ngijx-jgij-1·∏l=1l≠jnx-lgij·∏k=n+1n+mxkgij·δxj-βi∑j=1nhijx-jhij-1·∏l=1l≠jnx-lhij·∏k=n+1n+mxkhij·δxj=∑j=1ngijx-j-1·v-i+·δxj-∑j=1nhijx-j-1·v-i-·δxj.Let zi denote the perturbed variable (z=δxi) and ghij=gij-hij. In equilibrium the influx v-i+ equals the efflux v-i- for all institutes. The linearized system becomes(10)ddtzit=∑j=1nghij·v-i+x-j-1·zit,i=1,…,n.We further rewrite the linearized equation as a matrix form,(11)Z˙t=Ad∘EZt,where E=V-+⊗X-^d is a square matrix of dimensionn with the element eij=v-i+/x-j. The linear system in (11) is exponentially stable if all of the eigenvalues of the system matrix (Ad∘E) have negative real parts. This completes the proof.Theorem 2.
IfAd is positive definite,E is positive semidefinite (Ad>0,E≧0), and none of the element of X-d is zero (x-i≠0,i=1,…,n), then the system matrix Ad∘E is positive definite. The S-system in (1) is unstable.Proof.
We shall derive Theorem2 through the commutative property of the Hadamard product (denoted by “o”) and the following theorem [22].“If A and B are positive semi-definite, then so isA∘B. In addition, If B is positive and A has no diagonal entry equal to 0, thenA∘Bis positive definite.”
In the case of the system parameterAd>0, the system in (1) is unstable if the positive semidefinite matrixE has no diagonal entry which is equal to zero. The diagonal entry eii=v-i+/x-i because that E=V-+⊗X-^d. The concentration of thedependent variable (xi) varies with time during the entire experiment. At any time instant, xi depends on the net minus result of the influx vi+ and the efflux vi-. The equilibrium exists when a balance between the influx and efflux is achieved. vi+ is the products of the rate constants αi (positive) and the power-law functions of the involved state variables (the concentration of the involved institutes is positive). Therefore, the influx vi+ is always positive. So each diagonal element ofE is nonzero (eii≠0,i=1,…,n) at the equilibrium state x-i≠0,i=1,…,n. We then infer that the system matrix (Ad∘E) is positive definite. In other words, all eigenvalues of the system matrix (Ad∘E) of the linearized system in (11) are located at the right-hand plane. The respective S-system is unstable.Theorem 3.
IfAd is positive definite andE is negative definite or Ad is negative definite and E is positive definite, then the S-system in (1) is stable.Proof.
We knowAd∘E=[ghij·eij] and Ad∘(-E)=[-ghij·eij]. For all i=1⋯n we have(12)λiAd∘E=-λiAd∘-E.Ad>0 and E<0 implies Ad>0 and (-E)>0. We have Ad∘(-E)>0 through the theorem in [23].“If both A and B are positive definite, then so isA∘B.”
Through (12) we have λiAd∘E<0 for all i=1⋯n. The linearized system in (11) is stable and the respective S-system in (1) is also stable through the converse theorem in [21]. The same is for the case Ad<0 and E>0. This completes the proof.Corollary. If Ad is positive semidefinite andE is positive definite, and all of the self-interactions are not balanced off for all of the dependent variables (gii≠hii,i=1,…,n), then the system matrix Ad∘E is positive definite. The S-system in (1) is unstable.Proof.
Forgii≠hii,i=1,…,n, we have ghii≠0. In other words, each diagonal element of Ad is nonzero. The system in (1) is unstable because Ad∘E is positive definite through the theorem in [22], as shown in Proof of Theorem 2.Now we shall use the derived theorems to discuss the stability of our systems. The kinetic-order parameter of the branch system isA d = - 05 0 - 0.8 0 0.5 - 0.75 0 0 0 0.75 - 0.5 - 0.2 0.5 0 0 - 0.8. The eigenvalues ofAd is -0.2705 ± 0.5152i, -1.1902, -0.8188. Therefore, Ad is neither positive nor negative definite (semidefinite). Only Theorem 1 can be used to examine the stability of the branch system. We have to estimate the eigenvalues of Ad∘E instead of Ad only. The same phenomena exist in the cascade pathway and the small scale network. Table 4 lists the estimated eigenvalues of the system matrix (Ad∘E) with respect to various values of the independent variables (experimental environments) used in Tables 1–3. Table 4 shows that for each system all of the eigenvalues of the linearized system are located in the left-hand plane. The linearized system is stable. By using Theorem 1 we know that the respective S-systems are all stable (the eigenvalues of Cases 4 to 8 of Table 1 are not shown in Table 4).Table 4
Stability analysis for the three biological systems. The independent variables in thebranch pathway is set at x5 = 0.6 (Case 1), 0.3 (Case 2), and 0.9 (Case 3). The independent variable in the cascade pathway is set at x4 = 0.75 (Case 1), 0.9 (Case 2), and 0.6 (Case 3). The independent variables in the small-scale network are set at (x6, x7, x8) = (1, 0.75, 0.75) for Case 1, (1, 1, 1.25) for Case 2, and (1.25, 1.25, 1.25) for Case 3.
cases a branch pathway (4-gene) A cascade pathway (3-gene) A small-scale network (5-gene) λ i A d ∘ E stable λ i A d ∘ E stable λ i ( A d ∘ E ) stable  1 -1.0599 ±1.7013i yes -0.5543 yes -13.0199 yes -7.9502 ± 1.6192i -1.6468 -15.0341 -4.6200 -26.6061 -16.5117 ± 11.5849i  2 -1.2373±1.9465i yes -0.4818 yes -14.9535 yes -9.7292 ± 1.0894i -1.4313 -14.9535 -4.0154 -29.9070 -20 ± 14.1421i  3 -0.9651 ± 1.5704i yes -0.6581 yes -17.3520 yes -7.1134 ± 1.5568i -1.9552 -15.5201 -5.4851 -25.7730 -21.9829 ± 15.4717iFigure10 shows the dynamic simulation of the S-type branch pathway system at three different experimental environments: the independent variables are x5=0.6, 0.3 and 0.9. The initial conditions are set at 1.61.60.41.6 (black solid lines), 2.221.62.3 (red dashed lines), and 21.22.30.8 (blue dashed lines) for the case x5=0.6. The initial conditions are set at 1.61.60.41.6 (red solid lines) and 2.71.620.4 (blue dashed lines) for the case x5=0.3. The initial conditions are set at 2.70.40.41.6 (black solid lines), 1.61.60.41.6 (red dashed lines), and 0.42.70.42.7 (blue dashed lines) for the case x5=0.9. Simulation results demonstrate that the system always achieves equilibrium states even evolution starts at different initial conditions. This coincides with the prediction of Theorem 1. Figure 10 also implies that the system will be stabilized to another equilibrium state in response to environmental change from a nominal situation (Case 1) to another stress situation (Cases 2 or 3). As for other systems readers can estimate the kinetic-order parameter Ad first. If Ad is positive or negative definite (Ad>0 or Ad<0) then Theorems 2 and 3 or Corollary can be used for stability analysis.Figure 10
Simulation results for abranch pathway system under different experiment environments (independent variables x5 = 0.6, 0.3 and 0.9) with various initial conditions: the initial conditions are set at [1.61.60.41.6] (black solid lines), [2.221.62.3] (red dashed lines), and [21.22.30.8] (blue dashed lines) for the case x5 = 0.6, at [1.61.60.41.6] (red solid lines) and [2.71.620.4] (blue dashed lines) for the case x5 = 0.3, and at [2.70.40.41.6] (black solid lines), [1.61.60.41.6] (red dashed lines), and [0.42.70.42.7] (blue dashed lines) for the case x5 = 0.9.


## 4. Conclusion

The generalized mass action model (another popular BST-based model in biological systems) is composed of a set of coupled nonlinearly differential equations. However, each equation includes more than two power-law functions. For a generalized mass action model withLi influxes and Ri effluxes, i=1,…,n (n is the number of nonlinear differential equations), the number of potential phenotypes T=∏i=1nLi·Ri. Two of the T terms are chosen as dominate terms to get a phenotype for stability analysis. In S-systems the change rate of the expression level of variables (genes) equals to thenet influx minus thenet efflux. Under a fixed set of dependent variablesonly one phenotype exists when the kinetic-order-related matrix AD=[gij-hij]i,j=1i,j=n is nonsingular.Therefore, we do not need to construct a design space as researchers did in generalized mass action models. What we need to do is to examine the eigenvalues of the Hadamard product of a system matrix and the steady state parameters (the flux and the states). Through the block diagrams in Simulink environments we are able to simulate the system response to various experimental situations. In the future we shall further discuss the dynamic response to various perturbations (initial states, initial time, independent variables, and dependent variables) and the second order sensitivity (showing the parameter sloppiness) at any time instant or around the steady state.



---

*Source: 8075051-2018-12-02.xml*