Control Design for Signal Transduction Networks

Abstract

Signal transduction networks of biological systems are highly complex. How to mathematically describe a signal transduction network by systematic approaches to further develop an appropriate and effective control strategy is attractive to control engineers. In this paper, the synergism and saturation system (S-systems) representations are used to describe signal transduction networks and a control design idea is presented. For constructing mathematical models, a cascaded analysis model is first proposed. Dynamic analysis and controller design are simulated and verified.

Introduction

Recently, interdisciplinary studies are becoming popular; one of the most attractive studies is to combine both biological systems and control engineering.^1,2,3 Signal transduction networks of biological systems are characterized as high complexity because they are composed of many biochemical reactions. Typical modeling process for that kind of systems may involve the following issues: modeling, dynamic analysis and steady-state analysis.^4,5 However, the complexity of cellular signal transduction network is in general incomprehensible. Thus, an effective method to develop a mathematically equivalent model of the biochemical networks is highly desirable.

The synergism and saturation system (S-system) ⁶ has been a well-studied approach in modeling biochemical networks which characterizes the signal transduction networks. It was shown that the S-system representation in terms of ordinary differential equations (ODEs) is capable of capturing the behavior of biochemical dynamics. Applying logarithm on the state variables linearizes the state equations of the S-system at steady state. Based on the linearized S-system, it is possible to analyze and predict the S-system behavior rather than directly resorting to the original nonlinear model. In, ¹ a linearized S-system was derived and the robust stability analysis was conducted. In ⁷ , the properties of cascaded signal transduction pathway, robust and optimal design of system circuit and rule of gene regulation were introduced. Based on the steady-state analyses of the S-system model, a robust control method is proposed for biochemical networks via feedback and feedforward biochemical circuits was proposed by.^8,9 The proposed robust circuit design schemes provide a systematic method with applications in synthetic circuit design for biotechnological purpose.

The purpose of this paper is to model and analyze the signal transduction networks in biological systems and transform the mathematical model from the uncontrollable and nonlinear form to a controllable and linear form. We adopt the S-system and the Michaelis-Menten rate law to drive a mathematical model for signaling transduction networks. For simplifying the construction of the mathematical model, a method called ‘cascaded analysis model’ is proposed. The model is intended to be used for constructing a simplified mathematical model in the form of S-systems. This method avoids directly solving the entire mathematical model, which could be extremely complicated in structure. Rather, the problem can be broken down into smaller partitions to lessen computational burden.

Stability analysis is also presented. The Lyapunov stability theory is applied to determine stability of the S-system under Taylor expansion which could be used to estimate stability margin of the biological system without control. A new control design idea for the nonlinear S-systems is then attempted. By applying feedback linearization, ¹⁰ a proper coordinate transformation can be derived to establish an effective method for control of the nonlinear systems with smooth nonlinearities, which are nonlinear in their state variables but linear in their control variables. For S-systems, dependent variables can be chosen as state variables and independent variables as control inputs. On the basis of, ¹¹ one can convert the nonlinear system into a linearized controllable system. By choosing an appropriate controller, it's shown that the baseline steady state can be driven to the desired level in a given period.

Methods

Modeling of signal transduction networks

A signal transduction network includes many scaffolds which can be bound with molecules. Since the entire pathways which would influence reactions are, in general, too large to be conducted, an effective method for constructing a simplified mathematical model is highly desirable. To this aim, a method is attempted here to construct a simplified model.

First, consider the scaffold protein with n binding domains and each domain can be bound with one molecule. Also, define all states and pathways of the scaffold protein as follows

$$S=\sum _{i=0}^n\left(\begin{array}{c}n\\ i\end{array}\right)$$ (1) $$P=\sum _{i=0}^{n-i}\left(\begin{array}{c}n\\ i\end{array}\right)\left(\sum _{j=1}^{n-i}{2}^{j-1}\right)$$ (2)

with

$$\left(\begin{array}{c}n\\ i\end{array}\right)=\frac{n}{i!\left(n-i\right)!}$$ (3)

where S, P, and n denote, respectively, states, pathways and binding domains.

Second, consider the case where scaffold protein can combine with one molecule each time. Under this situation, one can neglect the redundant pathways. To simplify the overly complicated structure, the new pathways can be written as follows

$$P=\sum _{i=0}^n\left(\begin{array}{c}n\\ i\end{array}\right)\left(n-i\right)$$ (4)

Third, define states and molecules as the state variables x and reject the reactions on the pathways. Adopting the Michaelis-Menten rate law, ¹² each reaction can be represented as an ODE. Consequently, the general equations which describe the temporal changes in the biochemical system ¹³ can be formulated as

$$\begin{array}{ll}{\dot{x}}_i=V_i^{+}-V_i^{-},\hfill & i=1,2,\mathrm{\dots },n\hfill \\ V_i^{+}={\alpha }_i\prod _{j=1}^{n+m}{x}_j^{g_{ij}},\hfill & V_i^{-}={\beta }_i\prod _{j=1}^{n+m}{x}_j^{h_{ij}}\hfill \end{array}$$ (5)

where V_i⁺ =V_i⁺ (x₁, x₂, …, x_n, x_{n + 1}, …, x_{n + m}) and V_i^– = V_i^– (x₁, x₂, …, x _n , x_{n + 1}, …, x_{n + m}) are general functions of the dependent variables x₁, x₂, …, x_n and independent variables x_{n + 1}, x_{n + 2}, …, x_{n + m} α _i and β_i are rate constants; g_ij and h_ij are the kinetic orders.

Considering the steady state of the system, all rate constants and variables in (5) are given as nonzero and take the logarithm in (5). Defining y_i = ln x_i and arranging all terms for y_i to one side and other terms to another side gives

