Optimization-based Inference for Temporally Evolving Networks with Applications in Biology

2.1. Dynamic graph (linear time varying system)

The state X_k evolves along with time and constitutes the following linear-time varying system: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}X_k = G_k X_{k - 1} \tag{4}\end{align*} \end{document}

where G_k = g(GR,X_k|X_k−1) is a function of both the connectivity map and time series data. Note that G_k describes how the edge activities evolve over time. For example, for given connectivity map GR, we allow the change of strength of connection to drive our dynamic model consistent with biological system or experimental data. At each time step, only a few edges may evolve based on the relationship between X_k and X_k−1. If all the interactions between each component are properly identified, we can reconstruct the map G_k in terms of the connectivity and strength. For instance, G_k(i, j) = 0.5 represents that node j activates node i with strength 0.5. The strength might be related to the reaction rate and the concentration of other species, demonstrated by the Jacobian of a mass-action kinetics model.

The goal of system identification of biological systems is to infer each G_k for 1 ≤ k ≤ N consistent with both a biological data set χ and a priori information GR. In general, a gene regulation network (GRN) has the following characteristics (Marc et al., 2010):

1) Directionality: Regulatory control is directed from regulators to regulated genes.

Since GRNs have a sparse structure with combinatorial control, we should reconstruct the sparsest graph consistent with experimental observations. We can construct an optimization problem as follows: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} & \min_{G_k} \mid \mid {X}_{k} - G_k {X} _{k - 1} \mid \mid + \gamma \mid \mid {A_k} \mid \mid \\ & subject \ to \ G_k = g ( GR , X_k \mid X_{k - 1} ) \tag{5} \end{align*} \end{document}

where the second term in the cost function penalizes the cost of adding edges in order to avoid heavy combinatoric computation, A_k is defined as follows: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}A_k = G_k \cap ( G_k \cap G_{k - 1} ) ^c = G_k - G_{k - 1} \tag{6}\end{align*} \end{document}

and γ is a positive constant. Therefore, A_k enforces the network to be sparse, and thus, the cost function represents a trade-off between reconstruction error and sparsity. Here we define the function g as shown below: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}G_k = g ( GR , X_k \mid X_{k - 1} ) = GR \oplus MAP = {\rm Proj}_{MAP} GR \tag{7}\end{align*} \end{document}

where ⨁ is defined as a projection operator onto \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$MAP \in { \bf R}^{n_x \times n_x}$$ \end{document} , whose i-th column is a column vector, the components of which are all one if X_k−1(i) is active, which means the state of the i-th element is over the threshold. On the other hand, if X_k−1(i) is nonactive, then the i-th column of MAP is a zero-column vector. Therefore, this projection gives us all possible candidate edges based on both X_k−1 and GR. For example, if x_i at the (k − 1)th step is active, then the i-th column of GR contains the candidate edges. On the other hand, if x_i at (k − 1)th step is not active, we cannot use the i-th column of GR as candidate edges. By using Equation (7), our method generates a sparse network representation without any Lasso-type regressions in Equation (5).

If we implement the optimization problem for every single step as shown in Equation (5), the penalty term for sparsity does not play the role of generating a sparse network but instead uses all possible edges. In other words, distributing signal to all possible nodes (dense network) gives us a lower cost than distributing signal to only a few nodes (sparse network) in our formulation (Eq. 5). We can think about this by considering dynamic programming in optimization. The main idea behind dynamic programming is that, to solve a given problem, we need to solve different parts of the problem (subproblems), then combine the solutions of the subproblem to reach an overall solution in a recursive manner. Similarly, in order to find the sparsest smoothly evolving graph, we need to have a certain connection between every subproblem. For example, when we consider the overall time horizon in Equation (8), the penalty term for sparsity can play a key role in generating a sparse graph structure by connecting the discrete time dynamics at each time step with those at different time steps. Then, we can construct a convex optimization problem for the proposed identification problem as shown below: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} & \min_{G_1 , \ldots , G_{N}} \sum_{k = 1}^{N} \mid \mid {X}_{k} - G_k {X}_{k - 1} \mid \mid _2 + \gamma \left\{\sum_{k = 2}^{N} \mid \mid G_{k} - G_{k - 1} \mid \mid _F + \mid \mid G_1 \mid \mid _F + \mid \mid G_N \mid \mid _F \right\} \\ & subject \ to \ G_k = g ( GR , X_k \mid X_{k - 1} ) \tag{8}\end{align*} \end{document}

