New Measurement Methods of Network Robustness and Response Ability via Microarray Data

“Robustness”, the network ability to maintain systematic performance in the face of intrinsic perturbations, and “response ability”, the network ability to respond to external stimuli or transduce them to downstream regulators, are two important complementary system characteristics that must be considered when discussing biological system performance. However, at present, these features cannot be measured directly for all network components in an experimental procedure. Therefore, we present two novel systematic measurement methods – Network Robustness Measurement (NRM) and Response Ability Measurement (RAM) – to estimate the network robustness and response ability of a gene regulatory network (GRN) or protein-protein interaction network (PPIN) based on the dynamic network model constructed by the corresponding microarray data. We demonstrate the efficiency of NRM and RAM in analyzing GRNs and PPINs, respectively, by considering aging- and cancer-related datasets. When applied to an aging-related GRN, our results indicate that such a network is more robust to intrinsic perturbations in the elderly than in the young, and is therefore less responsive to external stimuli. When applied to a PPIN of fibroblast and HeLa cells, we observe that the network of cancer cells possesses better robustness than that of normal cells. Moreover, the response ability of the PPIN calculated from the cancer cells is lower than that from healthy cells. Accordingly, we propose that generalized NRM and RAM methods represent effective tools for exploring and analyzing different systems-level dynamical properties via microarray data. Making use of such properties can facilitate prediction and application, providing useful information on clinical strategy, drug target selection, and design specifications of synthetic biology from a systems biology perspective.


Introduction
Biological processes such as development, differentiation, tumorigenesis, and aging are increasingly being described in terms of temporal ordering of highly orchestrated transcriptional programs [1]. Although these processes can be analyzed using high-throughput gene expression microarray data at multiple time points, the computational methods available only identify which genes vary significantly, and how they vary across some or all of the time points measured [2,3]. Despite the technological limitations, system-level properties such as robustness increasingly attracted serious scientific interest in systems biology. Robustness -the ability to maintain continued performance against intrinsic perturbations and uncertaintyhas long been recognized as a key property of living systems and has been reviewed extensively elsewhere [4,5]. Moreover, this fundamental and ubiquitously observed system-level phenomenon cannot be understood by focusing on the individual components, even though individual components of a system may or may not be robust themselves.
The term robustness is encountered widely in very different scientific fields, from engineering and control theory to dynamic systems [5] and biology [6,7,8]. It is important to note that robustness describes a relative property, not an absolute one, because no system can maintain stability in all functions when it encounters a perturbation. In other words, robustness is not immutable. If a system cannot alter itself, it cannot adapt to intrinsic perturbations or changes to its internal environment. Most robust biological systems possess a set of certain mechanisms to achieve this: in particular, both positive and negative feedback control are ubiquitous. Negative feedback, the main mechanism used to achieve a robust response to perturbations, promotes the restoration of an initial condition of a system. Positive feedback amplifies sensitivity to changes and perturbations, a necessary feature for a cell that needs to make decisions robustly [9,10]. Therefore, feedback control is coupled with redundancy to enhance robustness. A signal transduction system is an additional network system that can be considered: it also plays an important role in responding to upstream signals or transducing environmental stimuli to downstream transcription factors in intracellular communication. From the signal energy point of view, the network response ability can be considered as system gain. Nevertheless, there is still no systematic method to measure the response ability of a network system at present.
In the last decade, with advances in experimental techniques, many researchers have utilized high-throughput data from DNA microarrays, yeast two-hybrid assays, co-immunoprecipitation, and ChIP-on-chip to study many kinds of bio-molecular networks. To this end, these kinds of data are often integrated to construct gene regulatory networks (GRNs) and proteinprotein interaction networks (PPINs). These molecular networks have been demonstrated to have great potential for revealing the basic functions and essential mechanisms of various biological phenomena; this is accomplished by understanding biological systems on a system-wide, holistic level rather than in terms of their individual components [11].
Although the concept of robustness tradeoffs has been proposed [12], few studies have quantitatively and systematically calculated the robustness and response ability of a network system. The qualitative notion that a biological system that is highly tolerant to perturbations is robust is largely uncontested; yet, there is no generally accepted approach using microarray data that can analyze robustness quantitatively and systematically from a systems biology perspective. Moreover, in a GRN or PPIN, each component can transduce the effects of external stimuli downstream, but the response ability of the network from the global point of view cannot be experimentally measured for all network components. To solve this problem, we propose estimating the response ability of a GRN or PPIN via microarray data computations. In this study, we present a new measurement method, the Network Robustness Measurement (NRM), to estimate the relative robustness of a GRN or PPIN at different time stages of biological processes and under distinct biological conditions, checking their ability to tolerate intrinsic perturbations and uncertainties based on Lyapunov stability theory. Moreover, we also propose another novel measurement method, the Response Ability Measurement (RAM), to estimate the ability of a network system to respond to and then transduce external stimuli to downstream regulators.
In addition, we illustrate our approaches by applying them to two distinct network systems: a GRN and a PPIN. First, the NRM and RAM are applied to an aging-related GRN constructed via microarray data. We show that a GRN is more robust to intrinsic perturbations in the elderly than in the young, and as such possesses poor response ability to external stimuli. We then applied these methods to PPINs of fibroblast and HeLa cells under oxidative stress conditions. We observe that the PPIN of cancer cells possesses better robustness than that of normal cells under oxidative stress. Moreover, the response ability of the PPIN calculated from the cancer cells is lower than that of the healthy cells. Since ''robustness'' and ''response ability'' are two complementary antagonistic characteristics of a network system in the context of system performance (i.e., a more robust system will be less responsive, and vice versa [13]), measurements of these systematic properties could provide verification of the system characteristics of biological networks constructed from microarray data. Measuring the properties of such trade-offs between the robustness and response ability of a network system may provide insight into systemic changes that occur under different biological conditions. These generalized NRM and RAM methods represent effective tools for exploring and analyzing different systems-level dynamical properties. Making use of such properties can facilitate not only predictions and applications in gene therapy and drug target selection from the perspective of systems biology, but also help calibrate design specifications of robust synthetic biological networks for practical applications [14,15,16]. Therefore, the proposed methods may have great potential for analyzing and designing the next generation of synthetic GRNs or PPINs.