$$In\left(\frac{{\beta }_i}{{\alpha }_i}\right)=\sum _{j=1}^{n+m}{g}_{ij}{y}_j-\sum _{j=1}^{n+m}{h}_{ij}\hspace{0.17em}{y}_j,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,2,\mathrm{\dots },n$$ (6)

Next, define $b_i=\text{In}\left(\frac{{\beta }_i}{{\alpha }_i}\right)$ and $a_{ij}=g_{ij}-h_{ij}$ . and a_ij = g_ij – h_ij.

A general S-system with n dependent variables and m independent variables can then be characterized in the matrix form as

$$A_D{\overrightarrow{y}}_D=\overrightarrow{b}-A_I{\overrightarrow{y}}_I$$ (7)

where

$$\begin{array}{l}{A}_D=\left[\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ a_{n_1}& a_{n_2}& \cdots & a_{n_n}\end{array}\right],\hspace{0.17em}{\overrightarrow{y}}_D=\left[\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_n\end{array}\right],\hspace{0.17em}\overrightarrow{b}=\left[\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_n\end{array}\right],\hfill \\ A_I=\left[\begin{array}{cccc}{a}_{1,n+1}& a_{1,n+2}& \cdots & a_{1,n+m}\\ a_{2,n+1}& a_{2,n+2}& \cdots & a_{2,n+m}\\ ⋮& ⋮& \ddots & ⋮\\ a_{n,n+1}& a_{n,n+2}& \cdots & a_{n,n+m}\end{array}\right],{\overrightarrow{y}}_I=\left[\begin{array}{c}{y}_{n+1}\\ y_{n+2}\\ ⋮\\ y_{n+m}\end{array}\right]\hfill \end{array}$$

with the subscripts D and I meant dependent and independent variables respectively.

According to (7), one can obtain the steady states y_i, i =1, …, n given by

$${\overrightarrow{y}}_D=A_D^{-1}\left(\overrightarrow{b}-A_I{\overrightarrow{y}}_I\right)$$ (8)

using the pre-described procedures, the originally complicated system could be transformation into an analyzable form.

The Michaelis-Menten equation for biological systems has been applied to investigate the concentration change of metabolites in each pathway of biochemical networks, and the concentration change equations are further expressed as ordinary differential equations. However, it may cost significant computation time to analyze all cellular signal reactions and interactions which are not all important or critical to the signal transduction networks. How to remove the redundant parts is an issue.

To simplify the analysis, we propose a cascaded analysis model to analyze the system with a simpler structure. A molecule that combined with a scaffold protein is a basic reaction in the mathematical model. This reaction can be described, for example, as a signal transduction pathway in Figure 1. After estimating all parameters of the S-system, one can compute the output concentration x₁ at the steady state. One then cascades the output concentration with a new molecule to generate a new signal transduction pathway. For the same reason, one can cascade molecules to construct a complete mathematical model as shown in Figure 2. Repeating the steps, one can construct a mathematical model which is easier than constructing the model at a time.

Figure 1.

Basic reaction in biological system.

Figure 2.

Complete cascaded analysis model; the small block is the 1st cascaded layer; the large block denotes the 2nd cascaded layer.

The method is demonstrated by a signal transduction network model with one scaffold protein and two binding domains in Figure 3. From (1) and (2), the number of states and pathways are 4 and 5, respectively. We neglect the redundant pathways described by (4) to simplify the complete model. The number of pathways of the new model becomes 4. We define states (S) and molecules as state variables x_i and implement the pathways to reactions. We can then modify reactions and construct new signal transduction pathways as shown in Figure 4.

Figure 3.

Reduced model with two binding domains.

Figure 4.

Constructing the signal transduction network by signal transduction pathways.

On the basis of the signal transduction pathways in Figure 4, three independent variables (x₇, x₈, x₉), which are outside signal of the signal transduction networks, are introduced to construct the analyzable model (Fig. 5). Applying the cascaded analysis model and considering the top part in Figure 5, we separate the pathways into two parts and define a new variable z_i to substitute x_i as indicated in Figure 6.

Figure 5.

Adding independent variables to modify the signal transduction network

Figure 6.

Example of cascaded analysis model.

Consider part 1 in Figure 6, the system includes three dependent variables (z₁, z₂, z₃) and two independent variables (z₄, z₅), and the fluxes contain variables (V⁺, V^–). The S-system is built as

$$\begin{array}{l}{\dot{z}}_1={\alpha }_1{z}_2^{g_{12}}{z}_3^{g_{13}}-{\beta }_1{z}_1^{h_{11}},\\ {\dot{z}}_2={\alpha }_2{z}_4^{g_{24}}-{\beta }_2{z}_2^{h_{22}}{z}_3^{h_{23}},\\ {\dot{z}}_3={\alpha }_3{z}_5^{g_{35}}-{\beta }_3{z}_2^{h_{32}}{z}_3^{h_{33}},\\ z_4,z_5=\text{constant}\end{array}$$ (9)

From which one can get the steady states of all state variables by taking logarithm on (9):

$$\begin{array}{l}{g}_{12}\text{In}{z}_2+g_{13}\text{In}{z}_3=\text{In}\left(\raisebox{1ex}{${\beta }_1$}\!\left/ \!\raisebox{-1ex}{${\alpha }_1$}\right.\right)+h_{11}\hspace{0.17em}\text{In}\hspace{0.17em}{z}_1,\\ g_{24}\hspace{0.17em}\text{In}{z}_4=\text{In}\hspace{0.17em}\left(\raisebox{1ex}{${\beta }_2$}\!\left/ \!\raisebox{-1ex}{${\alpha }_2$}\right.\right)+h_{22}\hspace{0.17em}\hspace{0.17em}\text{In}\hspace{0.17em}{z}_2+h_{23}\hspace{0.17em}\text{In}\hspace{0.17em}{z}_3,\\ g_{35}\hspace{0.17em}\text{In}{z}_5=\text{In}\left(\raisebox{1ex}{${\beta }_3$}\!\left/ \!\raisebox{-1ex}{${\alpha }_3$}\right.\right)+h_{32}\hspace{0.17em}\text{In}\hspace{0.17em}{z}_2+h_{33}\hspace{0.17em}\text{In}\hspace{0.17em}{z}_3\end{array}$$ (10)