Note that the first term of Equation (8), the summation of ||X_k − G_kX_k−1||, forces the minimization of the reconstruction error for a given dynamical network at time k for 1 ≤ k ≤ N. Also, the second term, the summation of ||G_k − G_k−1||_F, plays the role of realizing a smooth evolution and minimizes the change in network evolution. Finally, with the penalty term ||G₁|| _F  + ||G_N||_F, which acts as a boundary constraint, we can find the sparsest dynamic graph. We can also arrange and reformulate Equation (8) as follows: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \begin{split} & \min_{ {\cal G}} \mid \mid {\cal X} - {\cal GX} \mid \mid _2 + \gamma \mid \mid ( { \cal G}^T - { \cal G} ) \times W \mid \mid _F \\ & subject \ to \ given \ { \cal X} , W \\ & G_k^{act} \quad \ge 0 , G_k^{inhib} \ \le 0 , G_k^{others} = 0\end{split} \tag{9}\end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}{ \rm where} \ { \cal G} = \left[ \begin{matrix}O_{n_x} & G_N & O_{n_x} & \cdots & O_{n_x} \\ O_{n_x} & O_{n_x} & G_{N - 1} & \cdots & O_{n_x} \\ \cdots & \cdots & \cdots & \cdots & \cdots \\ O_{n_x} & O_{n_x} & \cdots & O_{n_x} & G_1 \\ O_{n_x} & O_{n_x} & \cdots & O_{n_x} & I_{n_x}\end{matrix} \right] , { \cal X} = \left[ \begin{matrix}X_N \\ X_{N - 1} \\ \cdots \\ X_1 \\ X_0\end{matrix} \right] { \rm and} \ W = \left[ \begin{matrix}I_{n_x} & O_{n_x} \\ O_{n_x} & I_{n_x} \\ I_{n_x} & O_{n_x} \\ O_{n_x} & I_{n_x} \\ \cdots & \cdots \\ \cdots & \cdots\end{matrix} \right] { \rm }\end{align*} \end{document}

where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \cal X} \in { \bf R} ^{n_x \cdot ( N + 1 ) \times 1} , { \cal G} \in { \bf R} ^{n_x \cdot ( N + 1 ) \times n_x \cdot ( N + 1 ) }$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$W \in { \bf R} ^{n_x \cdot ( N + 1 ) \times 2 \cdot n_x}$$ \end{document} for 1 ≤ k ≥ N. Note that the first term of the cost function in Equation (9) is a reconstruction error cost, and the second term plays the role of connecting each discrete system with another and realizing a smooth evolution of the network by selecting effective edges with inequality constraints.

2.2. Static graph (linear time invariant system)

If we assume that the graph model does not evolve with time (static graph G_k = G), such as with a linear-time invariant system (Han et al., 2007), we can modify the structure of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \cal G}$$ \end{document} and constraints as shown below for a fixed pattern graph: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \begin{split} & { \cal G} = \left[\begin{matrix}O_{n_x} & G & O_{n_x} & \cdots & O_{n_x} \\ O_{n_x} & O_{n_x} & G & \cdots & O_{n_x} \\ \cdots & \cdots & \cdots & \cdots & \cdots \\ O_{n_x} & O_{n_x} & O_{n_x} & O_{n_x} & G \\ O_{n_x} & O_{n_x} & O_{n_x} & O_{n_x} & I_{n_x}\end{matrix} \right] \\ & {G^{act} \ge 0 , G^{inhib} \le 0 , G^{others} = 0} \end{split} \tag{10}\end{align*} \end{document}

where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$G = \bar{g} ( GR )$$ \end{document} does not depend on time [compared with G_k = g(GR,X_k|X_k−1) for a linear time varying system]. Note that for a fixed graph structure, the optimal solution represents the average connectivity map (Han et al., 2007).

2.3. Dynamic and static graph

We can compare the dynamic graph and static graph method: the main difference in cost function from dynamic and static graph is the penalty for sparsity as follows: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\mid \mid ( { \cal G}^T - { \cal G} ) W ) \mid \mid _{F , dynamic} = \left[ \begin{matrix}O_{n_x} & - G_N \\ \Delta G_{N - 1} & O_{n_x} \\ O_{n_x} & \Delta G_{N - 2} \\ \cdots & \cdots \\ O_{n_x} & \Delta G_{1} \\ G_1 & O_{n_x}\end{matrix} \right] \tag{11}\end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\mid \mid ( { \cal G}^T - { \cal G} ) W ) \mid \mid _{F , static} = \left[ \begin{matrix}O_{n_x} & - G \\ O_{n_x} & O_{n_x} \\ O_{n_x} & O_{n_x} \\ \cdots & \cdots \\ O_{n_x} & O_{n_x} \\ G & O_{n_x}\end{matrix} \right] \tag{12}\end{align*} \end{document}