Materials and Methods
Since we provide test cases of measuring the robustness and response ability of a GRN and a PPIN in the Results section, here we will focus on how to construct dynamic system models for a GRN and a PPIN from microarray data and other databases. Then using these dynamic models, we will propose two methods to measure their robustness and response ability.

Constructing a Dynamic Model for a Gene Regulatory Network (GRN)
It is assumed that the network of genes of interest consists of Z genes, and the time-series microarray data used for constructing the GRN contains K time points. In order to take into account the delayed effects of molecular trafficking and cellular signal transduction, the GRN can be represented by the following linear discrete-time dynamic system:  (1) can be considered as the linearized system at the operation point of interest and H will be related to the location of the operation point of linearization: e.g., for a nonlinear gene network X ½tz1~f (X ½t)zn½t for some nonlinear interaction f (X ), the linearized gene network around equilibrium point of interest X 0 can be represented as X ½tz1~AX ½tzCzDf (X )zn½t, where A~L f (X ½t) LX ½t D X ½t~X 0 and Df (X ) denotes the high order non-linear term. In this case, H~CzDf (X ) represents the basal level of the linearized system at this equilibrium point of interest. A case of a more general, nonlinear GRN will be discussed later.
Remark 1: The network parameters of system matrix A and basal level H of GRN in (1) are identified by microarray data via recursive least square estimation method [17] in Text S1.
In the following subsections, we would like to measure robustness and response ability of a PPIN for a given biological phenomenon. Therefore, we need to construct the rough PPIN by data mining and then prune it using microarray data to obtain the refined PPIN for the phenomenon.

Constructing the Rough Protein-protein Interaction Network (PPIN)
In order to obtain a PPIN for a biological phenomenon and investigate its network robustness and response ability, we first need to reconstruct the rough PPIN from microarray data and protein-protein interaction (PPI) information: in this case, our data was taken from STRING (Search Tool for the Retrieval of Interacting Genes/Proteins). STRING, a database of known and predicted protein interactions, comprises direct (physical) and indirect (functional) associations mined from four sources: genomic context, high-throughput experiments, conserved co-expression and previous knowledge. The database currently covers 5,214,234 proteins from 1,133 organisms. As a preliminary measure, we first perform a Bonferroni correction to select those genes with significantly differential expression (the corrected P-value = 0.01/ n, where n = 25,802 denotes the number of genes) from microarray data for inclusion in the target pool. Bonferroni correction, one of the strictest and most conservative multiple testing corrections, adjusts the individual P-value for each gene to keep the overall error rate less than or equal to the user-specified P-cutoff value when performing a statistical test on a group of genes [18,19]. Subsequently, this target pool and the PPI database (STRING) are integrated to indicate the possible PPIs between two selected candidates, and the rough PPIN can then be constructed based on the selected candidates in the protein pool and the potential PPIs among them.

