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.
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 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. 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 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 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) , , is the net influx minus the net efflux ( ). For a fully connected network the prototype of S-systems 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 S-system 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 S-systems 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 S-system 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 Moore-Penrose psuedoinverse to get the best fit (least square) solution: and are kinetic-order-related parameters, and b is a rate-constant-related 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 left-hand plane, then the S-system 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 S-system 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 S-system 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 semi-definite, 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 power-law 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 right-hand plane. The respective S-system is unstable.
Theorem 3. If is positive definite and E is negative definite or is negative definite and is positive definite, then the S-system 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 S-system 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 self-interactions are not balanced off for all of the dependent variables (), then the system matrix is positive definite. The S-system 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 kinetic-order parameter of the branch system is
. The eigenvalues of Ad 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 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).
Figure 10 shows the dynamic simulation of the S-type 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 kinetic-order 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 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 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 S-systems 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 kinetic-order-related 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 107-2221-E-212-013.