where ΔG_k = G_k+1 − G_k. Also, if we modify the constraint for a dynamic graph similar to the static graph approach, the dynamic graph approach gives us a lower cost than the static graph approach because the structural constraint restricts the degrees of freedom in choosing edges: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \begin{split} \min_{ \cal G} \mid \mid { \cal X} - {\cal GX} \mid \mid _2 + \gamma \mid \mid ( { \cal G}^T - { \cal G} ) \times W \mid \mid _F \\ \le \min_{ \bar { \cal G}} \mid \mid { \cal X} - { \bar { \cal G}} { \cal X} \mid \mid _2 + \gamma \mid \mid ( { \bar { \cal G}}^T - { \bar { \cal G}} ) \times W \mid \mid _F\end{split} \tag{13}\end{align*} \end{document}

where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \cal G}$$ \end{document} represents the optimal solution of the dynamic graph approach, and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\bar{{ \cal G}}$$ \end{document} represents the optimal solution of the static graph approach.

2.4. Inhibition edges

Based on our formulation of the optimization problem, we can find the optimal solution that satisfies a trade-off between representation of data (dynamics) and sparsity and smooth evolution. However, the optimal solution does not include any inhibition edges (⊣) because it is not necessary according to our optimization problem as shown in Figure 1. For example, if X is active (1) and Y is not active (0), then there might be two possible cases: X inhibits Y (X ⊣ Y connected, Figure 1 [left]) or no connection between X and Y (Figure 1 [right]). Having no connection would give the lower cost. However, we can handle inhibition edges using Boolean logic as an algebraic constraint as shown below: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}{ \rm Y} = { \rm not} \ { \bar {\rm Y}} ( = \sim {\bar {\rm Y}} ) \tag{14}\end{align*} \end{document}

OPEN IN VIEWER

FIG. 1.

Possible cases for inhibition edge (dash end), where X and Y represent different genes or states; 1 represents activated state and 0 represents deactivated state: (a) inhibition reaction is triggered and (b) inhibition reaction has not occurred.

Also, we extend this algebraic constraint to a normalized state as shown below: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}{ \rm Y} + { \bar {\rm Y}} = 1 \tag{15}\end{align*} \end{document}

Consider the simple case shown in Figure 2, in which state X inhibits state Y. Using an algebraic constraint (Eq. 14), we can change the inhibition edge to an activation edge with the new state \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\bar{\rm Y}$$ \end{document} as shown below: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}{ \rm X} \dashv { \rm Y} \Leftrightarrow { \rm X} \dashv ( \sim \bar { \rm Y} ) \Leftrightarrow { \rm X} \rightarrow { \bar {\rm Y}} - { \rm Y} \tag{16}\end{align*} \end{document}

OPEN IN VIEWER

FIG. 2.

(a) Inhibition edge ( \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${\rm Y}$$ \end{document} inhibits \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${\rm Y}$$ \end{document} ), and (b) modified edge ( \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${\rm X}$$ \end{document} activates \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\bar{\rm Y}$$ \end{document} and the relation between \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${\rm Y}$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${\rm Y}$$ \end{document} is defined by Boolean logic or algebraic constraint).

Hence, we extend states if there are inhibition input edges and introduce a diagonal weighting matrix \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \cal M}$$ \end{document} , which makes all species have the same penalty as shown below: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\min_{ \cal G} \mid \mid { \cal M} \times ( { \bar { \cal X}} - { \cal G}{ \bar { \cal X}} ) \mid \mid _2 + \gamma \mid \mid ( { \cal G}^T - { \cal G} ) \times W \mid \mid _F \tag{17}\end{align*} \end{document}

where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\tilde{ \cal X}$$ \end{document} represents extended states and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$ { \cal M } _ { ii } = \ { 1 , \frac { 1 } { \sqrt 2 } \ } $$ \end{document} . If there exists \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$x , \bar x$$ \end{document} for a specific state, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$ { \cal M } _ { ii } = \frac { 1 } { \sqrt { 2 } } $$ \end{document} , and otherwise, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \cal M}_{ii} = 1$$ \end{document} .

2.5. Numerical examples

In this section, we consider simple examples to illustrate the proposed inference scheme.

2.5.1. Simple gene network

We first consider a toy example composed of four genes. The a priori information and the snapshot of gene expression are shown in Figure 3. Here, we do not consider state extension for inhibition edges, which means the optimal solution does not include any inhibition edges. By varying the parameter γ, we can sweep out the optimal trade-off curve between the reconstruction error and the sparsity of a solution as shown in Figure 4. We can choose the optimal parameter γ* by the graphical representation: the extreme point γ* on the trade-off between the sparsity and the reconstruction error. Once we fix the parameter γ*, we solve the constrained convex optimization problem (Eq. 29) using a MATLAB-based modeling system for convex optimization (CVX) (Boyd and Vandenberghe, 2004). Figure 5 shows the dynamics of the connectivity graph. We can capture the temporal graph not only in terms of connection but also by strength of the edge. From the optimal graph representations, we could extract how the signaling pathway evolves over time with a systems point of view. Also, we can compare the two approaches: dynamic and static graph approach with average of dynamic graph.

OPEN IN VIEWER

FIG. 3.