Using the notations defined in (7), one can determine the steady state ${\overrightarrow{y}}_D$ by solving

$$\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ 0& a_{22}& a_{23}\\ 0& a_{32}& a_{33}\end{array}\right]\hspace{0.17em}\left[\begin{array}{c}{y}_1\\ y_2\\ y_3\end{array}\right]=\left[\begin{array}{c}{b}_1\\ b_2-a_{24}{y}_4\\ b_3-a_{35}{y}_5\end{array}\right]$$ (11)

Similarly, the S-system model for part 2 in Figure. 6 can be constructed as follows

$$\begin{array}{l}{\dot{z}}_6={\alpha }_6{z}_7^{g_{67}}{z}_8^{g_{68}}-{\beta }_6{z}_6^{h_{66}},\\ {\dot{z}}_7={\alpha }_7{z}_9^{g_{79}}-{\beta }_7{z}_7^{h_{77}}{z}_8^{h_{78}},\\ {\dot{z}}_8={\alpha }_8{z}_{10}^{g_{8,10}}-{\beta }_8{z}_7^{h_{87}}{z}_8^{h_{88}},\\ z_9,z_{10}=\text{constant}\end{array}$$ (12)

The response time of each stage in the cascaded analysis model is governed by the degradation rate α_i and β_i of the protein at the stage of the cascaded analysis model. Using the cascaded analysis model, the original model can be replaced by a simplified one, which would be useful while constructing the signal transduction networks for analysis purpose.

Stability analysis

Most chemical reactions in biological systems operate at a steady-state level, and normal concentrations in biological systems are maintained by regulatory mechanisms that stabilize the steady states. Effective regulation makes the concentrations return to steady states after being effected by external stimulation.

In general, the power-law representation in biological systems can be considered as a canonical nonlinear system. By performing Taylor expansion, the power-law representation can be employed as a piecewise expression. It provides a global representation of which validity and accuracy can be governed. ¹⁴ On the other hand, the power-law representation can be employed as a local representation. Its accuracy within a neighborhood can be justified by investigating the effect resulting from the residual dynamics.

Consider, for instance, an S-system with two variables given by

$$\begin{array}{l}{\dot{x}}_1=F_1\left(x_1,x_2\right)={\alpha }_1{x}_2^{g_{12}}-{\beta }_1{x}_1^{h_{11}},\\ {\dot{x}}_2=F_2\left(x_1{x}_2\right)={\alpha }_2{x}_1^{g_{21}}-{\beta }_2{x}_2^{h_{22}}\end{array}$$ (13)

By applying the second-order Taylor expansion around the operation point (x₁₀, x₂₀) gives

$$F_1\left(x_1,x_2\right)=a_{11}{x}_1+a_{12}{x}_2+b_1+\left(c_1{x}_1^2+d_1{x}_2^2\right)$$

where

$$\begin{array}{l}{a}_{11}=h_{11}\left(h_{11}-2\right){\beta }_1{x}_{10}^{h_{{11}^{-1}}},\\ a_{12}=-g_{12}\left(g_{12}-2\right){\alpha }_1{x}_{20}^{g_{{12}^{-1}}},\\ b_1=F\left(x_{10},x_{20}\right)+\frac{1}{2}\left[-h_{11}\left(h_{11}-3\right)\times {\beta }_1{x}_{10}^{h_{11}}-g_{12}\left(g_{12}-3\right){\alpha }_1{x}_{20}^{g_{12}}\right],\\ c_1=\frac{-h_{11}\left(h_{11}-1\right){\beta }_1{x}_{10}^{h_{{11}^{-2}}}}{2},\\ d_1=\frac{g_{12}\left(g_{12}-1\right){\alpha }_1{x}_{20}^{g_{{12}^{-2}}}}{2}\end{array}$$

Similarly,

$$F_2\left(x_1,x_2\right)=a_{21}{x}_1+a_{22}{x}_2+b_2+\left(c_2{x}_1^2+d_2{x}_2^2\right)$$

where

$$\begin{array}{l}{a}_{21}=h_{22}\left(h_{22}-2\right){\beta }_2{x}_{10}^{h_{{22}^{-1}}},\\ a_{22}=-g_{21}\left(g_{21}-2\right){\alpha }_2{x}_{20}^{g_{21-1}},\\ b_2=F_2\left(x_{10},x_{20}\right)+\frac{1}{2}\left[-h_{22}\left(h_{22}-3\right)\times {\beta }_2{x}_{10}^{h_{22}}-g_{21}\left(g_{21}-3\right){\alpha }_2{x}_{20}^{g_{21}}\right],\\ c_2=\frac{-h_{22}\left(h_{22}-1\right){\beta }_2{x}_{10}^{h_{{22}^{-2}}}}{2},\\ d_2=\frac{g_{21}\left(g_{21}-1\right){\alpha }_2{x}_{20}^{g_{{21}^{-2}}}}{2}\end{array}$$

The linearized S-system becomes

$$\dot{x}\approx Ax+\Delta A\left(x\right)+b$$ (1.1)

where

$$\begin{array}{cc}A=\left[\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right],& \Delta A\left(x\right)=\left[\begin{array}{l}{c}_1{x}_1^2+d_1{x}_2^2\hfill \\ c_2{x}_1^2+d_2{x}_2^2\hfill \end{array}\right]\end{array}$$

and x = [x₁ x₂] ^T and b = [b₁ b₂] ^T ; ΔA(x) denotes the residual error.