Pruning the Rough PPIN via a Dynamic Model
Modern large-scale studies of PPIs permit the development of protein interaction networks, but all large-scale experiments and databases to date have returned high rates of false-positives [20]. Previous evidence has clearly demonstrated that utilizing multiple functional databases not only allows for better identification of PPIs, but also leads to better prediction of the function of unknown proteins [21]. The rough PPIN constructed here comprehensively details the protein interactions possible under all kinds of biological situations and experimental conditions. As such, these interactions must be further narrowed down by microarray data so that they represent only the protein interactions appropriate for the conditions under consideration. Therefore, we use a dynamic PPI model and model order selection method, Akaike Information Criterion (AIC) [22], together to prune the rough PPIN using time series microarray data, deleting false positive PPIs in our PPIN. Here, a dynamic PPI model for a rough PPIN can be represented as follows: where v i ½t denotes the expression profile of the ith target protein v i ½t in the rough PPIN at time point t; v i [t] denotes the model residue and measurement noise; a ii denotes the influence of the ith target protein on itself at the next time point; b ij denotes the individual interactive ability of protein j in the rough PPIN with target protein v i ½t when i=j, Vi,j~1 Á Á Á Q and the basal interaction k i represents unknown PPIs or other influences such as mRNA-protein interactions or protein synthesis. Therefore, the PPIN with dynamic protein interactions in (2)   where R denotes the interactive matrix of the PPIN, K denotes the basal levels of the PPIs, and V[t] denotes the modeling and measurement errors.

Remark 2:
The protein interaction model in (2) can be used to describe any target protein in the rough PPIN. Therefore, the whole PPIN can be described by the protein interaction model described in (3).
In equation (2), v i ½t denotes the expression profile of the ith target protein at time point t, which is estimated from mRNA expression profiles through the following translational sigmoid function [23]: where r denotes the transition rate of the sigmoid function, M i denotes the mean mRNA expression level of the corresponding protein i, and x i ½t represents the corresponding mRNA expression profiles of the i-th target gene in the rough PPIN.

Remark 3:
The network parameters R and K of PPIN in (3) are also identified by the recursive least square estimation method [17]. Further a system order detection method, i.e. Akaike Information Criterion (AIC), is employed to detect the number Q of protein interactions in interactive matrix R in (3), i.e. use AIC to prune the false positive protein interactions in R (see Text S2).

Network Robustness Measurement (NRM) of a Dynamic Network System
After constructing the dynamic models of the GRN in (1) and the PPIN in (3) from time series microarray data by the recursive least square estimation method as detailed above, we will measure their network robustness based on these dynamic models. The steady state X s of the dynamic system in (1) is obtained as: To simplify analysis of ''robustness'' of the steady state (phenotype) of the GRN, the origin of the dynamic system is shifted to the steady state X s (i.e., X ½t~X X ½tzX s ), achieving the following shifted dynamic system by subtracting equation (1) from equation (5) [24].X X ½tz1~AX X ½t ð 6Þ Therefore, the robustness of the network's steady state becomes the robustness of the shifted network in (6) at the originX X ½t:0: In the shifted dynamic system in (6), if the eigenvalues of A are all inside the unit circle DZDƒ1 in the complex Z domain, theñ X X ½t?0, and X ½t?X s , or approaches the steady state asymptotically: that is, the phenotype will still be maintained after a small perturbationX X ½t. Since robustness of the dynamic system in (1) at the steady state X s is equivalent to robustness of the dynamic system in (6) at the originX X ½t:0,the shifted dynamic system in (6) will simplify the robustness analysis procedure of a network system [24,25].