(a) A priori connectivity map, where the arrows indicate activation and blunted lines denote inhibition. (b) Snapshots of gene expression from time k = 1 to k = 4 (red or 1: activated states; green or 0: deactivated states).

OPEN IN VIEWER

FIG. 4.

Trade-off curve between the model fitting ( \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\mid \mid { \cal X} - { \cal GX} \mid \mid _2$$ \end{document} ) and the sparsity ( \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\mid \mid ( { \cal G}^T - { \cal G} ) \times W \mid \mid _F$$ \end{document} ) with turning parameter γ.

OPEN IN VIEWER

FIG. 5.

The optimal solution for the example in Section 2.5.1.: the magnitude of each edge represents strength of connection. Also, we compare the results from dynamic and static graph approaches with average of dynamic graph.

2.5.2. Simple gene network with different structure

Here, we add one edge that connects vertex 3 to vertex 2 as shown in Figure 6 and solve the optimization problem again. We can see the difference of the strength of edge 12 (e₁₂) compared with Section 2.5.1 above. Basically, for Section 2.5.1, the optimal graph shows the robust pathway distributing power evenly (e₁₂ and e₁₃ in Figure 5), because both pathways are effective equally. However, for Section 2.5.2, an additional pathway (e₃₂) changes the topology of the graph, which makes the optimal graph choose the more effective or economical path (e₁₃ − e₃₂) in Figure 7. In other words, the optimal solution shows that the strength of e₁₂ decreases because there exists a more effective pathway (e₁₃ − e₃₂).

OPEN IN VIEWER

FIG. 6.

A priori connectivity map with an additional edge (e₃₂), which connects from node 3 to node 2.

OPEN IN VIEWER

FIG 7.

An optimal solution for the example in Section 2.5.2.: an additional pathway (e₃₂) changes the topology of the graph, which makes the optimal graph choose a more effective path (e₁₃ − e₃₂) rather than e₁₂ (for the example in Section 2.5.1., the optimal graph shows distributing power evenly, e₁₂ = e₁₃).

2.6. Extension for the continuous dynamics case

In this article, we infer discrete time dynamics from a set of data with a priori graph structure information and focus on how to find the most reasonable model among many possible configurations. We can extend this idea for the continuous time case. The identification problem has led to a linear quadratic (LQ) optimal control problem with two main penalty functions by which we can match the experimental data with a sparse representation using a priori information of structure (described in the Appendix).

3.1. p53 Signal pathway

Aswani et al. (2009) proposed a graph-theoretic topological control applied to the p53 signaling pathway. We apply our approach to understand how the controller affects the biological pathway and capture the evolution of the signaling pathway. We define X = [x₁, x₂, x₃, x₄]^T = [x_MDM, x_p53, x_cyclinG, x_c], where x_c is a virtual state that represents the proposed control scheme (actually removing the edge in Aswani et al. 2009): \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \begin{split} &\rm {x}_{ c} = { x}_{{ MDM}} \ \ { if \ controller \ Off} \\ &\rm { x}_{ c} = 0 \ \qquad { if \ controller \ On}\end{split} \tag{18}\end{align*} \end{document}

Hence, by introducing this virtual state, we have an abstract model of abnormal p53 signaling pathway with controller in Figure 8 (right). Also, we can define GR as follows, based on Figure 9 with extending states, including the state extension due to incorporating the inhibition edges: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}GR = \overbrace{ \left[ \begin{matrix}0 & 0 & e_{21} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & e_{3 \bar1} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & e_{4 \bar2} \\ 0 & 0 & e_{23} & 0 & 0 & 0 \\ e_{14} & 0 & 0 & 0 & 0 & 0\end{matrix} \right] } ^{ \begin{matrix}x_1 & \bar x_1 & x_2 & \bar x_2 & x_3 & x_4\end{matrix}} = \left[ \begin{matrix}0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0\end{matrix} \right] \tag{19}\end{align*} \end{document}

OPEN IN VIEWER

FIG. 8.

(a) An abnormal p53 pathway in Figure 3c in Aswani et al., 2009, and (b) the abstract model that includes the effect of controller.

OPEN IN VIEWER

FIG. 9.

(a,b) Normalized time course plots for the abnormal p53 pathway with controller in Figure 4c in Aswani et al., 2009, and (c,d) dynamic evolution of each edge of abnormal p53 pathway with the controller: the p53 regulates MDM2 similar to the normal p53 pathway by increasing strength of inhibition edge ([p53]-[cyclin G]-[MDM2]).