Consider the stability of the linearized S-system, the residual error ΔA(x) satisfies

$$\begin{array}{cc}{\Vert \Delta A\left(x\right)\Vert }_2\le \alpha {\Vert x\Vert }_2,& \Delta A\left(x\right)\in {ℝ}^n\end{array}$$ (15)

where ${\Vert \cdot \Vert }_2$ denotes the Euclidean norm. Given an arbitrarily chosen Q = Q^T > 0 and A is stable, by the Lyapunov stability theory, there will exist a solution P = P^T > 0 to the following matrix equality:

$$A^{T}P+PA=-2Q$$ (16)

Using the Rayleigh principle it can be shown that the whole system of (14), without the bias term, would be asymptotically stable provided that: ¹⁵

$$\alpha <\frac{{\lambda }_{\min }\left(Q\right)}{{\lambda }_{\max }\left(P\right)}$$ (17)

It should be noted that Taylor expansion will not change stability of the system. The purpose here is to estimate stability margin of the biological system without control.

Control design

Feedback linearization is popular for nonlinear control designs. ¹⁶ As it was shown by ¹¹ that a nonlinear system can be casted into a linearized controllable system. On the basis of a proper coordinate transformation, feedback linearization establishes a convenient tool for the control design of the nonlinear systems. This form is applicable for the biochemical systems given as follows

$$\dot{x}\left(t\right)=f\left(x\right)+l\left(x\right)u$$ (18)

where x ∈ ℝ ⁿ and u ∈ ℝ ^m .

One can formulate the feedback linearization problem as follows. Define the nonlinear feedback control law as:

$$u=A\left(x\right)+B\left(x\right)v$$ (19)

where A(x) is an m-dimensional vector field, B(x) is an m × m matrix and v is an m-dimensional vector.

Define a coordinate mapping as

$$\psi :z=\psi \left(x\right)$$ (20)

such that the affine nonlinear control system (18) is transformed into a linear controllable system. For instance, by considering the nominal system (18), there exists a coordinate transformation as follows

$$\left[\begin{array}{c}{z}_1\\ z_2\\ ⋮\\ z_n\end{array}\right]=\left[\begin{array}{c}{\psi }_1\left(x\right)\\ {\psi }_2\left(x\right)\\ ⋮\\ {\psi }_n\left(x\right)\end{array}\right]=\left[\begin{array}{c}h\left(x\right)\\ L_{f}h\left(x\right)\\ ⋮\\ L_f^{n-1}h\left(x\right)\end{array}\right]$$ (21)

where the Lie derivative $L_{f}h\left(x\right)=\frac{\partial h\left(x\right)}{\partial x}f$ is the direction derivative of h along the direction of f and one can find h(x) satisfying

$$\begin{array}{l}\begin{array}{cc}{L}_1{L}_f^{k-1}h\left(x\right)=0,& k=1,2,\mathrm{\dots },n-1\end{array}\\ L_1{L}_f^{n-1}h\left(x\right)\ne 0\end{array}$$ (22)

so that (18) can be transformed into a linear form as follows

$$\dot{z}=Az+Bv$$ (23)

where

$$A=\begin{array}{cc}\left[\begin{array}{ccccc}0& 1& 0& \cdots & 0\\ 0& 0& 1& \cdots & 0\\ ⋮& ⋮& \cdots & \ddots & ⋮\\ 0& 0& 0& \cdots & 1\\ 0& 0& 0& \cdots & 0\end{array}\right],& B=\left[\begin{array}{c}0\\ 0\\ ⋮\\ 0\\ 1\end{array}\right]\end{array}$$

and one can choose a control law as follows

$$u=\frac{-L_f^{n}h\left(x\right)}{L_l{L}_f^{n-1}h\left(x\right)}+\frac{1}{L_l{L}_f^{n-1}h\left(x\right)}v$$ (24)

where we have used the following definition ¹⁷

$$\begin{array}{cc}{L}_l{L}_f^{i-1}h\left(x\right)=\sum _{j=1}^n\frac{\partial L_f^{i-1}}{\partial x_j}l\left(x\right),& i=1,2,\mathrm{\dots },n\end{array}$$

and

$$\begin{array}{l}{L}_f^{0}h\left(x\right)=h\left(x\right),\\ \begin{array}{cc}{L}_f^{k}h\left(x\right)=L_f{L}_f^{k-1}h\left(x\right)=〈d{L}_f^{k-1},f\left(x\right)〉,& k\ge 1\end{array}\end{array}$$

Next, consider the following generalized representation of the nonlinear system with uncertainties as follows

$$\dot{x}=f\left(x\right)+\Delta f\left(x\right)+l\left(x\right)u$$ (25)

There exists smooth function Δf*(x) in ℝ ⁿ such that the uncertainties in (18), for all x ∈ ℝ ⁿ , satisfy the matching condition:

$$\Delta F\left(x\right)=l\left(x\right)\Delta f^{\ast }\left(x\right)$$ (26)

Substituting the matched uncertainties (26) into (18), by the feedback linearization, it can be transformed into

$$\begin{array}{l}{\dot{z}}_1=z_2+\left(L_{l}h\right)\Delta f^{\ast }+\left(L_{l}h\right)u,\\ {\dot{z}}_2=z_3+\left(L_l{L}_{f}h\right)\Delta f^{\ast }+\left(L_l{L}_{f}h\right)u,\\ ⋮\\ {\dot{z}}_n=L_f^{n}h+\left(L_l{L}_f^{n-1}h\right)\Delta f^{\ast }+\left(L_l{L}_f^{n-1}h\right)u\end{array}$$ (27)

Substituting (22) and (24) into (27), one can finally transform (27) into the following form:

$$\begin{array}{l}{\dot{z}}_1=z_2,\\ {\dot{z}}_2=z_3,\\ ⋮\\ {\dot{z}}_n=v+\left(L_l{L}_f^{n-1}h\right)\Delta f^{\ast }\end{array}$$ (1.1)

Assume that v is given by

$$v=-Kz$$ (29)

where K is a constant row vector. Substituting (29) into (23) yields

$$\dot{z}=A_{c}z+B\Delta B\left(z\right)$$ (30)

where

$$\begin{array}{cc}{A}_c=\left[\begin{array}{ccccc}0& 1& 0& \cdots & 0\\ 0& 0& 1& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \cdots & 1\\ -k_1& -k_2& -k_3& \cdots & -k_n\end{array}\right],& B=\left[\begin{array}{c}0\\ 0\\ ⋮\\ 0\\ 1\end{array}\right]\end{array}$$

and $\Delta B=\left(\left(L_l{L}_f^{n-1}h\right)\Delta f^{\ast }\right)\circ {\psi }^{-1}\left(z\right)$ where ○ denotes the composition of functions. Thus the function $\left(L_l{L}_f^{n-1}h\right)\Delta f^{\ast }$ can be transformed from x-domain to z-domain.

In the first stage of control design, the gain vector K should be chosen so that A_c would be stable. Stability analysis of the uncertain control system is next carried out by using Lyapunov stability theory.

To proceed, a Lyapunov candidate function is defined as follows

$$V\left(z\right)=z^{T}Pz$$ (31)

where P = P^T > 0.

Taking derivative with respect to time gives

$$\dot{V}\left(z\right)=z^T\left(A_c^{T}P+P{A}_c\right)z+2z^{T}PB\Delta B\left(z\right)$$ (32)

As A_c has been a stable matrix, for any Q = Q^T > 0, there is a unique symmetric positive definite solution P to the following Lyapunov matrix equation:

$$A_c^{T}P+P{A}_c=-Q$$ (33)

Substituting (33) into (32) gives

$$\dot{V}\le -{\lambda }_{\min }\left(Q\right){\Vert z\Vert }_2^2+2{\Vert z^{T}PB\Delta B\left(z\right)\Vert }_2$$ (34)

If there exists a constant α with λ_min(Q) > α such that

$$2{\Vert z^T\hspace{0.17em}PB\Delta B\left(z\right)\Vert }_2\le \alpha {\Vert z\Vert }_2^2$$ (35)

Then (34) becomes

$$\dot{V}\le -\left[{\lambda }_{\min }\left(Q\right)-\alpha \right]{\Vert z\Vert }_2^2$$ (36)

Clearly, one can have $\dot{V}\left(x\right)<0$ , i.e. asymptotic stability of the uncertain system under feedback control, by the suitably chosen matrix Q.

Demonstrative Example

Cascaded analysis model

According to Figure 3 and (9), we set all parameters to construct the S-system as follows

$$\begin{array}{l}{\dot{z}}_1=2z_2^{0.5}{z}_3^{0.5}-2z_1^1,\\ {\dot{z}}_2=2z_4^{0.5}-2z_2^{0.6}{z}_3^{0.4},\\ {\dot{z}}_3=2z_5^{0.5}-2z_2^{0.4}{z}_3^{0.6},\\ z_4=0.5,\\ z_5=0.5\end{array}$$ (37)

The dynamic analysis was simulated by using the program: Power Law Analysis and Simulation (PLAS). The output concentration of the S-systems is obtained as z₁ = 0.701069. On the basis of the output concentration and (12), we have the parameters for the next stage model as follows

$$\begin{array}{l}{\dot{z}}_6=2z_7^{0.5}{z}_8^{0.5}-2z_6^1,\\ {\dot{z}}_7=1.68z_9^{0.5}-1.68z_7^{0.6}{z}_8^{0.4},\\ {\dot{z}}_8=2z_{10}^{0.5}-2z_7^{0.4}{z}_8^{0.6},\\ z_9=0.701,\\ z_{10}=0.5\end{array}$$ (38)

The output concentration of the S-systems is z₆ = 0.7696588.

Consider the complete S-system with its signal transduction pathways illustrated as in Figure 7. The S-system can be represented as

Figure 7.

Reference signal transduction network for the complete model.

$$\begin{array}{l}{\dot{\widehat{z}}}_1=2{\widehat{z}}_6^{0.5}-2{\widehat{z}}_1^{0.6}{\widehat{z}}_2^{0.4},\\ {\dot{\widehat{z}}}_2=2{\widehat{z}}_7^{0.5}-2{\widehat{z}}_1^{0.4}{\widehat{z}}_2^{0.6},\\ {\dot{\widehat{z}}}_3=2{\widehat{z}}_1^{0.5}{\widehat{z}}_2^{0.5}-2{\widehat{z}}_3^{0.6}{\widehat{z}}_4^{0.4}\\ {\dot{\widehat{z}}}_4=2{\widehat{z}}_8^{0.5}-2{\widehat{z}}_3^{0.4}{\widehat{z}}_4^{0.6},\\ {\dot{\widehat{z}}}_5=2{\widehat{z}}_3^{0.5}{\widehat{z}}_4^{0.5}-2{\widehat{z}}_5^1,\\ {\widehat{z}}_6={\widehat{z}}_7={\widehat{z}}_8=1\end{array}$$ (39)

The steady output concentration of the S-systems is ẑ₅ = 0.7071418. Compared the cascaded analysis model (38) with the complete model (39), the dynamic behavior of both models are quite similar as displayed in Figure 8, where the steady-state error between the two cases is less than 6%.

Figure 8.

Comparison of the output concentration. Dashed line denotes the output concentration of the cascaded analysis model, solid line denotes the output concentration of the reference system.

Stability

The reference S-system for the signal transduction network is modeled as