Remark 4:
The steady state Ws of the PPIN in (3) is obtained as: The origin of the dynamic PPIN is shifted to the steady state W s (i.e., W ½t~W W ½tzWs), achieving the following shifted dynamic system [24].W W ½tz1~RW W ½t ð 7Þ Therefore, the robustness of the PPIN phenotype becomes the robustness of the above shifted PPIN at the originW W ½t:0: We then use this shifted PPIN in (7) to estimate its robustness and response ability in the same way as the GRN in (6) by following the same network robustness measurement (NRM) and response ability measurement (RAM) methods. Bearing this in mind, we only discuss the NRM and RAM of the GRN in (6). The measurements of NRM and RAM of a PPIN follow the same procedure.
The robustness of the shifted network in (6) is a measure of its ability to tolerate intrinsic molecular-level perturbation [13]. If the linear network system suffers from these kinds of perturbations, whether due to process noise, thermal fluctuations or genetic mutations, the interactive matrix A is perturbed as A(1zg), where g denotes the ratio of intrinsic perturbation, and the corresponding additional system perturbation DA~gA [13]. Then the network system with intrinsic perturbation considerations factored in can be represented by:X A higher g value means more intrinsic perturbations, and the larger the g value that a network can tolerate, the greater robustness the network possesses. Since state-dependent intrinsic perturbations can influence the stability of a network system, we must further discuss the robustness of the perturbative network system in (8). According to quadratic stability theory [25,26], if the Lyapunov equation of energy function for the perturbative network system in (8) is chosen as VX X t ½ À Á~X X T t ½ PX X t ½ for a positive symmetric matrix P~P T w0, then the system in (8) is : that is, if the energy of the GRN is not increased by intrinsic perturbations. Based on this idea, we are able to obtain the following robust stability principle for a GRN of interest that is subject to intrinsic perturbations.
Proposition 1: Suppose a network system suffers from intrinsic perturbations as described in (8). The perturbative network is robustly stable if the following inequality has a positive definite solution P~P T w0: Proof: See Text S3. Remark 5: (i) If a network system is free of intrinsic perturbations as in (6) (i.e., g~0), the stable condition is reduced to the following matrix inequality, A T PAƒP , which has a positive definite solution Pw0 [13]. In order to guarantee the stability of a perturbation-free network system, eigenvalues of the system interaction matrix A should be inside the unit circle DZDƒ1 in the Z complex domain: i.e., A T Av1. From (9), it can be seen that when g is increased gradually, eigenvalues of system interactive matrix A of the perturbation system in (8) should be nearer the origin of the complex domain Z than those from the perturbation-free system in (6) in order to tolerate more intrinsic perturbation g. Some eigenvalues of system interactive matrix A near the unit circle are more easily perturbed outside the unit circle by intrinsic fluctuations, making the perturbative network system unstable. Therefore, the distance between eigenvalue locations of A to the unit circle represents a measurement of robustness for a linear network system. (ii) Based on the above analysis, the robustness g o of our network system is the maximum perturbation tolerance allowed by the network, and it can be measured by solving the following constrained optimization: i.e., the maximum perturbation ratio g o tolerable by the network system while maintaining stability. The constrained optimization problem in (10) can be solved by increasing g until no positive solution Pw0 exists in (9): i.e., to the highest possible g without violating the robust stability in (9). A positive definite solution Pw0 in (9) can be easily obtained by using the Linear Matrix Inequality (LMI) Toolbox of Matlab.
We can solve the above constrained optimization problem to measure the robustness of the network system simply by increasing g until no positive definite solution Pw0 can be found in (9), and then the maximal perturbation g o tolerated by the network system represents the robustness under intrinsic perturbations. Hence, we are able to obtain the relative robustnesses of network systems we would like to investigate under many different biological circumstances. When comparing g o for each network, the higher the value of g o , the greater the robustness estimate for its respective network.

Response Ability Measurement (RAM) of a Dynamic Network System
In this subsection, we turn to how to estimate the response ability of a network system to external stimuli after its robustness has been calculated. Assuming a network system responds to external stimuli U t ½ , which might include upstream regulatory signals and external signals from factors such as carcinogens, oxidative stress and ambient pro-inflammatory molecules, the dynamic network (6) is then modified as follows [27].
where U t ½ ~½u 1 ½tu 2 ½t . . . u Z ½t T represents external stimuli and Y t ½ denotes the output signal responses of the network system of interest. For example, if the output signal response of all genes within the network system is analyzed, then C~I, an identity matrix. If only the response of the last gene to U t ½ is analyzed, the : i.e., all elements of C are zero except for the last element. In a network system, each component may amplify or attenuate the effect of external stimuli when communicating downstream. The effect of input signals U t ½ on output signals Y ½t is less than or equal to a positive value d, if the following inequality [25,26]: holds for all possible bounded input signals U½t=0 andX X ½0~0.
The physical meaning attributable to (12) is that the effect of all possible external stimuli U½t on Y ½t is less than or equal to d 2 from the energy point of view: i.e., d denotes the upper bound of the effect of U½t on Y ½t. IfX X ½0=0, then inequality (12) should be modified as the following [25,26]: to account for the effect ofX X ½0 on Y ½t. d o , the smallest upper bound of d in (12), is called the ''network response ability'' of the system to all bounded stimulus signals U½t. From d o , we obtain a more systematic insight into the ability of a network to respond to external stimuli, which is useful when considering signal transduction or intracellular communication. Via this method, we can measure and compare network response abilities d o of various network systems under different biological conditions and states. Proposition 2: The response of the network system in (11) has an upper bound d in (12), if there exists a positive definite P~P T w0 solution to the following linear matrix inequality (LMI) [28]: i.e., if the above LMI holds for some Pw0, then the effect of U t ½ on Y ½t must be less than or equal to d: i.e., (12) or (13) hold for the network system in (11).
Proof: See Text S4. Remark 6: (i) The response ability d o of a network system to external stimuli can be obtained by minimizing its upper bound d, achieved by solving the following constrained optimization problem: subject to dw0,Pw0, and (14) ð15Þ This constrained optimization can be solved by decreasing the upper bound d in (12) until no Pw0 exists in (14) using Matlab LMI toolbox [27].
(ii) d o , the smallest upper bound of d in (12), is called the network response ability of a network system to external stimuli. If d o v1, it means external stimuli are attenuated (buffered) by the network system; if d o w1, external stimuli are amplified. (iii) Since the output/input signal energy ratio (i.e., the response) in (12) or (13) is considered for all possible external stimuli, the measure of network response ability in (14) and (15) is dependent more on system matrix A than on the external stimuli. This is analogous to the way a lowpass filter is dependent more upon the characteristics of the filter than on noise. Therefore, it is possible to measure the network response ability of a network system from its system characteristics without knowledge of what the external stimuli are, with the caveat that these signals must be bounded. The system matrix A can be estimated from timeseries microarray data as in the previous subsections. Therefore, employing the Matlab LMI toolbox, we can use (15) to measure network response abilities using microarray data under different specific experimental conditions.