Here, we normalize the data then apply our algorithm. We can capture the dynamic evolution of the graph in Figure 9. The controller causes the p53 concentrations to increase to higher levels by regulation edge from the murine double minute, an important negative regulator of the p53 tumor suppressor (MDM2) to p53 (also known as protein 53 or tumor protein 53) and causes increased strength of inhibition edge ([p53]-[cyclin G]-[MDM2]). In other words, p53 regulates MDM2 similar to the normal p53 pathway (Aswani et al., 2009). If the controller is not applied again, the strength of edge ([p53]-[cycle G]-[MDM2]) decreases, and the strength of activation edge [p53]-[Controller]-[MDM2] increases. This causes MDM concentrations to increase to higher levels which cause regulation p53 by inhibition edge ([MDM2]-[Controller]-[p53]) similar to the abnormal p53 pathway (Aswani et al., 2009).

We can also apply our algorithm for the normal p53 pathway in order to compare with the abnormal p53 pathway with controller. In a normal p53 pathway, we can consider all possible combinations of both Ras and L26 as two input signals. Here, the basic assumption is that the inhibition reaction is stronger than the activation reaction. Then, we find that the p19 (alternate reading frame of the INK4a/ARF locus (ARF)) mainly regulates MDM2, and it cannot affect MDM2 from p53 through p19 ARF as shown in Figure 10. Hence, we can use the same abstract model in Figure 8 for a normal p53 signaling pathway. The optimal solution shows that the normal p53 cell uses mainly inhibition edges from p53 to MDM2 through cyclin G, which means p53 regulates MDM2, as shown in Figure 11. Therefore, the controller drives the abnormal p53 cell to the normal p53 cell by removing the inhibition edge from MDM2 to p53 as Aswani et al. (2009) proposed.

OPEN IN VIEWER

FIG. 10.

Possible cases for a normal p53 signaling pathway with different combinations of both [Ras] and [L26], where H represents an active state and L represents a nonactive state (Aswani et al., 2009): the p19 ARF mainly regulates MDM2 and it cannot affect MDM2 from p53 through p19 ARF.

OPEN IN VIEWER

FIG. 11.

(a,b) Normalized time course plots for the normal p53 pathway in Figure 4a (Aswani et al., 2009), and (c,d) the dynamic evolution of each edge of the normal p53 pathway: the edge activity shows that the normal p53 cell uses mainly inhibition edges from p53 to MDM2 through cyclin G.

3.2. HER2 overexpressed breast cancer

We are interested in HER2 overexpressed breast cancer, which represents 20–30% of breast cancers. The experimental studies were done for investigating the effects of tyrosine kinase inhibitors (TKIs) on the BT474 and SKBr3 cell lines (Sergina et al., 2007). In this work, short-term effects and long-term effects of applying gefitinib (a TKI) to those cell lines were studied, and important effects of how the cancer cells overcome or escape from the inhibitory effects of TKIs were discovered. The authors in Sergina et al. (2007) showed that HER3 is recruited from the cytoplasm to the cell membrane by vesicular trafficking to increase the triggering signal in order to escape from HER2 inhibition. Also, they tested the effects of vesicular trafficking: when vesicular trafficking was stopped, phospho-HER3 and phospho-Akt did not survive the inhibition of HER2.

We suspect there might be short-term and long-term topological changes because the TKI can inhibit and regulate downstream effectively in the short term, but it cannot regulate for the long term. Therefore, we hypothesize that during the short term, there might be positive negative (PN) feedback (Kim et al., 2007), so the TKI inhibits HER3 effectively. However, for long-term behavior, even a small triggering signal could amplify the phospho-Akt signal, because of positive positive (PP) feedback, which is similar to vesicular trafficking. On the other hand, if the topology does not change, TKI should be able to regulate downstream over the long term even though HER3 is recruited by vesicular trafficking.

We define the a priori map from biological information (Sergina et al., 2007; Amin et al., 2010; Itani et al., 2009), where we include a nucleus model to capture this possible topology change. The behaviors of the nucleus are not yet understood; however, we abstract it with the switch as shown in Figure 12. Basically, there is a fail-safe mechanism, HER2-HER3 signaling, which is buffered so that it is protected against an inhibition of HER2 catalytic activity, and it is driven by the negative regulation of HER3 by Akt (Amin et al., 2010). Also, there is a compensatory mechanism by cross-talk between MAPK and Akt that results in robust activation of this buffering. However, the compensatory buffering prevents apoptotic tumor cell death from occurring as a result of the combined loss of MAPK and Akt signaling (Amin et al., 2010). For example, once a signal is triggered and either MAPK or Akt is high, then the nucleus stays active so MAPK and/or Akt are trying to keep the compensatory buffering. However, once both MAPK and Akt are downregulated, the nucleus is deactivated for all time.

OPEN IN VIEWER

FIG. 12.

Bifan motif of nucleus, which is two-layered graph with edges from nodes in top- to bottom-layer: there is a fail-safe mechanism, HER2-HER3 signaling, which is buffered so that it is protected against an inhibition of HER2 catalytic activity, and a compensatory mechanism by cross-talk between MAPK and Akt, which results in robust activation of buffering.