$$\begin{array}{l}\begin{array}{cc}{\dot{x}}_1=2x_2-1.2x_1^{0.5}{x}_3^{-1},& x_1\left(0\right)=2\end{array}\\ \begin{array}{cc}{\dot{x}}_2=2x_1^{0.1}{x}_3^{-1}{x}_4^{0.5}-2x_2,& x_2\left(0\right)=0.1\end{array}\\ x_3=0.5\\ x_4=1\end{array}$$ (40)

and the permissible ranges of the state variables are 1.3456 ≤ x₁ ≤ 3.5861 and 0.1 ≤ x₂ ≤ 2.2724.

Performing the second-order Taylor expansion for the S-system (40) around the steady state gives the linearized S-system:

$$\begin{array}{l}{\dot{x}}_1=-0.9505x_1+2x_2+7.3853+0.0442x_1^2,\\ {\dot{x}}_2=-2x_1-0.3631x_2+7.8019-0.0756x_1^2\end{array}$$ (41)

Consider (41). in the compact matrix form as

$$\left[\begin{array}{c}{\dot{x}}_1\\ {\dot{x}}_2\end{array}\right]=\left[\begin{array}{cc}-0.9505& 2\\ -2& -0.3631\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]+\left[\begin{array}{c}7.3853\\ 7.8019\end{array}\right]+\left[\begin{array}{c}0.0442x_1^2\\ -0.0378x_1^2\end{array}\right]$$ (42)

where the system eigenvalues are –0.6052 and –2.3453 and the residual errors satisfy

$$\begin{array}{l}\Vert \left[\begin{array}{c}0.0442x_1^2\\ -0.0378x_1^2\end{array}\right]\Vert \le \alpha {\Vert \left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]\Vert }_2\\ \Rightarrow 0.0038x_1^4\le {\alpha }^2\left(x_1^2+x_2^2\right)\end{array}$$ (43)

For the extreme case, one can find the permissible lower bound of α is 0.0827.

To discuss the robust stability of (42), we choose Q = I₂. According to the Lyapunov stability theory, there is a unique P = P^T > 0 to (16) as

$$P=\left[\begin{array}{cc}1.4852& -0.2058\\ -0.2058& 1.6204\end{array}\right]$$ (44)

From (17), the permissible range of α ensuring where stability is given by

$$\alpha <\frac{{\lambda }_{\min }\left(Q\right)}{{\lambda }_{\max }\left(P\right)}=0.6588$$ (45)

That is, the system would remain its stability provided that 0.0827 ≤ α ≤ 0.6588.

Control design

Consider the signal transduction network illustrated in Figure 9. Let one of the independent variable x₃ be the control variable and the dependent variables x₁ and x₂ be the state variables. In order to discern the control variable and state variables, we define the control variable as u. Then the S-system can be written as follows

$$\begin{array}{l}{\dot{x}}_1=2x_2-2.4x_1^{0.5},\\ {\dot{x}}_2=4x_1^{0.1}{x}_3-2x_2,\\ x_3=u\end{array}$$ (46)

Figure 9.

Reference system for the control design.

Now one can compute the linearization as the following steps. First, the system can be expressed in the control format as

$$\dot{x}=f\left(x\right)+l\left(x\right)u$$ (47)

where

$$\begin{array}{cc}f\left(x\right)=\left[\begin{array}{c}2x_2-2.4x_1^{0.5}\\ -2x_2\end{array}\right],& l\left(x\right)=\left[\begin{array}{c}0\\ 4x_1^{0.1}\end{array}\right]\end{array}$$ (48)

For the nominal system (47), a coordinate transformation exists as

$$\left[\begin{array}{c}{z}_1\\ z_2\end{array}\right]=\left[\begin{array}{c}{\psi }_1\left(x\right)\\ {\psi }_2\left(x\right)\end{array}\right]=\left[\begin{array}{c}h\left(x\right)\\ L_{f}h\left(x\right)\end{array}\right]$$ (49)

where h(x) = x₁, which satisfies $\frac{\partial h}{\partial x}l=0,\frac{\partial \left(L_{f}h\right)}{\partial x}l\ne 0$ . That is, setting

$$\left[\begin{array}{c}{z}_1\\ z_2\end{array}\right]=\left[\begin{array}{c}{x}_1\\ 2x_2-2.4x_1^{0.5}\end{array}\right]$$ (50)

transforms the original state-space representation into

$$\left[\begin{array}{c}{\dot{z}}_1\\ {\dot{z}}_2\end{array}\right]=\left[\begin{array}{cc}0& 1\\ 0& 0\end{array}\right]\left[\begin{array}{c}{z}_1\\ z_2\end{array}\right]+\left[\begin{array}{c}0\\ 1\end{array}\right]v$$ (51)

Next, select the new control input

$$v=-\left[\begin{array}{cc}{k}_1& k_2\end{array}\right]z$$ (52)

and substitute this into (51) to give

$$\dot{z}=A_{c}z$$ (53)

where

$$A_c=\left[\begin{array}{cc}0& 1\\ -k_1& -k_2\end{array}\right]$$ (54)

On the basis of (51), the control law for the original can be chosen as

$$u=\frac{v-L_f^{2}h\left(x\right)}{L_l{L}_{f}h\left(x\right)}$$ (55)

The S-system with control can then be expressed in the closed-loop configuration as

$$\begin{array}{l}{\dot{x}}_1=2x_2-2.4x_1^{0.5},\\ {\dot{x}}_2=4x_1^{0.1}u-2x_2,\\ u=\frac{-k_1{x}_1-k_2\left(2x_2-2.4x_1^{0.5}\right)+2.4x_1^{0.5}{x}_2+4x-2.88}{8x_1^{0.1}}\end{array}$$ (56)

Letting K = [1 100] yields the transient responses illustrated as in Figure 10. The result shows the effect of control input for signal transduction.

Figure 10.