Extension of Network Robustness and Response ability Measurements to Nonlinear Network Systems
Suppose the GRN or PPIN could be constructed by microarray data for the following nonlinear network system: Suppose we are interested in the phenomenon around the equilibrium point Xei.e., Xe~f (Xe) -and letX X ½t~X ½t{Xe. Then, we get the shifted nonlinear network system as: In this system,X X ½t:0 is the equilibrium point of interest. Suppose the nonlinear network in (17) could be globally linearized as described in [27] for allX X ½t~X ½t{Xe as: where V[R nxn denotes the system parameter set of these globally linearized systems. Suppose V could be described by a list of its vertices via the following convex hull: By the global linearization technique described in [27], all the globally linearized systems in V of the nonlinear system in (16) or (17) could be interpolated by the m linear system X X ½tz1~AiX X ½t,i~1,2, Á Á Á m at the vertices of convex hull of the polytope: where the interpolation functions h i (X X ), i~1, Á Á Á m have the following properties: i.e., the trajectory of the nonlinear system in (16) or (17) could be interpolated by the trajectories of m linear systems in (20) via the nonlinear interpolation functions h i (X X ), i~1, Á Á Á m Suppose the nonlinear network system suffers from the following intrinsic perturbations: Based on the robust stability principle, we deduce the following. Proposition 3: Suppose the nonlinear network system suffers from intrinsic perturbations as in (22). The perturbative nonlinear network system is robustly stable if the following LMIs have a positive definition P~P T w0 : Proof: See Text S5. The physical meaning of the LMIs in (23) is that the perturbation g is tolerable within all linearized network systems. Therefore, based on the global linearization method, the robustness g o of the nonlinear network system in (16) can be measured by solving the following constrained optimization: g o~m ax g subject to Pw0 and LMIs in (23) i.e., the maximum perturbation g o , which is tolerable by all linearized network systems, is a global extension of (15) from a linear network system to a nonlinear network system. Similarly, if the shifted nonlinear network system with output Y ½t suffers from external stimuli as in: then, based on the global linearization technique detailed in [27], the shifted nonlinear network system of (25) with output Y ½t can be represented by the following interpolated system: Then we get the following result: Proposition 4: The network response of the nonlinear network system in (26) has an upper bound d in (12), if there exists a positive definite solution P~P T w0 to the following LMIs: i.e., if the above LMIs hold for some Pw0, then the effect of U½t on Y ½t must be less than or equal to d: That is, (12) or (13) holds for the nonlinear network system (25) or (26).
Proof: See Text S6. Therefore, based on the global linearization method, the network response ability d o of the nonlinear network system in (25) or (26) can be measured by solving the following constrained optimization:

Remark 7:
If the nonlinear network system is constructed by the system identification method from microarray or other experimental data, then the network robustness g o and response ability d o can be measured by (24) and (28), respectively. However, to approximate the nonlinear network system via interpolation, these systems must undergo the global linearization in (18) and (19) to find the vertices of the convex hull of the polytope in order to interpolate this linearized system at the vertices (as in the system in (26)). This is generally a very complex task. If we only want to measure the network robustness and response ability of a network system near an equilibrium point (e.g. phenotype), then the measurement methods for network robustness and response ability by (10) and (15) based on a linearized system at its equilibrium point are sufficient. In the following section, we give several practical examples based on linear network systems for measuring the network robustness and response ability of a GRN and PPIN constructed by time-series microarray data.

Results/Discussion
In this paper, we have presented two novel computational methods, Network Robustness Measurement (NRM) and Response Ability Measurement (RAM), which estimate the robust-ness and response ability of network systems using time-series microarray data. Briefly, the basic methodologies consist of integrating time-series microarray data with either experimental literature or protein-protein interaction information to construct the dynamic system of a gene regulatory network (GRN) or protein-protein interaction network (PPIN). Both of these systems correspond to real phenomena in living organisms and can be used to further evaluate network robustness and response ability (see Figure 1). Details have been described in the preceding Materials & Methods section. These methods are equally applicable to and useful for GRNs and PPINs. In order to demonstrate the potential applications of NRM and RAM for estimating system characteristics exhibited under different biological conditions, we elaborate on the usage of the above methods in two separate cases and explain the analytic results within the following subsection.

Case 1: Estimating the Robustness and Response Ability of a Gene Regulatory Network (GRN) at Different Stages of Life
Aging is an extremely complex and system-level process, and it has attracted much attention in medical research, especially since chronic diseases are quite prevalent in elderly populations. These illnesses may be the result of both gene mutations that lead to intrinsic perturbations and environmental changes that stimulate signaling in the body. In addition, aging may result in reduced responses to environmental stimuli such as carcinogens, oxidative stress, and pro-inflammatory molecules. Therefore, measuring robustness to intrinsic perturbations and response ability to external stimuli of an aging-related GRN could provide key insights into the system-level mechanisms involved in aging.
given GRN, the more robust the GRN is.
Construction of a gene regulatory network associated with aging-related pathophysiological phenotypes. Combing data from several references in the literature, we construct a GRN with multiple regulatory loops that are highly associated with agingrelated pathophysiological phenotypes, consisting of the following sixteen genes: FOXOs, NF-kB, p53, SIRT1, HIC1, Mdm2, Arf1, PTEN, PI3K, Akt, JNK, IKKs, IkB, BTG3, E2F, and ATM ( Figure 2). All genes of interest are critical for regulating cellular functions related to longevity, including detoxification of reactive oxygen species (ROS), cell cycle arrest, repair of damaged DNA, apoptosis, and pro-inflammation, as age increases [29,30]. The GRN can be represented by the linear discrete-time dynamic system shown in Table 1.
In this case, we want to estimate the network robustness and response ability of this GRN, and also to investigate whether they might be altered with increasing age. Furthermore, those altered system characteristics associated with a high incidence of agingrelated disease, such as cancer, are prime targets for discussion. The high-throughput microarray data assembled in [31], which details the age-related effects on gene expression in the thymus and spinal cord dissected from male C57BL/6 mice of ages 1, 6, 16, and 24 months, are used as our gene expression profiles. In order to make this data more useful for identifying the system parameters of our dynamic GRN, we first rescale the original gene expression time profiles. For young mice, those aged 17d and 19d are categorized as 3-week-old animals; 40d, 43d, and 44d as 6-week-old animals; 174d, 180d, and 186d as 25-week-old animals; and, 193d and 194d as 28-week-old animals. For aged mice, those aged 476d, 481d, and 485d are categorized as 68week-old animals; 495d and 498d as 71-week-old animals; 714d and 719d as 102-week-old animals; and, 730d, 733d, and 743d as 105-week-old animals. (All animals had been placed on ad libitum diets.) All of the above expression level data are calculated in log2 scale, and then the Z-score method is used to normalize the gene array expression data as has been previously described [32].
We construct multiple regulatory loops of a GRN containing sixteen genes ( Figure 2) associated with aging-related pathophysiological phenotypes. Using the recursive least square parameter estimation method [17] and the microarray data from the thymus and spinal cord [31], a dynamic model of the sixteen genes within the GRN in equation (1) is constructed (shown in Table 1) with estimated parameters at different stages of life (summarized in Table 2). With equation (1) and Tables 1-2, we calculate sixteen eigenvalues for interactive matrix A at the young and aged stages. The result (Figure 3) shows some eigenvalues of the interactive matrix A clustered together in similar regions near the point Z = 1, demonstrating that some modes of the same network will more readily approach some constant steady state than others [12,33].
Estimating network robustness and response ability of GRNs at different stages of life by NRM and RAM methods. It has previously been shown that gene networks associated with aging-related pathophysiological phenotypes suffers from intrinsic perturbations mainly due to noise processes like molecular fluctuation or genetic mutation. The robustness of the GRN is defined as its ability to tolerate intrinsic perturbations and plays a principal role as a fail-safe mechanism during evolutionary processes [12]. Therefore, comparisons of robustness g o for different GRNs pertaining to young and aged stages may quantitatively describe their relative ability to tolerate these intrinsic perturbations. By using our NRM method, we evaluate the network robustness as shown in Table 3 to represent the maximum intensity of intrinsic perturbations that the network system can tolerate. In other words, the higher the g o calculated from a Concerning the RAM, if a GRN also responds to external stimuli U½t, such as upstream signals or biological molecules outside the network, as in equation (11), the network response ability measured should be based on the output/input energy ratio within a time interval as detailed in equation (12). Supposing the ratio of the effect of input signals U t ½ on Y ½t is upper-bounded by a positive value d , then we can calculate the smallest value d o in equation (12) (i.e., the minimum upper bound), and call it the ''network response ability'' to external stimuli. This method, RAM, is based on system gain theory, and incorporates principles behind the external input signal in system theory perspectives [27]. In general, a GRN with greater robustness always has lower response ability, and vice versa. This is an inherent characteristic of physical systems [27]. The trade-off between robustness and response ability is demonstrably visible in Table 3 and constitutes secondary proof of our proposed method.
Comparing the GRNs of the thymus and spinal cord data between the ''young'' group and the ''aged'' group, we show that both tissues evaluated have a higher robustness g o value in the aged group than the young group. These results are shown in Table 3. The aging-related GRN examined in the first case study is less robust to intrinsic perturbations in the young group than in the aged group. In other words, the GRN at the aged stage can function in spite of the accumulation of DNA damage or genetic mutations accompanying the aging process. In addition, response ability (d o ) of the GRN calculated from the aged group is lower than that from the young group (Table 3A). In order to tolerate and survive the accumulated intrinsic perturbations (e.g., genetic mutations) associated with the increase in age, the elderly GRN becomes more robust than the young GRN. Therefore, elderly GRNs are less able to respond to external stimuli or transduce them downstream. Hence, the response ability protecting organismal function against external stimuli (e.g., pro-inflammatory molecules) may become worse in the elderly GRN. To confirm the power of the proposed methods, the proposed network robustness and response ability is validated by computer simulation in Figure  S1-S3 in Text S7.
In addition, a more robust GRN could harbor more accumulated genetic mutations, which might provide more evolutionary paths to other GRN phenotypes via random drift, and even lead to aging-related diseases like cancers, metabolic disorders, and arthritis [34,35]. On the other hand, the young aging-related Table 1. Ordinary differential equations of gene regulatory networks in Figure 2 for sixteen genes associated with aging-related pathophysiological phenotypes.

Gene
Dynamic Equations of Sixteen Genes of Interest GRN is less robust, with a greater response to external stimuli. From the above data, we can conclude that the GRNs of the young group more efficiently respond to external stimuli. These values show that while the GRN of the young group is less robust to intrinsic perturbations, external stimuli elicit a strong response from it, suggesting that gene expression may be readily reprogrammed to mediate downstream genes or regulators. This also allows for modulation of gene expression in response to external stimuli such as exposure to oxidative stress, carcinogens, and pro-inflammatory molecules.

Case 2: Estimating the Robustness and Response Ability of a Protein-protein Interaction Network (PPIN) between Normal and Cancer Cells
Robustness, an intrinsic and systematic property of complex biological systems, enables living organisms to maintain their functions in the face of various perturbations to assure their survival [4,12,36]. Moreover, most genes code for sensors, actuators, and the complex regulatory networks that control them, and thus confer robustness to the cell under various challenging circumstances rather than providing merely the basic functionality required for survival under a single, ideal environment. Normal as well as cancer cells utilize the same processes. At present, although researchers have contended that cancer is a highly robust disease, and clinical observations suggest that inefficiencies in cancer therapy may be due to the according drug resistance, a useful computational method to evaluate the network robustness from microarray data and compare it between cancer and normal cells remains lacking. Therefore, in the second case, the NRM method we have proposed is applied to calculate the PPIN robustness of HeLa cells versus normal fibroblasts under menadione treatment (i.e., exposed to oxidative stress) from corresponding microarray data and protein interaction information. The network response ability is also evaluated and discussed.

Constructing refined PPINs of normal fibroblast cells and
HeLa cancer cells under menadione treatment. We use microarray data downloaded from the Gene Expression Omnibus (GEO) database on the NCBI website in this study to investigate the network robustness and response ability of a PPIN in normal cells versus cancer cells under oxidative stress [37]. Normal human diploid lung fibroblasts and HeLa S3 cervical carcinoma cells are treated with 10 mM menadione bisulfate for 0.5, 1, 3, 4, 6, 8, 12, 24, and 36 hours and 0. 5,1,2,4,8,12,24, and 32 hours, respectively. For the oxidative stress condition, the temporal responses at several time points over an interval of 24-36 hours are measured, allowing stress effects to be observed over an interval longer than a normal cell cycle duration.
As shown in the flowchart in Figure 1, the microarray data are integrated with database information concerning protein interactions as the input for our proposed method (see Materials & Methods). The candidate protein interactions are taken from all experimental conditions in these data to construct the rough PPIN. To proceed further, the corresponding microarray data are used to identify the parameters in equation (8) via the recursive least square estimation method [17], and then to prune down the rough PPINs based on these identified parameters, thus obtaining the refined PPINs of both fibroblast and HeLa cells under menadione treatment via the AIC (shown in Figure 4) [22]. After obtaining the refined PPINs, we can see there are still 161 candidate proteins with 371 interactions in the fibroblast PPIN and 222 candidate proteins with 607 interactions in the HeLa PPIN under oxidative stress conditions. Furthermore, the respective system interactive matrices R of the fibroblast and HeLa cell PPINs can be acquired and then applied to the sequential evaluation of network robustness by equation (10) and response ability by equation (15).  Table 2. Estimated parameters of gene regulatory networks in Table 1 with sixteen genes in the thymus (A) and spinal cord (B) at the young and aged stages.  robustness g o value than the PPIN of fibroblasts using our NRM method. These results are shown in Table 3B. The PPIN is more robust to such perturbation in the cancer cells than in the normal cells. Moreover, the response ability (d o ) of the PPIN calculated from the HeLa cells is lower than that from the fibroblasts (Table 3B). Therefore, the PPIN of cancer cells is less efficient in transducing or responding to external stimuli when compared with that of normal cells. The analytic results support previous contentions that cancer cells are more robust than normal cells [38]. Accordingly, compelling evidence has previously suggested that the level of reactive oxygen species (ROS) is higher in many types of cancer cells than in their normal counterparts [39,40,41]. Although a moderate increase in ROS can promote cell proliferation and differentiation, excessive amounts can cause oxidative damage to DNA, protein, and lipids. However, during tumorigenesis or advanced malignant transformation, oncogenic signals can elevate cellular levels of ROS to the extent that numerous genes are altered, accumulating changes via both genetic and epigenetic mechanisms [42,43]. Therefore, cancer cells themselves actively generate high levels of endogenous ROS and are persistently exposed to them, stimulating cell proliferation through redox-sensitive transcription factors and promoting antioxidant adaptive mechanisms to minimize oxidative damage [44,45]. Therefore, the intrinsic oxidative stress generated during tumor initiation or malignant transformation may exert selective pressure on the cell population and enrich those capable of stress adaptation: i.e., cancer cells. These cells that survive oxidative stress may have adapted mechanisms to counteract the cytotoxic effects of increased ROS and to promote pathways essential to cellular survival [45]. All of these alterations confer cellular resistance and system robustness to subsequent intrinsic or extrinsic oxidative stress, even at high levels. In clinical cancer therapy, one chemotherapeutic strategy is to elevate excess cellular oxidative stress in order to deplete the antioxidant system and overwhelm the cancer system's ability to neutralize and eliminate such insults, inducing apoptosis and destroying the cancer network system. However, the efficiency of such chemo-or radiotherapy is still unsatisfactory. Moreover, previous evidence has suggested that tumors are highly robust and maintain their proliferative potential against a wide range of anti-cancer therapies. Heterogeneities have been regarded as a vital property of cancer cells because they facilitate robustness through redundancy, and so subsystems killed by chemotherapy can be functionally replaced to ensure tumor proliferation and survival [38]. Hence, the notion of the PPIN of HeLa cells having higher robustness, as evaluated by our proposed NRM method, is consistent with previous contentions and experimental observations. Utilizing microarray data in systemlevel computations, we have sufficiently demonstrated that cancer cells possess better robustness and less response ability under oxidative stress when compared with normal cells.

Conclusions
Our proposed methods, NRM and RAM, consist of multiple analytical steps, from constructing dynamic GRNs and PPINs via microarray data to estimating network robustness and response ability based on their dynamic network systems. However, each individual step can be employed separately, depending on the specific application in question. For example, in the first step we select a set of differentially expressed genes in aging (Case 1) or cancer (Case 2) as our target genes. However, another user could replace this set of genes with other sets of genes of interest to facilitate targeted studies over a wide range of issues. In addition, our proposed methods can evaluate the network robustness and response ability not only for GRNs (Case 1) but also for PPINs (Case 2) depending on the needs of the user. The generalized NRM and RAM methods are effective tools for exploring and analyzing different systems-level dynamical properties. Making use of such network properties can facilitate biomedical predictions and applications, providing useful information for clinical strategy, drug target selection, and design specifications of synthetic biology from the perspective of systems biology [10,14,15,16].

Supporting Information
Text S1 Parameter Identification of a Gene Regulatory Network by Microarray Data.