We apply the proposed optimization technique and the result is shown in Figure 13 and 14. Here, we use the generated data based on biological experimental data (western blot (Sergina et al., 2007; Amin et al., 2010). By applying the proposed method, we find that there are three main phases: before TKI is introduced (triggering network), right after TKI is introduced (short-term), and long-term behavior after TKI is introduced. We can capture the topology change of the biological network: for the initial stage (Fig. 14a), the signal is triggered and propagated along activation edges. After TKI is introduced, downstream components such as phospho-HER3, PI3K, Akt, and MAPK are regulated because TKI inhibits and regulates downstream components. Moreover, the biological network shows PN feedback, which effectively modulates signal responses. Finally, for long-term behavior, even if a small triggering signal is introduced (because of TKI inhibition, step 17–step 20), the downstream components are not regulated but are activated because the biological network evolves to PP feedback, which induces a slower but amplified signal response and enhances bi-stability.

OPEN IN VIEWER

FIG. 13.

The upper two panels show the normalized biological data and the assumed nucleus level. The lower panels show the strength of the downstream edges. For example, the edge connecting HER23 to MAPK (middle panel) is activated from step 4 to step 9 but deactivated from step 9 to 18.

OPEN IN VIEWER

FIG. 14.

(a) Signal is triggered and propagated along activation edges; (b) after TKI is introduced (short term), downstream components such as phospo- HER3, PI3K, Akt, and MAPK are regulated because TKI inhibits and regulates downstream components (positive negative feedback); and (c) for long-term behavior, even though a small triggering signal is introduced, the downstream components are not regulated but are activated because the biological network evolves to positive positive feedback (gray: not triggered edge; red: activation edge; blue: inhibition edge; light red/blue: deactivated edges after once activated).

5.1. Problem statement

We define a state vector \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$x ( t ) = [ x_1 ( t ) , \ldots , x_n ( t ) ] ^T$$ \end{document} , the components of which represent concentrations of proteins or states in a biological network. It is assumed that the state of the network evolves over time and this evolution of state x(t) can usually be modeled using an ODE: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\dot{x} = f ( x , p ) \tag{20}\end{align*} \end{document}

where p is a parameter set. As we mentioned in this article, many studies in systems biology impose a structure on f (ċ), such as mass-action kinetics or Hill function dynamics, and identify parameters using least-square criteria. However, we are interested in finding the discrete time-varying influence map that can be formulated as a discrete time-varying linear system. Here, we basically extend this idea for the continuous case. The nonlinear dynamic system (Eq. 20) can be approximated by a time-varying linear system based on forming the Jacobian around steady states as shown below: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\left[\begin {matrix} \delta \dot {x} _1 (t) \\ \delta \dot {x} _2 (t) \\ \ldots \\ \delta \dot {x} _n (t) \end {matrix} \right] = \left[ \begin {matrix} \frac {\partial f_1} {\partial x_1} & \frac {\partial f_1}{\partial x_2} & \ldots & \frac {\partial f_1} {\partial x_n} \\ \ldots & \ldots & \ldots & \ldots \\ \ldots & \ldots & \ldots & \dots \\ \frac {\partial f_n} {\partial x_1} & \frac {\partial f_n} {\partial x_2} & \ldots & \frac { \partial f_n } { \partial x_n } \end {matrix} \right] \left[\begin {matrix} \delta x_1 \\ \delta x_2 \\ \ldots \\ \delta x_n\end {matrix} \right] = G ( t ) \delta x ( t ) \tag {21} \end{align*} \end{document}

where we assume there is no parameter variation (δp = 0). A system in the form of Equation (21) can be considered as a temporally evolving weighted directed graph. Then, G(t) is a time-varying adjacency matrix, or influence matrix, of dimension n × n, which describes the temporal evolution of the edges with strength change. In general, G(t) is a sparse matrix (Marc et al., 2010): \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}G_ { i , j } ( t ) = \frac { \partial f_i } { \partial x_j } = \begin{cases} \neq 0 \hbox { if node } j \hbox { can affect node } i \hbox { directly } \\ \\ = 0 \hbox { otherwise } \end{cases} \tag { 22 } \end{align*} \end{document}

where G_i,j(t) is nonzero if there exists a direct connection between node j (input node) and node i (output node). Otherwise, G_i,j(t) is zero.

Definition. Let G(t) be a time-varying adjacency matrix that represents a dynamic graph with n nodes and k edges, where k is the number of candidate edges from a priori information. The component of G(t), denoted e(t) = comp(G(t)), is a k × 1 vector whose elements are the nonzero entries of G(t).

Example. Consider the dynamic graph shown in above Figure A1. Following the conventions introduced above, the corresponding adjacency matrix G(t) has the form: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}G ( t ) = \left[ \begin{matrix}0 & e_{21} & 0 & e_{41} \\ 0 & 0 & 0 & 0 \\ e_{13} & 0 & 0 & e_{43} \\ 0 & e_{24} & 0 & 0\end{matrix} \right] \tag{23}\end{align*} \end{document}

