Research Article  Open Access
ShinqJen Wu, ChengTao Wu, "SimulinkBased Analysis for Coupled Metabolic Systems", Applied Computational Intelligence and Soft Computing, vol. 2018, Article ID 8075051, 10 pages, 2018. https://doi.org/10.1155/2018/8075051
SimulinkBased Analysis for Coupled Metabolic Systems
Abstract
Stability analysis and dynamic simulation are important for researchers to capture the performance and the properties of underling systems. Ssystems 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 Ssystems to various situations. In this study, we first set up the respective block diagrams of Ssystems for modulebased simulation. We then derive reasonable theorems to examine the stability of Ssystems 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.
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. Ssystems [1, 2] which are composed of a set of coupled powerlawbased ODEs are popular and have potential for large scale systems because Ssystems are able (a) to capture dynamic behaviour of various biological systems such as gene regulation networks, metabolic pathways, and signal transduction cascades, (b) to identify net 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. Onestage 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 Ssystem modelling. A selfinteractive 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 eventinduced operations [5]. Seedinginspired chemotaxis genetic algorithms (SCGA) were designed to be attracted to an attracter purposely and then jump from it successfully through seedinginspired strategies and winnerchemotaxis population migration [6]. Kimura and collaborators introduced SVMbased linear programming classification and extracted gene interactive information from the classifiers through a genetic local search with a distanceindependent diversity control [7–9].
Some researchers divided Ssystem modelling into two stages (structure identification and parameter estimation). In this way a multiobjective optimization problem becomes two singleobjective optimization problems (twostage system identification). Ssystem 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 Ssystems 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 Ssystems, 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 Ssystems are available the research point is to theoretically analyse Ssystems 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 Ssystems to various environmental situations. We then theoretically derive the stability conditions of Ssystems and examine them with the proposed graph models.
2. SimulinkBased Dynamic Behaviour Prediction
The steady state values and the sensitivity at equilibrium points of Ssystems were theoretically derived through algebraic equations because of the special powerlaw structures [14, 15]. Voit and coworkers discussed the operating principles of Ssystems from a nominal equilibrium point to a new given 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 Ssystems through SoftCADS (software for sensitivity calculation) [17]. In this section we shall develop block diagrams of Ssystems 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.
Ssystems derived from biochemical system theory (BST) are used to describe interactive behaviour of metabolites. The rate change of constitutes (metabolites, proteins or genes) , , is the net influx minus the net efflux ( ). For a fully connected network the prototype of Ssystems is described aswhere m is the number of independent constitutes, α_{i} and β_{i} are the rate constants, and and denote the interactive strength from on (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 source and three dependent constituents () in Figure 1. The constitute is generated from and the production is inhibited by both and . The generated induces ’s generation and further the generation of . Based on this pathway, we write the respective Ssystem asWe 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 and an experimental condition (Case 8 in Table 1).

Case 2 (genetic branch pathways (4 genes) [19]). The branch pathway in Figure 4 has one constant source and four dependent constituents (). The constitute is generated from and the production is inhibited by . further generates both and . The comes from and the degradation is excited by. We have the respective Ssystems in 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 and the experimental condition (Case 3 in Table 2).

Case 3 (a small scale genetic network [20]). The genetic system in Figure 7 have three constant sources () and five dependent constitutes (). 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 with an initial condition set at

3. Dynamic Behaviour Analysis
The Ssystem in (1) is a highly nonlinear system. Nontrivial steady state values (equilibrium points) are obtained by setting the equilibrium flux Let the steady state of the dependent variables (), the independent variables (), the kinetic orders , and the rate constants . We havewhere denotes in equilibrium, with and , with and , with , with , and with . The solution uniquely exists when is nonsingular:In the case of is singular the system has multiple solutions. We can use MoorePenrose psuedoinverse to get the best fit (least square) solution: and are kineticorderrelated parameters, and b is a rateconstantrelated parameter. In other words, the steady state () depends on the reaction parameters () and the independent variable (). 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 () should be changed such that the steady states of dependent variables () 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 into dynamic behaviour directly. We concern with what 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 product are located at the lefthand plane, then the Ssystem in (1) is stable, wherein the steady state parameter E is the outer product of the equilibrium flux vector and equilibrium state inverse vector ; i.e., where the vector with , and denotes the dependent variable in equilibrium.
Proof. By the converse theorem [21], we know the stability of the nonlinear Ssystem at the neighbourhood of concurs with that of the respective linearized system. So we linearize the system in (1) around the equilibrium state vector as follows. The perturbed independent variable , varies with time byLet denote the perturbed variable () and . In equilibrium the influx equals the efflux for all institutes. The linearized system becomesWe further rewrite the linearized equation as a matrix form,where is a square matrix of dimension n with the element /. The linear system in (11) is exponentially stable if all of the eigenvalues of the system matrix have negative real parts. This completes the proof.
Theorem 2. If is positive definite, E is positive semidefinite (), and none of the element of is zero (), then the system matrix is positive definite. The Ssystem in (1) is unstable.
Proof. We shall derive Theorem 2 through the commutative property of the Hadamard product (denoted by “o”) and the following theorem [22]. “If A and B are positive semidefinite, then so is . In addition, If B is positive and A has no diagonal entry equal to 0, then is positive definite.” In the case of the system parameter , the system in (1) is unstable if the positive semidefinite matrix E has no diagonal entry which is equal to zero. The diagonal entry because that The concentration of the dependent variable () varies with time during the entire experiment. At any time instant, depends on the net minus result of the influx and the efflux . The equilibrium exists when a balance between the influx and efflux is achieved. is the products of the rate constants (positive) and the powerlaw functions of the involved state variables (the concentration of the involved institutes is positive). Therefore, the influx is always positive. So each diagonal element of E is nonzero () at the equilibrium state We then infer that the system matrix () is positive definite. In other words, all eigenvalues of the system matrix () of the linearized system in (11) are located at the righthand plane. The respective Ssystem is unstable.
Theorem 3. If is positive definite and E is negative definite or is negative definite and is positive definite, then the Ssystem in (1) is stable.
Proof. We know and . For all we have and implies and . We have through the theorem in [23]. “If both A and B are positive definite, then so is .” Through (12) we have for all . The linearized system in (11) is stable and the respective Ssystem in (1) is also stable through the converse theorem in [21]. The same is for the case and . This completes the proof.
Corollary. If is positive semidefinite and E is positive definite, and all of the selfinteractions are not balanced off for all of the dependent variables (), then the system matrix is positive definite. The Ssystem in (1) is unstable.
Proof. For , we have . In other words, each diagonal element of is nonzero. The system in (1) is unstable because 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 kineticorder parameter of the branch system is
. The eigenvalues of A_{d} is 0.2705 ± 0.5152i, 1.1902, 0.8188. Therefore, 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 instead of only. The same phenomena exist in the cascade pathway and the small scale network. Table 4 lists the estimated eigenvalues of the system matrix 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 lefthand plane. The linearized system is stable. By using Theorem 1 we know that the respective Ssystems are all stable (the eigenvalues of Cases 4 to 8 of Table 1 are not shown in Table 4).

Figure 10 shows the dynamic simulation of the Stype branch pathway system at three different experimental environments: the independent variables are , 0.3 and 0.9. The initial conditions are set at (black solid lines), (red dashed lines), and (blue dashed lines) for the case . The initial conditions are set at (red solid lines) and (blue dashed lines) for the case . The initial conditions are set at (black solid lines), (red dashed lines), and (blue dashed lines) for the case . 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 kineticorder parameter first. If is positive or negative definite ( or ) then Theorems 2 and 3 or Corollary can be used for stability analysis.
4. Conclusion
The generalized mass action model (another popular BSTbased model in biological systems) is composed of a set of coupled nonlinearly differential equations. However, each equation includes more than two powerlaw functions. For a generalized mass action model with influxes and effluxes, ( is the number of nonlinear differential equations), the number of potential phenotypes Two of the T terms are chosen as dominate terms to get a phenotype for stability analysis. In Ssystems the change rate of the expression level of variables (genes) equals to the net influx minus the net efflux. Under a fixed set of dependent variables only one phenotype exists when the kineticorderrelated matrix 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.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research was supported by the Ministry of Science and Technology of Taiwan, Grant no. MOST 1072221E212013.
References
 M. A. Savageau, Biochemical systems analysis: A study of function and design in molecular biology, Addison Wesley, Reading, Mass, USA, 1976.
 E. O. Voit, Computational analysis of biochemical systems: A practical guide for biochemists and molecular biologists, Cambridge University Press, Cambridge, UK, 2000.
 S.J. Wu, C.T. Wu, and J.Y. Chang, “Fuzzybased selfinteractive multiobjective evolution optimization for reverse engineering of biological networks,” IEEE Transactions on Fuzzy Systems, vol. 20, no. 5, pp. 865–882, 2012. View at: Publisher Site  Google Scholar
 S.J. Wu and C.T. Wu, “Computational optimization for Stype biological systems: cockroach genetic algorithm,” Mathematical Biosciences, vol. 245, no. 2, pp. 299–313, 2013. View at: Publisher Site  Google Scholar  MathSciNet
 S.J. Wu and C.T. Wu, “A bioinspired optimization for inferring interactive networks: Cockroach swarm evolution,” Expert Systems with Applications, vol. 42, no. 6, pp. 3253–3267, 2015. View at: Publisher Site  Google Scholar
 S.J. Wu and C.T. Wu, “Seedinginspired chemotaxis genetic algorithm for the inference of biological systems,” Computational Biology and Chemistry, vol. 53, pp. 292–307, 2014. View at: Publisher Site  Google Scholar
 S. Kimura, S. Nakayama, and M. Hatakeyama, “Genetic network inference as a series of discrimination tasks,” Bioinformatics, vol. 25, no. 7, pp. 918–925, 2009. View at: Publisher Site  Google Scholar
 S. Kimura, Y. Amano, K. Matsumura, and M. OkadaHatakeyama, “Effective parameter estimation for Ssystem models using LPMs and evolutionary algorithms,” in Proceedings of the 2010 6th IEEE World Congress on Computational Intelligence, WCCI 2010  2010 IEEE Congress on Evolutionary Computation, CEC 2010, Spain, July 2010. View at: Google Scholar
 S. Kimura, K. Matsumura, and M. OkadaHatakeyama, “Efficient parameter estimation for the inference of ssystem models of genetic networks: Proposition of further problem decomposition and alternate function optimization,” ChemBio Informatics Journal, vol. 11, no. 1, pp. 24–40, 2011. View at: Publisher Site  Google Scholar
 J. Sun, J. M. Garibaldi, and C. Hodgman, “Parameter estimation using metaheuristics in systems biology: A comprehensive review,” IEEE Transactions on Computational Biology and Bioinformatics, vol. 9, no. 1, pp. 185–202, 2012. View at: Publisher Site  Google Scholar
 E. O. Voit, “Biochemical Systems Theory: A Review,” ISRN biomathematics, vol. 2013, Article ID 897658, 53 pages, 2013. View at: Publisher Site  Google Scholar
 A. R. Chowdhury, M. Chetty, and N. X. Vinh, “Incorporating timedelays in SSystem model for reverse engineering genetic networks,” BMC Bioinformatics, vol. 14, no. 1, article no. 196, 2013. View at: Publisher Site  Google Scholar
 A. R. Chowdhury, M. Chetty, and R. Evans, “Stochastic Ssystem modeling of gene regulatory network,” Cognitive Neurodynamics, vol. 9, no. 5, pp. 535–547, 2015. View at: Publisher Site  Google Scholar
 C. T. Wu, S. J. Wu, and J. Y. Chang, “Computational analysis of Stype biological systems,” International Journal of Engineering Research and Applications, vol. 3, no. 1, pp. 1976–1987, 2013. View at: Google Scholar
 F. Shiraishi, E. Yoshida, and E. O. Voit, “An efficient and very accurate method for calculating steadystate sensitivities in metabolic reaction systems,” IEEE Transactions on Computational Biology and Bioinformatics, vol. 11, no. 6, pp. 1077–1086, 2014. View at: Publisher Site  Google Scholar
 Y. Lee, P.W. Chen, and E. O. Voit, “Analysis of operating principles with Ssystem models,” Mathematical Biosciences, vol. 231, no. 1, pp. 49–60, 2011. View at: Publisher Site  Google Scholar  MathSciNet
 K. Sriyudthsak, H. Uno, R. Gunawan, and F. Shiraishi, “Using dynamic sensitivities to characterize metabolic reaction systems,” Mathematical Biosciences, vol. 269, pp. 153–163, 2015. View at: Publisher Site  Google Scholar  MathSciNet
 K.Y. Tsai and F.S. Wang, “Evolutionary optimization with data collocation for reverse engineering of biological networks,” Bioinformatics, vol. 21, no. 7, pp. 1180–1188, 2005. View at: Publisher Site  Google Scholar
 O. R. Gonzalez, C. Küper, K. Jung, P. C. Naval Jr., and E. Mendoza, “Parameter estimation using simulated annealing for Ssystem models of biochemical networks,” Bioinformatics, vol. 23, no. 4, pp. 480–486, 2007. View at: Publisher Site  Google Scholar
 W. S. Hlavacek and M. A. Savageau, “Rules for coupled expression of regulator and effector genes in inducible circuits,” Journal of Molecular Biology, vol. 255, no. 1, pp. 121–139, 1996. View at: Publisher Site  Google Scholar
 M. Vidyasagar, Nonlinear Systems Analysis, Englewood Cliffs, NJ, USA, PrenticeHall, 2nd edition, 1993.
 R. A. Horn and C. R. Johnson, Topics in Matrix Analysis, Cambridge University Press, Cambridge, UK, 1991. View at: Publisher Site  MathSciNet
 R. A. Horn and C. R. Johnson, Matrix Analysis, Cambridge University Press, Cambridge, UK, 1985. View at: Publisher Site  MathSciNet
Copyright
Copyright © 2018 ShinqJen Wu and ChengTao Wu. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.