Dynamic simulation for the case of K = [1 100].

Next, consider the system (47) with uncertainties and Δf^* = 0.1. On the basis of (25)–(30) and (48)–(54), one can derive a linearizd system as follows

$$\dot{z}=\left[\begin{array}{cc}0& 1\\ -k_1& -k_2\end{array}\right]z+\left[\begin{array}{c}0\\ 1\end{array}\right]\Delta B\left(z\right)$$ (57)

where ΔB(z) = (8x₁^0.1 Δf*) ○ ψ^–1 (z). Select K = [10 100] then (57) becomes

$$\dot{z}=\left[\begin{array}{cc}0& 1\\ -10& -100\end{array}\right]z+\left[\begin{array}{c}0\\ 0.8z_1^{0.1}\end{array}\right]$$ (58)

We now choose Q = I₂ and solve P from (33) as follows

$$P=\left[\begin{array}{cc}5.0550& 0.0500\\ 0.0500& 0.0055\end{array}\right]$$ (59)

From (35), we have

$$\frac{0.08z_1^{1.1}+0.0088z_1^{0.1}{z}_2}{z_1^2+z_2^2}\le {\alpha }_1$$ (60)

and from (36), we have

$$\alpha <{\lambda }_{\min }\left(Q\right)\Rightarrow {\alpha }_1<1$$ (61)

That is, if α₁ satisfies (60) and (61), then the system would be robustly stable.

Second, we choose Q = 0.5I₂ and solve P as

$$P=\left[\begin{array}{cc}2.52750& 0.02500\\ 0.02500& 0.00275\end{array}\right]$$ (62)

From (35), we have

$$\frac{0.04z_1^{1.1}+0.0044z_1^{0.1}{z}_2}{z_1^2+z_2^2}\le {\alpha }_2$$ (63)

and from (36), we have

$$\alpha <{\lambda }_{\min }\left(Q\right)\Rightarrow {\alpha }_2<0.5$$ (64)

We examine the value of Q to the permissible range of α. Let the permissible ranges of z₁ and z₂ are 0.0604 ≤ z₁ ≤ 0.9141 and –0.0836 ≤ z₂ ≤ 0, respectively. Considering the extreme case of (60), one can find

$$0.0852\le {\alpha }_1$$ (65)

According to (61) and (65), the permissible range of α₁ is then given by

$$0.0852\le {\alpha }_1<1$$ (66)

Similarly, the permissible range of α₂ is

$$0.0426\le {\alpha }_2<0.5$$ (67)

On the basis of (66) and (67), one can easily find that the permissible range of α is proportional to the magnitude of Q.

Control design for cascaded analysis model

Consider the complete S-system as follows and the signal transduction pathways is shown in Figure 11.

$$\begin{array}{l}{\dot{\tilde{x}}}_9=2{\tilde{x}}_{11}-2{\tilde{x}}_8^{1.01}{\tilde{x}}_9,\\ {\dot{\tilde{x}}}_8=2{\tilde{x}}_{10}-2{\tilde{x}}_8^{0.5}{\tilde{x}}_9^{0.5},\\ {\dot{\tilde{x}}}_6=2{\tilde{x}}_7-2{\tilde{x}}_8^{1.01}{\tilde{x}}_9,\\ {\dot{\tilde{x}}}_5=2{\tilde{x}}_8^{0.5}{\tilde{x}}_9^{0.5}-2{\tilde{x}}_5^{0.5}{\tilde{x}}_6^{0.5},\\ {\dot{\tilde{x}}}_3=2{\tilde{x}}_4-2{\tilde{x}}_2^{1.01}{\tilde{x}}_3,\\ {\dot{\tilde{x}}}_2=2{\tilde{x}}_5^{0.5}{\tilde{x}}_6^{0.5}-2{\tilde{x}}_2^{0.5}{x}_3^{0.5},\\ {\dot{\tilde{x}}}_1=2{\tilde{x}}_2^{0.5}{\tilde{x}}_3^{0.5}-2{\tilde{x}}_1,\\ {\tilde{x}}_4={\tilde{x}}_7={\tilde{x}}_{10}={\tilde{x}}_{11}=1\end{array}$$

Figure 11.

Reference S-system in control design for cascaded analysis model.

The S-system can be further expressed in a three-layered cascaded analysis model shown as in Figure 12. We construct the first layered S-system as follows

$$\begin{array}{l}{\dot{x}}_{13}=2x_{15}-2x_1^{1.01}{x}_2,\\ {\dot{x}}_{12}=0.5x_{14}-0.5x_{12}^{0.5}{x}_{13}^{0.5},\\ {\dot{x}}_{11}=2x_{12}^{0.5}{x}_{13}^{0.5}-2x_{11}\end{array}$$

Figure 12.

Three-layered cascaded analysis model.

The output concentration of the S-systems is x₁₁ = 1. On the basis of the output concentration, we proceed to construct the second layer as

$$\begin{array}{l}{\dot{x}}_8=2x_9-2x_7^{1.01}{x}_8,\\ {\dot{x}}_7=x_{10}-x_7^{0.5}{x}_8^{0.5},\\ {\dot{x}}_6=2x_7^{0.5}{x}_8^{0.5}-2x_6\end{array}$$

The output concentration of the S-systems is x₆ = 1. On the basis of this model, we construct the third layer as

$$\begin{array}{l}{\dot{x}}_1=2x_4-2x_1^{1.01}{x}_2,\hfill \\ {\dot{x}}_2=2x_5-2x_1^{0.5}{x}_2^{0.5},\hfill \\ {\dot{x}}_3=2x_1^{0.5}{x}_2^{0.5}-2x_3\hfill \end{array}$$

Now select the independent variable x₄ be the control variable, other independent variable be a constant and x₁, x₂, x₃ be the state variables. ¹¹ We set all parameters then the S-system form can be written as follows