OPEN IN VIEWER

FIG. A1.

A simple graph model.

Its component e(t) is constructed by extracting the nonzero elements from each column, which produces the vector: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}e ( t ) = [ e_{13} , e_{21} , e_{24} , e_{41} , e_{43} ] ^T = : [ e_1 , e_2 , e_3 , e_4 , e_5 ] ^T \tag{24}\end{align*} \end{document}

Using e(t), we can reformulate Equation (21) as follows (for example, n = 4, k = 5): \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \left[ \begin{matrix} \dot{x}_1 \\ \dot{x}_2 \\ \dot{x}_3 \\ \dot{x}_4\end{matrix} \right] _{n \times 1} & = \left[ \begin{matrix}0 & e_{21} & 0 & e_{41} \\ 0 & 0 & 0 & 0 \\ e_{13} & 0 & 0 & e_{43} \\ 0 & e_{24} & 0 & 0\end{matrix} \right] _{n \times n} \left[ \begin{matrix}x_1 \\ x_2 \\ x_3 \\ x_4\end{matrix} \right] _{n \times 1} = \left[ \begin{matrix}0 & x_2 & 0 & x_4 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ x_1 & 0 & 0 & 0 & x_4 \\ 0 & 0 & x_2 & 0 & 0\end{matrix} \right] _{n \times k} \left[ \begin{matrix}e_{13} ( t ) \\ e_{21} ( t ) \\ e_{24} ( t ) \\ e_{41} ( t ) \\ e_{43} ( t ) \end{matrix} \right] _{k \times 1} \\ & = A ( x ) e ( t ) \\ & = e_1 ( t ) \left[ \begin{matrix} 0 \\ 0 \\ x_1 \\ 0\end{matrix} \right] + e_2 ( t ) \left[ \begin{matrix}x_2 \\ 0 \\ 0 \\ 0\end{matrix} \right] + e_3 ( t ) \left[ \begin{matrix}0 \\ 0 \\ 0 \\ x_2\end{matrix} \right] + e_4 ( t ) \left[ \begin{matrix}x_4 \\ 0 \\ 0 \\ 0\end{matrix} \right] + e_5 ( t ) \left[ \begin{matrix}0 \\ 0 \\ x_4 \\ 0\end{matrix} \right] \tag{25} \end{align*} \end{document}

where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$A ( x ) \in R^{n \times k}$$ \end{document} is a linear function of x, which can be constructed from a priori information, representing possible influence modes of biological networks. For example, the first mode, [0 0 x₁ 0] ^T in Equation (25), shows that node 1 activates node 3 (i.e., \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$x_1 \rightarrow \dot{x}_3 \rightarrow x_3$$ \end{document} ). Also, each e_i(t) represents a time-varying coefficient or an activity of i-th mode \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\in R^{n \times 1}$$ \end{document} . Therefore, we can assign the network topology by adding edges, for example, if there are suspicious interactions among nodes.

In order to find all e_i(t), which drive our system dynamics with influence modes, we formulate a linear quadratic (LQ) optimal control problem, as a regulation problem (x(t) → x_d(t)) with control inputs both e(t) and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\dot{e} ( t )$$ \end{document} . Once we solve the LQ problem, e_i(t)* shows the optimal activity or sequence of each mode over time, which drives our dynamic system to match biological data.

5.2. Time varying linear system

In order to formulate the LQ optimal control problem, we define the controlled system as follows: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \frac { dx ( t ) } { dt } = A ( x ) e ( t ) \tag {26} \end{align*} \end{document}

and the optimal control is sought to minimize the quadratic performance index as follows: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} J = & \frac { 1 } { 2 } ( x ( t_f ) - x_d ( t_f ) ) ^T S_1 ( x ( t_f ) - x_d ( t_f ) ) \\ & + \frac { 1 } { 2 } \int_0^ { t_f } \ { \mid x ( t ) - x_d ( t ) ] ^T Q_1 [ x ( t ) - x_d ( t ) ] + \dot { e } ( t ) ^T R \dot { e } ( t ) + e ( t ) ^T Q_2 e ( t ) \ } dt \tag{27} \end{align*} \end{document}

where S₁, Q₁ and Q₂ are positive semidefinite matrices and R is a positive definite matrix. The LQ problem as formulated above is concerned with tracking of the desired trajectory (x_d(t), biological data). In the performance index J, the first term penalizes the deviation of x(t_f) from the desired trajectory at the final time. Inside the integral, the first term penalizes the transient deviation of x(t), from the desired trajectory x_d(t), which represents the error dynamics. The second penalizes the change of activity of edges (dynamic graph), which attempts to minimize the variation of activity of edges over time (smoothly evolving). Also, the third term penalizes the activities of edges. Therefore, the second and third term attempt to achieve a sparse and smoothly evolving biological network. In order to use a general LQ framework, first we define \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\dot{e} ( t ) = v ( t )$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\bar{x} ( t ) = x ( t ) - x_d ( t )$$ \end{document} . Here, we assume that we know x_d(t) and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\dot{x}_d ( t )$$ \end{document} for 0 ≤ t ≤ t_f because once we have x_d(t) then we can get \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\dot{x}_d ( t )$$ \end{document} by using the derivative of a polynomial fitting. We define an (n + k) × 1-dimensional state \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$X ( t ) = [ \bar{x} ( t ) ^T , e ( t ) ^T ] ^T$$ \end{document} . Then, the state equation for the enlarged state vector can be formulated as follows: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \begin{split} \frac {d} {dt} \left[ \begin{matrix} \bar x (t) \\ e (t) \end{matrix} \right] & = \left[ \begin{matrix} 0_ { n \times n } & A (x) _ { n \times k } \\ 0_ { k \times n } & 0_ { k \times k } \end{matrix} \right] \left[\begin{matrix} \bar x ( t ) \\ e ( t ) \end{matrix} \right] + \left[ \begin{matrix} - \dot x_d ( t ) \\ v ( t ) \end{matrix} \right] = { \cal A } ( x ) X ( t ) + \left[ \begin{matrix} - \dot x_d ( t ) \\ 0_ { k \times1 } \end{matrix} \right] + \left[ \begin{matrix} 0_ { n \times1} \\ v (t) \end{matrix} \right] \end{split} \\ & = {\cal A} (x) X (t) + W (t) + V (t) \tag {28} \end{align*} \end{document}

where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \cal A} ( x )$$ \end{document} is also a linear function of x. Note that the augmented system is still a linear system because there is no multiplication between A(x) and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\bar{x} ( t )$$ \end{document} . Also, the performance index (Eq. 27) can be written as follows: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} J & = \frac {1} {2} X ( t_f ) ^T \left[\begin{matrix}S_1 & 0 \\ 0 & 0 \end{matrix} \right] X ( t_f ) + \frac {1} {2} \int_0^{t_f} \left\{ X (t) ^T \left[ \begin{matrix}Q_1 & 0 \\ 0 & Q_2\end{matrix} \right] X ( t ) + V ( t ) ^T \left[ \begin{matrix}0 & 0 \\ 0 & R\end{matrix} \right] V ( t ) \right\} dt \\ & = \frac {1}{2} X ( t_f ) ^T { \cal S} X ( t_f ) + \frac {1} {2} \int_0^{t_f} \{ X (t) ^T {\cal Q} X (t) + V (t) ^T {\cal R}V (t) \} dt & \tag{29}\end{align*} \end{document}

The problem is now reformulated as a standard LQ problem with the exception of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \cal R}$$ \end{document} , which is a singular matrix. However, we are interested in v(t), so the solution of the continuous time LQ problem is given by the state feedback control law as shown below: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} V (t) & = - {\cal R} ^ {+} P (t) X (t) = - \left[ \begin{matrix} 0 & 0 \\ 0 & R^ {- 1} \end{matrix} \right] { P } ( t ) X ( t ) = \left[ \begin{matrix} 0 \\ v_ {opt} ^* ( t ) \end{matrix} \right] \tag {30}\end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} - \frac {d {P}(t)}{dt} = {\cal A} (x) ^T {P} (t) + {P} (t) {\cal A} (x) - {P} (t) ^T {\cal R} ^ {+} {P} (t) + {\cal Q} \tag{31} \end{align*} \end{document}

where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$P ( t_f ) = { \cal S}$$ \end{document} and Equation (31) is a Riccati equation (proof in Chang et al., 2011a).

Note that the Riccati Equation (31) includes A(x) term in \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \cal A}$$ \end{document} , yet we can handle this easily by replacing x by x_d. This trick is reasonable because our optimal control input, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$v_{opt}^* ( t )$$ \end{document} , drives x(t) to x_d(t) by choosing proper Q₁, Q₂, and R. Otherwise, we would have to solve a two point boundary value problem (TPBVP) by numerical iteration.

Proposition (Chang et al., 2011a). The Riccati equation (31) can be solved by replacing x by x_d, using \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \cal Q , R}$$ \end{document} , which drive x to x_d.

Also, we can evaluate the dynamic graph e(t) by integration: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}e^* ( t ) = \int_0^{t} v_{opt}^* ( \tau ) d \tau \tag{32}\end{align*} \end{document}

Therefore, this proposed LQ optimal control framework allows us to capture pivotal development events and dynamics of the temporally evolving system.

6.3. Numerical example

Numerical examples are illustrated in Chang et al. (2011a).