$$\begin{array}{l}{\dot{x}}_1=u-2x_1^{1.01}{x}_2,\hfill \\ {\dot{x}}_2=1-2x_1^{0.5}{x}_2^{0.5},\hfill \\ {\dot{x}}_3=2x_1^{0.5}{x}_2^{0.5}-2x_3\hfill \end{array}$$ (68)

Now one can compute the linearized model. First, the system (68) is expressed in the control format as

$$\dot{x}=f\left(x\right)+l\left(x\right)u$$ (69)

where

$$f\left(x\right)=\left[\begin{array}{c}-2x_1{x}_2\\ 1-2x_1^{0.5}{x}_2^{0.5}\\ 2x_1^{0.5}{x}_2^{0.5}-2x_3\end{array}\right],l\left(x\right)=\left[\begin{array}{c}1\\ 0\\ 0\end{array}\right]$$

For the nominal system (69), a coordinate transformation exists as follows

$$\left[\begin{array}{c}{z}_1\\ z_2\\ z_3\end{array}\right]=\left[\begin{array}{c}{\psi }_1\left(x\right)\\ {\psi }_2\left(x\right)\\ {\psi }_3\left(x\right)\end{array}\right]=\left[\begin{array}{c}h\left(x\right)\\ L_{f}h\left(x\right)\\ L_f^{2}h\left(x\right)\end{array}\right]$$ (70)

where h(x) = x₂ + x₃, which stratifies $\begin{array}{ccc}\frac{\partial h}{\partial x}l=0,& \frac{\partial \left(L_{f}h\right)}{\partial x}l=0,& \frac{\partial \left(L_f^{2}h\right)}{\partial x}l\ne 0\end{array}$ . That is, setting

$$\left[\begin{array}{c}{z}_1\\ z_2\\ z_3\end{array}\right]=\left[\begin{array}{c}{x}_2+x_3\\ 1-2x_3\\ -4x_1^{0.5}{x}_2^{0.5}+4x_3\end{array}\right]$$ (71)

transforms the original state-space representation into

$$\left[\begin{array}{c}{\dot{z}}_1\\ {\dot{z}}_2\\ {\dot{z}}_3\end{array}\right]=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ 0& 0& 0\end{array}\right]\left[\begin{array}{c}{z}_1\\ z_2\\ z_3\end{array}\right]+\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]v$$ (72)

Next, selecting the new control input v = –Kz and substitute this into (72) gives

$$\dot{z}=A_{c}z$$ (73)

where

$$A=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ -k_1& -k_2& -k_3\end{array}\right]$$

On the basis of (72) the control law for the original can be chosen as

$$u=\frac{v-L_f^{3}h\left(x\right)}{L_l{L}_f^{2}h\left(x\right)}$$ (74)

The S-system with control can then be expressed in the closed-loop configuration as

$$\begin{array}{l}{\dot{x}}_1=u-2x_1^{1.01}{x}_2,\\ {\dot{x}}_2=1-2x_1^{0.5}{x}_2^{0.5},\\ {\dot{x}}_3=2x_1^{0.5}{x}_2^{0.5}-2x_3,\\ u=\frac{-k_1\left(x_2+x_3\right)-k_2\left(1-2x_3\right)-k_3\left(-4x_1^{0.5}{x}_2^{0.5}+4x_3\right)-4x_1^{0.5}{x}_2^{1.5}+2x_1^{0.5}{x}_2^{-0.5}-4x_1-8x_1^{0.5}{x}_2^{0.5}+8x_3}{-2x_1^{-0.5}{x}_2^{0.5}}\end{array}$$ (75)

Letting K = [1 10 40] yields the transient responses illustrated as in Figure 13.

Figure 13.

Dynamic simulation for the case of K = [1 10 40].

For the demonstrative purpose, consider the situation of Δf^* = 0.1. On the basis of (25)–(30) and (73), one can derive the linearized system as follows

$$\dot{z}=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ -k_1& -k_2& -k_3\end{array}\right]z+\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]\Delta B\left(z\right)$$ (76)

where ΔB(z) = (–2x₁^–0.5x₂^0.5Δf*) ○ ψ^–1(z) Select K = [1 10 20], then (76) becomes

$$\dot{z}=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ -1& -10& -20\end{array}\right]z+\left[\begin{array}{c}0\\ 0\\ 0.4\left(\frac{2z_1+z_2-1}{z_3+2z_2-2}\right)\end{array}\right]$$

We choose Q = I₃ and compute P from solving (33) as

$$P=\left[\begin{array}{ccc}6.0578& 10.0779& 0.5000\\ 10.0779& 21.4347& 1.0578\\ 0.5000& 1.0578& 0.0779\end{array}\right]$$

From (35), we have

$$\frac{2\left(0.2z_1\left(\frac{2z_1+z_2-1}{z_3+2z_2-2}\right)+0.42312z_2\left(\frac{2z_1+z_2-1}{z_3+2z_2-2}\right)+0.03116z_3\left(\frac{2z_1+z_2-1}{z_3+2z_2-2}\right)\right)}{z_1^2+z_2^2+z_3^2}\le \alpha$$

and from (36), we have

$$\alpha <{\lambda }_{\min }\left(Q\right)\Rightarrow \alpha \le 1$$

That means the system under control would be robustly stable.

Conclusion

This paper proposes a method for constructing the dynamic model of signal transduction networks and a primary control design has been proposed for the system. A cascaded analysis model for constructing the signal transduction network model has been proposed. The advantage of cascaded analysis model is that the model preserved the major dynamic feature of the S-systems with a simplified model while avoiding much computation burden. The stability condition for the linearized S-system has been derived. A method for controlling the steady state to the desired level using the technique of feedback linearization has also been attempted. The development of this issue is undergoing theoretical investigation and experimental verification.