BioMed Research International

Volume 2013 (2013), Article ID 703849, 15 pages

http://dx.doi.org/10.1155/2013/703849

Research Article

A Systems’ Biology Approach to Study MicroRNA-Mediated Gene Regulatory Networks

Xin Lai,1,2 Animesh Bhattacharya,³ Ulf Schmitz,¹ Manfred Kunz,³ Julio Vera,² and Olaf Wolkenhauer1,4

¹Department of Systems Biology and Bioinformatics, University of Rostock, 18051 Rostock, Germany

²Laboratory of Systems Tumor Immunology, Department of Dermatology, Faculty of Medicine, University of Erlangen-Nuremberg, Ulmenweg 18, 91054 Erlangen, Germany

³Department of Dermatology, Venereology and Allergology, University of Leipzig, 04155 Leipzig, Germany

⁴Institute for Advanced Study (STIAS), Wallenberg Research Centre at Stellenbosch University, Stellenbosch 7600, South Africa

Received 31 July 2013; Revised 12 September 2013; Accepted 17 September 2013

Academic Editor: Tao Huang

Copyright © 2013 Xin Lai et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

MicroRNAs (miRNAs) are potent effectors in gene regulatory networks where aberrant miRNA expression can contribute to human diseases such as cancer. For a better understanding of the regulatory role of miRNAs in coordinating gene expression, we here present a systems biology approach combining data-driven modeling and model-driven experiments. Such an approach is characterized by an iterative process, including biological data acquisition and integration, network construction, mathematical modeling and experimental validation. To demonstrate the application of this approach, we adopt it to investigate mechanisms of collective repression on p21 by multiple miRNAs. We first construct a p21 regulatory network based on data from the literature and further expand it using algorithms that predict molecular interactions. Based on the network structure, a detailed mechanistic model is established and its parameter values are determined using data. Finally, the calibrated model is used to study the effect of different miRNA expression profiles and cooperative target regulation on p21 expression levels in different biological contexts.

1. Introduction

Although microRNAs (miRNAs) are physically small, they have been shown to play an important role in gene regulation [1]. Currently, an increasing number of studies are being carried out to deepen our understanding of miRNA regulatory mechanisms and functions. However, experimental approaches have limitations when dealing with complex biological systems composed of multiple layers of regulation such as the transcriptional and post-transcriptional regulation by transcription factors (TFs) and miRNAs [2]. Most experimental approaches focus on the identification of miRNA targets and the investigation of physiological consequences when perturbing miRNA expressions but are unsuited to provide a system-level interpretation for observed phenomena. Therefore, the introduction of a systematic approach, which can unravel the underlying mechanisms by which miRNAs exert their functions, becomes increasingly appealing.

The systems biology approach, combining data-driven modeling and model-driven experiments, provides a systematic and comprehensive perspective on the regulatory roles of miRNAs in gene regulatory networks [3–5]. To investigate a gene regulatory network, an iterative process of four steps is needed. (I) Biological network construction: a map is constructed that shows interactions among molecular entities (such as genes, proteins and RNAs), using information from literature and databases. (II) Model construction: depending on the biological problem investigated and experimental data available, the interaction map can be translated into a detailed mechanistic model that can simulate the temporal evolution of molecular entities. The values of parameters in this model can be determind from literature, databases or they are directly estimated from quantitative experimental data using optimization methods. (III) Computational experiments: once a model is established, it can be simulated and/or analyzed for its general behavior. (IV) Experimental validation: model predictions together with biological explanations are integrated to guide the design of new experiments, which in turn validate or falsify the model. If model predictions are in agreement with the experiments, the model justifies the biological hypotheses behind it. In turn, these hypotheses, which provide reasonable explanations for the biological phenomenon, lead to an enhanced understanding of the gene regulatory network. Otherwise, the structure of the mathematical model is modified to generate new hypotheses and suggest new experiments.

The application of the systems biology approach to the analysis of a gene regulatory network is demonstrated with a case study of the regulation of p21 by multiple miRNAs [4]. The network combining putative targets of TF and miRNA regulation with experimentally proven molecular interactions was constructed and visualized. Next, the network was translated into a detailed mechanistic model, which was characterized and validated with experimental data. Finally, the integration of quantitative data and modeling helped us to generate and validate hypotheses about mechanisms of collective miRNA repression on p21.

2. Results

2.1. The Systems Biology Approach to Study miRNA-Mediated Gene Regulatory Networks

2.1.1. Mathematical Modeling

The aim of our analysis is to unravel the complex mechanisms by which gene regulatory networks involving miRNAs are regulated. We iteratively integrate data from literature, experiments and biological databases into a detailed mechanistic model of a gene regulatory network. The model is then used to formulate and test hypotheses about mechanisms of miRNA target regulation and cellular process-related variability. The methodology includes four steps which are briefly summarized in Table 1. In the coming sections we present these steps in detail.

Table 1: Overview of the methodology. Key points in each step of the methodology and the main resources for constructing miRNA-mediated gene regulatory networks are given.

Step 1: data retrieval

Regulation types
Resources

Transcriptional gene regulation
TRED (http://rulai.cshl.edu/cgi-bin/TRED/tred.cgi?process=home): a database that provides an integrated repository for both cis- and transregulatory elements in mammals

TRANSFAC (http://www.gene-regulation.com/pub/databases.html): a database that collects eukaryotic transcriptional regulation, comprising data on TFs, their target genes, and binding sites

The UCSC table browser (http://genome.ucsc.edu/): a popular web-based tool for querying the UCSC Genome Browser annotation tables

HTRIdb (http://www.lbbc.ibb.unesp.br/htri/): an open-access database for experimentally verified human transcriptional regulation interactions

MIR@NT@N (http://maia.uni.lu/mironton.php/): an integrative resource based on a metaregulation network model including TFs, miRNAs, and genes

PuTmiR (http://www.isical.ac.in/~bioinfo_miu/TF-miRNA.php): a database of predicted TFs for human miRNAs

TransmiR (http://202.38.126.151/hmdd/mirna/tf/): a database of validated TF-miRNA interactions

miRGen 2.0 (http://diana.cslab.ece.ntua.gr/mirgen/): a database of miRNA genomic information and regulation

Posttranscriptional gene regulation
miRecords (http://mirecords.biolead.org/): a resource for animal miRNA-target interactions

Tarbase (http://www.microrna.gr/tarbase/): a database that stores detailed information for each miRNA-gene interaction, the experimental validation methodologies, and their outcomes

miRTarBase (http://mirtarbase.mbc.nctu.edu.tw/): a database that collects validated miRNA-target interactions by manually surveying the pertinent literature

miRWalk (http://www.umm.uni-heidelberg.de/apps/zmf/mirwalk/): a comprehensive database that provides information on miRNAs from human, mouse, and rat, on their predicted as well as validated binding sites in target genes

Protein-protein interaction
HPRD (http://www.hprd.org/): a centralized platform to visually depict and integrate information pertaining to domain architecture, posttranslational modifications, interaction networks, and disease association for each protein in the human proteome

STRING (http://string-db.org/): a database of known and predicted protein interactions. The interactions include direct (physical) and indirect (functional) associations

MPPI (http://mips.helmholtz-muenchen.de/proj/ppi/): a collection of manually curated high-quality PPI data collected from the scientific literature by expert curators

DIP (http://dip.doe-mbi.ucla.edu/dip/Main.cgi): a catalog of experimentally determined interactions between proteins

IntAct (http://www.ebi.ac.uk/intact/main.xhtml): a platform that provides a database system and analysis tools for molecular interaction data

Reactome (http://www.reactome.org/): an open-source, open access, manually curated, and peer-reviewed pathway database

GO annotation
Amigo GO (http://amigo.geneontology.org/cgi-bin/amigo/go.cgi): the official GO browser and search engine

miR2Disease (http://www.mir2disease.org/): a manually curated database that aims at providing a comprehensive resource of miRNA deregulation in various human diseases

miRCancer (http://mircancer.ecu.edu): a miRNA-cancer association database constructed by text mining on the literature

PhenomiR (http://mips.helmholtz-muenchen.de/phenomir/): a database that provides information about differentially expressed miRNAs in diseases and other biological processes

miRGator (http://mirgator.kobic.re.kr/): a novel database and navigator tool for functional interpretation of miRNAs

miRó (http://ferrolab.dmi.unict.it/miro): a web-based knowledge base that provides users with miRNA-phenotype associations in humans

Step 2: network construction and visualization

(i) Visualize regulatory interactions in platforms such as CellDesigner and Cytoscape that support standardized data formats

(ii) Calculate confidence scores for assessing reliability of interactions in gene regulatory networks

Step 3: model construction and calibration

(i) Formulate equations using rate equations

(ii) Fix parameter values using available biological information

(iii) Estimate the other unknown and immeasurable parameter values using optimization methods which can minimize the distance

between model simulations and experimental data such as time course qRT-PCR and western blot data

Step 4: model validation and analysis

(i) Design new experiments and generate new data to verify the calibrated model

(ii) Study complex properties and behavior of the system

( 1)  Data Retrieval. To construct a gene regulatory network composed by different levels of regulation, we collect information from different resources which are briefly described below, and more resources for data retrieval are introduced in Table 1.

(a)Transcriptional level regulators. Experimentally verified TFs for a gene can be extracted from literature or databases such as TRED, TRANSFAC, or HTRIdb [6–8]; the putative TFs, which are associated with conserved TF-binding sites residing in the promoter region of a gene, can for example, be extracted from the table of TFs with conserved binding sites in the UCSC genome browser or the TRANSFAC database [9]. miRGen 2.0 is a database that provides both predicted and experimentally verified information about miRNA regulation by TFs [10].(b)Post-transcriptional regulators. Databases such as miRecords, Tarbase and miRTarBase are repositories of experimentally validated miRNA:gene interactions [11–13]. Predictions of miRNA:gene interactions are accumulated in databases like miRWalk [14].(c)Protein-protein interactions. Both the Human Protein Reference Database (HPRD) [15] and the STRING database [16] document experimentally verified protein-protein interactions; besides, STRING also provides putative protein-protein interactions ranked by confidence scores. Further details about the exact mechanism of protein-protein interactions can be found in Reactome [17].

( 2)  Network Construction and Visualization. Based on the information collected, a gene regulatory network is constructed and visualized for providing an overview. For this purpose, we recommend CellDesigner which uses standardized symbols (Systems Biology Graphical Notation—SBGN) [18] for visualization and stores gene regulatory networks in the SBML format (Systems Biology Markup Language) [19]. CellDesigner also provides the possibility to simulate temporal dynamics of the gene regulatory network due to the integration of the SBML ODE (ordinary differential equation) solver. Besides, Cytoscape is another powerful tool for integration of biological networks and gene expression data [20].

For assessing the reliability of interactions considered in gene regulatory networks, confidence scores can be computed as being documented in our previous publication [4]. The factors that are used to determine the confidence score for molecular interactions can be: the number of publications reporting an interaction, experimental methods used to identify an interaction, interaction types and computational predictions. The computed confidence scores range from 0 to 1, where values towards 1 indicate higher confidence, whereas values towards 0 indicate lower confidence in a given interaction. For example, the confidence score for a miRNA:gene interaction can be calculated using the following equation:

where are weights that are assigned to the scores which account for the number of publications (), detection method () and the number of predicted binding sites (). The values of the weights can be assigned based on expert knowledge, and the higher the value of the weight is, the bigger impact it has on the confidence score for the interaction. The values of and can be calculated using the logarithmic equation: , where denotes the number of publications describing the miRNA:gene interaction or the number of binding sites that the miRNA has in the 3′ UTR (untranslated region) of the gene. The value of is a cut-off that represents the number of publications or binding sites required for to obtaining their maximum values. Various methods such as western blots, qRT-PCR and reporter assays can be applied to support the miRNA:gene interaction, but these methods provide experimentalists with different levels of confidence, thus differing confidences can be reflected using different values for based on the experience of experimentalists. Of note, although the confidence scores cannot be directly converted into a mathematical model, with the help of these scores we can discard non-reliable putative interactions to generate the ultimate version of a gene regulatory network. The final version of the network can be further analyzed to identify regulatory motifs like feedforward loops (FFLs), for example, with the help of the Cytoscape plugin NetDS [21]. Thereafter, the complete network or parts of it can be converted into a detailed mechanistic model which is described in detail in the following section.

( 3)  Model Construction and Calibration. After the construction, visualization and refinement of a gene regulatory network, it is converted into a detailed mechanistic model which enables the investigation of unanswered biological questions and validation of hypotheses. Ordinary differential equations (ODEs) describe how processes of synthesis, biochemical modification and/or degradation affect the temporal concentration profile of biochemical species like proteins, RNAs and metabolites. An ODE model can be constructed using appropriate kinetic laws such as the law of mass-action, which states that the rate of a chemical reaction is proportional to the probability that the reacting species collide. This collision probability is in turn proportional to the concentration of the reactants [22]. A general representation of ODE-based models using mass-action kinetics is given by the following equation

where represents state variables which denote the molar concentration of the th biochemical specie. Every biochemical reactionis described as a product of a rate constant () and biochemical species that are involved in this reaction. , the so-called stoichiometric coefficients, relate the number of reactant molecules consumed to the number of product molecules generated in the reaction . denotes kinetic orders which are equal to the number of species of involved in the biochemical reactionThe rate constants, kinetic orders and the initial conditions of state variables are defined as model parameters. Besides mass-action kinetics, other kinetic rate laws such as Michaelis-Menten kinetics, Hill equation and power-laws are also frequently used in ODE models.

After ODEs are formulated, the model requires to be calibrated, a process by which parameter values are adjusted in order to make model simulations match experimental observations as good as possible. To do so, there are two possible means: characterization of parameter values using available biological information or estimation of parameter values using optimization methods. Some parameter values can be directly measured or obtained from literature or databases. For example, the half-life () of some molecules (e.g., protein) can be measured in vitro via western blotting. This information can be used to characterize their degradation rate constants through the equation . The database SABIO-RK provides a platform for modelers of biochemical networks to assemble information about reactions and kinetic constants [27]. However, for most model parameters, whose values cannot be measured in laboratories or be accessed from literature or databases, parameter estimation is a necessary process to characterize their values. Before running parameter estimation, initial parameter values and boundaries should be set within physically plausible ranges. To do so, the database BioNumbers provides modelers with key numbers in molecular and cell biology, ranging from cell sizes to metabolite concentrations, from reaction rates to generation times, from genome sizes to the number of mitochondria in a cell [28]. After parameter estimation, unknown parameter values are determined using optimization methods which can minimize a cost function that measures the goodness of fit of the model with respect to given quantitative experimental data sets. Parameter estimation using optimization algorithms is an open research field, in which several methods have been developed according to the nature and numerical properties of biological data analyzed. The discussion for choosing proper optimization methods for parameter estimation is beyond the scope of this paper, but the interested reader is referred to the paper published by Chou and Voit [29].

( 4)  Model Validation and Analysis. Usually, the model simulations are compared with the experimental data used for the parameter estimation, but a good agreement between both is not enough to guarantee the predictive ability of the model. Therefore, it is necessary to validate the model with data sets that are not used during the parameter estimation. This process is called model validation and can ensure more reliable and accurate model predictions. To do so, the data generated in new experiments or extracted from literature are compared with model simulations, which are obtained after configuring the model according to the new experimental settings. Once a model is validated, it can be used to perform predictive simulations, which are helpful to study the dynamic properties of biochemical systems, guide the design of new experiments in the laboratory and formulate additional hypotheses. In addition, tools such as sensitivity and bifurcation analysis can be used to study complex properties and behavior of the modeled system. Sensitivity analysis is used to evaluate the influence of model parameters (e.g., initial concentrations of the state variables and rate constants) on model outputs, such as the temporal behavior of network components [30]. Bifurcation analysis is employed to detect control parameters (also known as bifurcation parameters) whose variations are able to change drastically the dynamical properties of the biochemical system, as well as the stability of its fixed points [31]. The application of these tools to mathematical modeling is beyond the scope of this paper, but the interested reader is referred to the publication of Zhou et al. [32] and Marino et al. [33].

2.1.2. Experiment Methods

As mentioned in the previous sections, after a model is established, it can be calibrated and validated using temporal experimental data. To do so, the data can either be derived from literature or generated to calibrate the model by own experiments. In case of ODE models, the most suitable data for model calibration is quantitative time-series data obtained from perturbation or quantitative dose-response experiments. The experiments, in which time-series data are measured for different regulators (such as miRNAs and TFs) of a gene regulatory network, can be obtained using the techniques described in the subsections below.

( 1)  qRT-PCR. Quantitative real time polymerase chain reaction (qRT-PCR) has been used to identify mRNAs regulated by overexpression or silencing of a specific miRNA [34, 35]. miRNAs typically exhibit their regulatory effects by associating with specific 3′ UTR regions of the mRNAs called miRNA seed regions [34]. This association can lead either to a temporary inhibition of translation or complete degradation of the mRNA in which case qRT-PCR mediated detection is beneficial. For this, the cells are transfected with the miRNA and non-targeting control oligonucleotides at an appropriate concentration using either a lipid based transfection reagent or nucleofection. The transfected cells are then incubated for the necessary time periods (e.g., 24 h, 48 h, 72 h, etc.) after which lysates are prepared and total RNA is extracted. 1000–2000 ng of the RNA is then converted into cDNA using a reverse transcription kit. Taqman qRT-PCR is performed using 10–20 ng of cDNA and primers labeled with fluorescence probes to detect the transcriptional levels of the target mRNA. Housekeeping genes like GAPDH and HPRT are used as endogenous controls for data normalization. Relative expression at different time points is determined by comparing the Ct values of the miRNA transfected cells with Ct values of non-targeting control transfected cells and expressed as relative expression [36]. However, qRT-PCR in spite of being highly sensitive is applicable specifically under conditions of complete or partial degradation of the target mRNA. miRNA-mediated inhibition in translation can be better demonstrated by techniques like immunoblotting.

( 2)  Western Blot/Immunblotting. Western blotting or protein immunoblotting is a technique to detect the expression of a gene at protein level. This technique is particularly useful in determining the regulatory effects of a miRNA on expression of a target gene which is temporarily inhibited. For this the cells are transfected with a miRNA or an antagomiR as mentioned above and cell lysates are prepared at appropriate time points using protein lysis buffers (e.g., Radio-Immunoprecipitation Assay buffer or RIPA) containing lysis agents like Dithiothreitol (DTT) and protease inhibitors. The proteins from each sample are quantitated using Bradford or BCA reagents and compared with bovine serum albumin (BSA) standards for accurate protein estimation. 20–40 μg of protein is then loaded and resolved on a sodium dodecyl sulfate polyacrylamide gel (SDS-PAGE) along with a pre-stained protein marker. The protein bands are then transferred onto a nitrocellulose membrane followed by incubation with the appropriate primary and secondary antibodies linked to fluorescent dyes or horse radish peroxidase enzyme (HRP). The protein expression is then analyzed using either fluorescence or chemiluminiscence HRP substrates in gel documentation systems (e.g., LI-COR Odyssey). Housekeeping genes like β-actin or β-tubulin are used for protein normalization. The time point for maximum target gene suppression generally varies depending on the number of miRNA binding sites at the 3′ UTR of the target gene and the extent of complementarity of the seed region [37]. Immunoblotting is a widely used technique to provide confirmatory evidence for the inhibitory effects of miRNA at the protein level, but it fails to explain the underlying interaction mechanisms.

( 3)  Reporter Gene Assay. As each miRNA can inhibit the expression of a large number of genes, regulation of a particular target gene may either be by direct interaction or be an indirect consequence of it. In direct regulation, a miRNA binds to the complementary sequences at the 3′ UTR of a target gene and thereby suppresses its expression. As a consequence of this, the expression levels of a number of downstream genes (indirect targets) are also dysregulated making it crucial to differentiate between primary and secondary miRNA targets. To determine the interaction specificity, a reporter construct (luciferase) with intact or mutated 3′ UTR of the target gene cloned at the 5′ end is co-transfected into the cells along with the miRNA. The regulatory effect of the miRNA on the target gene expression is then measured using the expression of a reporter gene. In the absence of the appropriate binding sequences (mutated 3′ UTR), the miRNA cannot suppress the reporter mRNA suggesting that the suppressive effect of the miRNA is mediated by a direct interaction. The reporter activity can be analyzed at different time points such as 24 hr, 48 hr and 72 hr to determine the time dependent suppression of a target gene expression by a miRNA.

2.2. Case Study: The Regulation of p21 by Multiple and Cooperative miRNAs

2.2.1. Construction and Visualization of p21 Regulatory Network

By using the approach described above, we investigated the regulation of p21 by its multiple targeting miRNAs. p21, also known as cyclin-dependent kinase inhibitor 1 (CDKN1A), is a transcriptional target of p53. It is required for proper cell cycle progression and plays a role in cell death, DNA repair, senescence and aging (reviewed in [38]). Interestingly, p21 was the first experimentally validated miRNA target hub, which is a gene that is simultaneously regulated by many miRNAs. This made it an ideal candidate for a case study of our approach [23]. To do so, we first constructed a p21 regulatory network using the following steps:

(a)We extracted miRNA-target interactions from the publication of Wu et al. [23] where a list of predicted p21-regulating miRNAs was subjected to experimental validation.(b)Experimentally verified TFs of p21 were extracted from literature and putative TFs having conserved binding sites in the 5 kb upstream region of the p21 open reading frame were extracted from UCSC table browser. A list of TFs controlling the expression of the miRNAs was constructed using information of experimentally proven TF-miRNA interactions extracted from TransmiR (release 1.0) [39]. In addition, we generated a list of putative TFs of miRNAs with binding sites in the 10 kb upstream region of the miRNA genes using information from the databases PuTmiR (release 1.0) [40] and MIR@NT@N (version 1.2.1) [41], and from the table of TFs with conserved binding sites in the UCSC genome browser (hg18) [9]. (c)Information about protein interactions was extracted from the Human Protein Reference Database (HPRD, release 9.0) [15] and the STRING database (release 9.0) [16]. Only the experimentally verified p21-protein interactions were used to construct the network.(d)Additionally, we associated the TFs in the network to nine biological processes based on the Gene Ontology (GO) [42]. The corresponding GO terms were cell proliferation, cell apoptosis, immune response, inflammatory response, cell cycle, DNA damage, cell senescence, DNA repair and cell migration.

Next, we visualized the network in CellDesigner and computed a confidence score for each interaction in the network (Figure 1). The confidence scores provide us with the reliability of the interactions considered in p21 regulatory network. With the help of these scores, we discarded non-reliable interactions and constructed the mechanistic model accounting for p21 regulation by its targeting miRNAs. Besides, the interested experimentalists can further use this information to choose reliable interacting candidates of p21 for designing relevant experiments. The scores for each interaction of p21 regulatory network are shown in Table 2.

Table 2: The confidence scores for the interactions in p21 regulatory network. The network consists of four types of interactions: miRNA-p21, TF-miRNA, TF-p21, and p21-protein.

(a) miRNA-p21 interaction scores

miRNA
miR-125a-5p
miR-132
miR-208a
miR-28-5p
miR-298
miR-299-5p
miR-345
miR-363

Score
0.73
0.73
0.73
0.80
0.73
0.80
0.73
0.73

miRNA
miR-423-3p
miR-515-3p
miR-572
miR-639
miR-654-3p
miR-657
miR-93

Score
0.85
0.73
0.73
0.73
0.73
0.80
0.95

(b) TF-miRNA interaction scores

Score

Score

Score

Score

Score

Score

miR-208

miR-132

miR-657

miR-125a-5p

miR-28-5p

miR-93

    MAX
0.27
    EGR1
0.93
    TFAP2A
0.34
    EGR1
0.86
    MZF1
0.28
    TFAP2A
0.34

    EGR1
0.29
    CREB1
0.86
    SP1
0.34
    PAX5
0.24
    E47
0.27
    SOX9
0.24

    ARNT
0.27
    AREB6
0.27
    MZF1
0.34
    NRSF
0.24
    FOXC1
0.31
    GATA3
0.27

    RFX1
0.24
    E47
0.24
    MAX
0.29
    ZID
0.24
    FOXI1
0.24
    MZF1
0.31

    ATF6
0.24
    ELK1
0.27
    USF1
0.27
    PPARG
0.27
    SRY
0.27
    SP1
0.34

    YY1
0.27
    EGR2
0.29
    EGR1
0.34
    GATA3
0.34
    JUN
0.29
    POU2F1
0.24

    JUNB
0.24
    FOS
0.24
    RELA
0.29
    MZF1
0.34
    POU2F1
0.24
    EGR1
0.29

    POU2F1
0.24
    NFIC
0.24
miR-654-3p

    TP53
0.34
    SRF
0.27
    CEBPA
0.24

    NFIC
0.27
    RELA
0.27
    YY1
0.27
    NFIC
0.34
    CEBPA
0.29
    NFYA
0.24

miR-345

miR-423-3p

    BACH1
0.24
miR-299-5p

    RORA
0.24
    RUNX1
0.24

    SP1
0.34
    BACH1
0.27
    FOXL1
0.29
    CREB1
0.24
miR-363

    STAT1
0.24

    RELA
0.27
    STAT1
0.27
    SRY
0.24
    SRY
0.27
    E2F1
0.62
    E2F1
0.85

    NFATC2
0.24
    POU3F2
0.24
    NFATC2
0.24
    SRF
0.29
    FOXC1
0.33
    MYC
0.84

    HNF4A
0.31
    STAT5A
0.27
    FOS
0.24
    RUNX1
0.27
    MZF1
0.34
miR-298

    NFYA
0.24
    SRF
0.27
    POU2F1
0.27
miR-572

miR-639

    JUNB
0.31

    POU2F1
0.24
    MEF2A
0.27
    NFB1
0.24
    FOXF2
0.24
    NFYA
0.24

    NFIC
0.31
    POU2F1
0.24
    HNF1A
0.24

(c) TF-p21 interaction scores

TF (verified)
SP1
SP3
E2F1
Runx1
Runx2
STAT1
E2A
STAT5
TP53
TP63
TP73

Score
1.00
1.00
0.77
0.77
0.77
0.72
0.72
0.67
0.67
0.67
0.67

TF (verified)
STAT3
CUX1
TFAP2A
BRCA1
VDR
RARA
C/EBPα
C/EBPβ
RB1
Tbx2
EHMT2

Score
0.66
0.66
0.60
0.60
0.60
0.60
0.60
0.60
0.60
0.60
0.60

TF (putative)
NFB1
RELA
STAT2
STAT6
SRF

Score
0.44
0.44
0.44
0.44
0.43

(d) p21-protein interaction scores

Protein
TP53
PCNA
CASP3
CCNA1
CCND1
SKP2
BCCIP
CCNA2
CCND2
CCNE2
AKT1
C1orf123
CCDC85B
CCNB1
CCNB2
CCND3

Score
1.00
0.92
0.87
0.87
0.87
0.87
0.81
0.81
0.81
0.81
0.72
0.72
0.72
0.72
0.72
0.72

Protein
CCNE1
CDC45
CDC5L
CDC6
CDC7
CDK1
CDK14
CDK2
CDK3
CDK4
CDK6
CEBPA
CIZ1
CSNK2A1
CSNK2B
DAPK3

Score
0.72
0.72
0.72
0.72
0.72
0.72
0.72
0.72
0.72
0.72
0.72
0.72
0.72
0.72
0.72
0.72

Protein
ESR1
GADD45A
GADD45B
GADD45G
GMNN
GNB2L1
HDAC11
ITGB1BP3
MAP3K5
MAPK8
MCM10
PARP1
PIM1
POLD2
PSMA3
RAB1A

Score
0.72
0.72
0.72
0.72
0.72
0.72
0.72
0.72
0.72
0.72
0.72
0.72
0.72
0.72
0.72
0.72

Protein
SET
SLC25A11
STAT3
TEX11
TK1
TSG101
TTLL5
XRCC6
CDC20
CDC27
CDK5
DHX9
MDM2
NR4A1
PRKCH
RPA1

Score
0.72
0.72
0.72
0.72
0.72
0.72
0.72
0.72
0.56
0.56
0.56
0.56
0.56
0.56
0.56
0.56

Verified: experimentally verified interaction; putative: predicted interaction.

Figure 1: p21 regulatory network. The network contains several layers of regulators of p21: TFs (light blue and red boxes), miRNAs (dark blue and green boxes), and proteins (grey boxes). In each big box, there are small boxes which represent individual components of this layer of regulation. miRNAs are classified into two groups according to the mechanisms by which the expression of p21 is repressed. One group causes p21 translation repression (green box). These miRNAs bind to p21 mRNA resulting in the repressed translation of p21 but unchanged mRNA expression level. The other group of miRNAs (dark blue box) decreases the stability of p21 mRNA by modifying its structure, leading to mRNA decay and finally the downregulation of p21. TFs are classified into three groups: p21 TFs (red), miRNA TFs (yellow) and their common TFs (light blue). The p21 interacting-proteins are framed in the grey boxes. The purple boxes represent nine processes, and the TFs associated with these processes are indicated in the ellipses above them using corresponding figures. This data is adapted from our previous publication [4].

2.2.2. Mechanistic Modeling of p21 Regulation by Its Targeting miRNAs

( 1)  Model Construction. After constructing the regulatory network, a detailed mechanistic model of ODEs, which describes the biochemical reactions underlying the regulation of p21 was established. We chose the ODE modeling approach, as it is a simple formalism for describing temporal dynamics of biochemical systems and a wide range of tools are available to explore their properties. Precisely, the model considered the mRNA (mp21; (3)) and protein (p21; (6)) of the miRNA target hub p21, the p21-targeting miRNAs (;; (5)), and the complexes formed by p21 mRNA and miRNA, (4). Altogether, the model is constituted by 32 state variables and 64 parameters:

For mp21, processes considered in the model were: (i) its synthesis () mediated by TFs (), (ii) its degradation (), and (iii) its association with a miRNA (). For each , processes considered were: (i) the synthesis () mediated by TFs (), (ii) the degradation (), and (iii) the association with the p21 mRNA target (). For each complex, processes considered were: (i) the formation of the complex by a and the p21 mRNA (), and (ii) the complex degradation (). For p21, processes considered were: (i) its synthesis (, and (ii) its degradation (). An additional algebraic equation accounting for the total measurable amount of p21 mRNA () was also included. The SBML file of the model is available for download at http://www.sbi.uni-rostock.de/uploads/tx_templavoila/p21_TargetHub_03092013.xml.

( 2)  Model Calibration and Validation. For model calibration, we fixed some parameter values using published data and estimated the other unknown parameters using the time-series data published by Wu et al. [23], in which the p21 mRNA (northern plot) and protein levels (western plot) were measured 48 hr after transfection of individual p21-targeting miRNAs into human embryonic kidney 293 cells. The unknown parameter values were estimated using an iterative method combining global (particle swarm pattern search) [43] and local (downhill simplex method in multi-dimensions) [44] optimization algorithms. For each considered in the model, the method minimizes the distance between model simulations and experimental data using the following cost function

whereandrepresent the simulated p21 mRNA and protein expression levels for each at time point .andrepresent the measured value for each at time point , and their standard deviations areand. Here, is the time point (48 hr) after overexpression of the individual miRNAs in embryonic kidney 293 cells at which the expression levels of the p21 and its mRNA were measured [23]. The model calibration results are shown in Figure 2(a) and the obtained parameter values are listed in Table 3.

Table 3: Initial concentrations of model variables and model parameter values. Based on the experimental data, the p21-targeting miRNAs verified by Wu et al. [23] were divided into two groups: the translation repression group (marked with asterisk) and the mRNA deadenylation group. A miRNA was classified into the mRNA deadenylation group if its overexpression can result in 20% or more downregulation of the p21 mRNA level (i.e., p21 mRNA level ≤ 0.8; the basal level is 1); otherwise, it was classified into the translation repression group. For the translation repression group, only was estimated and was fixed. For the other group, both and were estimated. The initial concentrations of p21 and mp21 were set to 1, and this value was used as their basal expression levels. During the parameter estimation, the initial concentrations of p21-targeting miRNAs were set to 100, because in the publication [23] the expression levels of p21  and mp21 were measured after the individual introduction of the p21-targeting miRNAs with amount of 100 nM. Due to the lack of biological information to characterize the transcriptional activation function () of p21 and its targeting miRNAs, the corresponding functions were assumed to be 1 for simplicity. The data is adapted from our previous publication [4].

Initial concentration of variables and TF functions

Variable
Description
Initial concentration (a.u.)

p21
p21 protein
1

mp21
p21 mRNA
1

p21-targeting miRNAs
100

Complexes formed by and mp21
0

()
p21's transcriptional activation function
1

The transcriptional activation function of
1

Fixed parameter values

Parameter
Description
Value (hr⁻¹)
Reference

Synthesis rate constant of mp21
0.1155
fixed

Degradation rate constant of mp21
0.1155
[24]

()
Synthesis rate constant of
0.0289
fixed

(5)
Degradation rate constant of
0.0289
[25]

Synthesis rate constant of p21
1.3863
fixed

Degradation rate constant of p21
1.3863
[26]

Estimated parameter values

miRNA (state variable)

(hr⁻¹)

(a.u. ⁻¹· hr⁻¹)

Experimental data of p21 (protein, mRNA ± SD)

miR-298 ()*
0.1155
0.0254

miR-208a ()*
0.1155
0.0041

miR-132 ()*
0.1155
0.0275

miR-28-5p ()*
0.1155
0.0119

miR-125-5p ()*
0.1155
0.0018

miR-299-5p ()*
0.1155
0.0080

miR-345 ()*
0.1155
0.0051

miR-93 ()
0.1564
0.0235

miR-423-3p ()
0.9118
0.0055

miR-515-3p ()
0.2098
0.0253

miR-363 ()
0.2261
0.0399

mR-657 ()
0.3465
0.0158

miR-639 ()
0.4305
0.0327

miR-572 ()
0.3039
0.0360

miR-654-3p ()
9.7485
0.0024

(a)

(b)

(c)

(d)

Figure 2: Model calibration and validation. (a) Model calibration. The figures show the relative change of the p21 mRNA and protein expression levels after overexpression of the indicated miRNAs (Model: model simulation; Data: experimental data). These data were normalized to the control group in which the p21 mRNA and protein expression levels were measured when the miRNAs were normally expressed (a.u.: arbitrary unit). (b) Experimental workflow. In the experiments, Sk-Mel-147 cells were seeded in six well plates. Then, mature miRNA mimics were transfected individually at a concentration of 100 nM (miR-572 and miR-93) or in combination at 50 nM each (miR-572 + miR-93). After 48 hr transfection with miRNA mimics, the cells were pulse treated with 250 nM doxorubicin for 1 hour after which normal growth medium was replenished. The immunblotting were performed to measure p21 expression at 0, 2, 4, 6, 8 and 24 hr post-doxorubicin treatment. (c) Temporal dynamics of p21 transcriptional function. Afterdoxorubicin treatment, the expression of p53, a TF of p21, was measured using immunblotting and these data were used to characterize the transcriptional function of p21 using MATLAB linear interpolation function. (d) Model validation. We measured the expression of p21 protein in response to genotoxic stress in the four scenarios as described in the main text. The measured data (Data) were compared with the model simulations (Model). The figures (a), (c) and (d) are adapted from our previous publication [4].

Experimental results showed that a stronger repression of the target gene can occur when two miRNA binding sites on the target mRNA are in close proximity [24, 45]. To test the consequences of this hypothesis, we predicted cooperative miRNA pairs for p21, with seed site distances between 13–35 nt in the p21 3′ UTR. To substantiate the cooperative effect associated with pairs of miRNAs, we introduced a group of new state variables () into the original model. These state variables account for the ternary complexes composed of p21 mRNA and two putatively cooperating miRNAs ( and ). For these new variables, processes considered are: (i) the association of p21 mRNA with and into a complex (), and (ii) the degradation of the complex (). After expansion, the corresponding modified and new ODEs are listed below:

To model stronger repression of the target gene by cooperating miRNAs, we assumed a stronger association rate constant for the complex which is equal to the sum of their individual association rate constants (). Similarly, the degradation rates of the complexes were assumed to be equal to the sum of degradation rate constants of single miRNA binding complexes (). However, it has to be noted that these added equations are an abstract description of miRNA cooperativity, because the details of this mechanism are not yet known.

To experimentally validate the capability of our model to predict the relative p21 concentrations regulated by cooperative miRNAs, we selected miR-572 and miR-93 as a case study. These two miRNAs were chosen, because their predicted target sites in the p21 3′ UTR are in close proximity to each other and thereby, they can induce cooperative repression on p21 as suggested in [24]. The experiments were performed as follows:

(i)Melanoma cells (Sk-Mel-147) were transfected with the mature miRNA mimics of the two miRNAs either individually (100 nM) or in combination (50 nM each), whereas untreated cells were used as control (Figure 2(b)).(ii)Next, the cells were treated with doxorubicin, a genotoxic-stress inducing agent. The agent can upregulate the expression of p53, which is a known TF of p21, and therefore it can result in the upregulation of p21 (Figure 2(c)).(iii)After doxorubicin treatment, the expression levels of p21 were measured by immunoblotting at different time points (0, 2, 4, 6, 8, 24 hr). The p21 expression values were normalized based on the p21 expression level in the control group measured at time point 0 hr (Figure 2(d)).

By doing so, we obtained the p21 response after genotoxic stress in four scenarios: (1) endogenous miRNA expression (Control); (2) overexpression of miR-572; (3) overexpression of miR-93; and (4) both miRNAs moderately overexpressed. Thereafter, we derived a model of seven ODEs based on the original equations which was configured according to the designed experiments, making the simulation results comparable with the experimental data. As shown in the Figure 2(c), the simulations are in good agreement with the experimental observations. The individual overexpression of miR-572 or miR-93 led to the reduction of the upregulation of p21 response after genotoxic stress induction. The two miRNAs cause different degrees of repression due to their different repression efficiencies on p21. Interestingly, the combined overexpression of both miRNAs induced the strongest downregulation of p21, and therefore verifying the hypothesis of their cooperative regulation of p21. Above all, the results not only validated the model but also demonstrated the ability of our method to identify cooperative miRNA pairs for p21.

( 3)  Predictive Simulations. As there are abundant of network motifs such as FFLs in p21 regulatory network and these network motifs are important for determining p21 dynamics, it is interesting to investigate the dynamics of p21 in network modules where FFLs are involved. To do so, two network modules including both miR-93 and miR-572, and their TFs were exemplified. In Figure 3(a), the network module contains an incoherent FFL composed by AF2α, miR-93 and p21, and a cascaded regulation in which p21 is repressed by FOXF2 via miR-572; in Figure 3(b), the same cascaded regulation together with a coherent FFL composed of TP53, MYC, miR-93 and p21 forms another regulatory module of p21. By modulating the transcriptional strengths of the two miRNAs by their TFs, two types of p21 dynamics were identified: saturation and pulse (Figure 3(c), top). In the former, the p21 expression increases and reaches its steady state at the highest level; in contrast, in the latter the p21 expression increases to a peak and thereafter drops to a steady state at lower level. For the two different network modules, various combinations of transcriptional strengths of the two miRNAs lead to different distributions of the two p21 dynamical patterns. For the network module of an incoherent FFL plus the cascaded regulation (Figure 3(c), bottom left), the saturation pattern appears in two distinct regions: one with weakened transcriptional strength of miR-93 and the other with enhanced transcriptional strength of miR-93; for the other module (Figure 3(c), bottom right), the saturation pattern only appears in the region, in which the transcriptional strength of miR-93 is enhanced. Taken together, the results showed that for the two different network motifs the dynamical pattern of p21 is changing according to different combinations of its upstream regulators, suggesting the adaptation of p21 dynamics for different biological contexts.

(a)

(b)

(c)

Figure 3: Different p21 dynamics in different network motifs. We ran simulations to show the different dynamics of p21 for two different network motifs (a) and (b). Through simulations, two dynamical patterns of p21 were identified: saturation and pulse ((c), top). For each network motif, the corresponding distributions of the two dynamical patterns were plotted ((c), bottom). For different combinations of the transcriptional strengths, the normalized distance () between peaks () and steady states () of p21 is determined by the equation . If , for the corresponding combination of transcriptional strengths the p21 dynamics is saturation, otherwise it is pulse. The regions showing different dynamical patterns of p21 are separated using the white lines.

Furthermore, we performed a number of simulations to show the influence of miRNA regulation on p21 expression levels in different cellular processes arranged in a consecutive manner. In this procedure, a cell is first in the process of cell cycle followed by proliferation (cell growth), then the cell responses to DNA damage and enters into the process of senescence, and finally apoptosis is initiated. As shown in Figure 4, during the cell cycle process, the p21 expression is low due to the activation of most of its targeting miRNAs; when the cell starts proliferating, the p21 expression declines to an even lower level because of more activated miRNAs in this process; after responding to DNA damage, the p21 expression soars to a high level, which is caused by the activation of its TFs like p53 and fewer expressed miRNAs.; although the p21 expression keeps at a similar level while the cell is undergoing senescence, the expressed miRNAs are different from the previous process; finally, the p21 expression decreases again to a low level due to the reemergence of most of its targeting miRNAs and the cell enters apoptosis. Interestingly, the model simulations are also consistent with experimental observations: under non-stressed condition the low expression level of p21 is needed for cell proliferation; the upregulation of p21 happens after response to DNA damage via p53 and the increased p21 expression further results in cell cycle arrest leading to senescence and apoptosis [25]. Besides, when considering the effect of cooperating miRNAs, the p21 expression levels were indistinguishable from the previous simulation of DNA damage response. However, for the other processes p21 expression levels were significantly lower compared to the simulations without considering the cooperative effect of the miRNAs. Above all, these results indicated that selective expression of cooperative miRNAs could be adopted by cells to ensure diverse expression levels of p21 to meet the requirements of different cellular processes.

Figure 4: p21 expression regulated by cooperative miRNAs for different cellular processes. The associations of the miRNAs with these cellular processes were derived from GO terms of their TFs. A miRNA was supposed to be expressed (in bold black font) in a cellular process only if its TF is related to the corresponding GO term of this process. The p21 expression levels are computed for each process with (w)/without (wo) considering the cooperative effect among the p21 targeting miRNAs.

3. Conclusions and Discussion

In this paper, we presented a systems biology approach, combining data-driven modeling and model-driven experiments, to investigate the role of miRNA-mediated repression in gene regulatory networks. This approach provides a systematic way to gain a deeper understanding of the regulation of target genes by mutiple and cooperative miRNAs. Using the regulation of p21 by multiple miRNAs as a case study, we showed how the ODE-based model, which is calibrated and validated by means of experimental data, is suitable for predicting the temporal dynamics of molecular concentrations involved in biochemical systems.

Provided there are sufficiently rich quantitative data sets avaliable to characterize the model, the use of the methodology here shown can be extended to more complex regulatory networks, involving multiple targets, cooperating TFs and miRNAs and signaling pathways displaying cross-talk via post-tranlational modifications. In this case, the critical element is the quality and quantitiy of the available data. Insufficiency and low quality of experimental data can cause errors in the process of model construction and overfitting in parameter estimation can lead to uncertainties in the model predictions. We believe that the quick development of quantitative high throughput techniques such as transcriptomics, proteomics and miRNomics will facilitate the construction and characterization of larger miRNA-mediated regulatory networks.

Other modeling frameworks than ODE-based models can be used to describe biological systems, such as probabilistic (e.g., Bayesian) or logical (e.g., Boolean) models. Importantly, different modeling frameworks have different properties and perform well regarding different perspectives and levels of mechanistic details of biochemical systems [46]. For example, Bayesian models are helpful in the construction of connections in signaling networks and can reveal the most likely underlying structure of the network in a probabilistic manner. Boolean models use binary values (0 and 1) and logical gates (AND, NOT, and OR) to describe activities of network components and the information flow among them. We believe that in the coming future, hybrid models, which consist of modeling framework and experimental technique specific sub-modules, will provide the necessary compromise between quantitative/qualitative accuracy and scalability for the investigation of large biochemical networks [47].

Disclosures

The authors declare that they have no competing financial interests.

Authors’ Contribution

Julio Vera and Olaf Wolkenhauer are equal contributers.

Acknowledgments

The authors would like to acknowledge the funding from the German Federal Ministry of Education and Research (BMBF): eBio-miRSys (0316175A to Xin Lai and Julio Vera), eBio-SysMet (0316171 to Olaf Wolkenhauer) and Gerontosys-ROSAge (0315892A to Olaf Wolkenhauer).

References

A. E. Pasquinelli, “MicroRNAs and their targets: recognition, regulation and an emerging reciprocal relationship,” Nature Reviews Genetics, vol. 13, no. 4, pp. 271–282, 2012.

J. Vera, X. Lai, U. Schmitz, and O. Wolkenhauer, “MicroRNA-regulated networks: the perfect storm for classical molecular biology, the ideal scenario for systems biology,” Advances in Experimental Medicine and Biology, vol. 774, pp. 55–76, 2013.

S. Nikolov, J. Vera, U. Schmitz, and O. Wolkenhauer, “A model-based strategy to investigate the role of microRNA regulation in cancer signalling networks,” Theory in Biosciences, vol. 130, no. 1, pp. 55–69, 2011.

X. Lai, U. Schmitz, S. K. Gupta et al., “Computational analysis of target hub gene repression regulated by multiple and cooperative miRNAs,” Nucleic Acids Research, vol. 40, no. 18, pp. 8818–8834, 2012.

X. Lai, O. Wolkenhauer, and J. Vera, “Modeling miRNA regulation in cancer signaling systems: miR-34a regulation of the p53/Sirt1 signaling module,” Methods in Molecular Biology, vol. 880, pp. 87–108, 2012.

C. Jiang, Z. Xuan, F. Zhao, and M. Q. Zhang, “TRED: a transcriptional regulatory element database, new entries and other development,” Nucleic Acids Research, vol. 35, supplement 1, pp. D137–D140, 2007.

L. A. Bovolenta, M. L. Acencio, and N. Lemke, “HTRIdb: an open-access database for experimentally verified human transcriptional regulation interactions,” BMC Genomics, vol. 13, article 405, 2012.

V. Matys, O. V. Kel-Margoulis, E. Fricke et al., “TRANSFAC and its module TRANSCompel: transcriptional gene regulation in eukaryotes,” Nucleic Acids Research, vol. 34, database issue, pp. D108–D110, 2006.

D. Karolchik, R. Baertsch, M. Diekhans et al., “The UCSC genome browser database,” Nucleic Acids Research, vol. 31, no. 1, pp. 51–54, 2003.

P. Alexiou, T. Vergoulis, M. Gleditzsch et al., “miRGen 2.0: a database of microRNA genomic information and regulation,” Nucleic Acids Research, vol. 38, supplement 1, pp. D137–D141, 2010.

F. Xiao, Z. Zuo, G. Cai, S. Kang, X. Gao, and T. Li, “miRecords: an integrated resource for microRNA-target interactions,” Nucleic Acids Research, vol. 37, supplement 1, pp. D105–D110, 2009.

S.-D. Hsu, F.-M. Lin, W.-Y. Wu et al., “miRTarBase: a database curates experimentally validated microRNA-target interactions,” Nucleic Acids Research, vol. 39, supplement 1, pp. D163–D169, 2011.

P. Sethupathy, B. Corda, and A. G. Hatzigeorgiou, “TarBase: a comprehensive database of experimentally supported animal microRNA targets,” RNA, vol. 12, no. 2, pp. 192–197, 2006.

H. Dweep, C. Sticht, P. Pandey, and N. Gretz, “miRWalk—database: prediction of possible miRNA binding sites by “walking” the genes of three genomes,” Journal of Biomedical Informatics, vol. 44, no. 5, pp. 839–847, 2011.

T. S. Keshava Prasad, R. Goel, K. Kandasamy et al., “Human protein reference database—2009 update,” Nucleic Acids Research, vol. 37, supplement 1, pp. D767–D772, 2009.

D. Szklarczyk, A. Franceschini, M. Kuhn et al., “The STRING database in 2011: functional interaction networks of proteins, globally integrated and scored,” Nucleic Acids Research, vol. 39, supplement 1, pp. D561–D568, 2011.

D. Croft, G. O'Kelly, G. Wu et al., “Reactome: a database of reactions, pathways and biological processes,” Nucleic Acids Research, vol. 39, supplement 1, pp. D691–D697, 2011.

N. le Novère, M. Hucka, H. Mi et al., “The systems biology graphical notation,” Nature Biotechnology, vol. 27, no. 8, pp. 735–741, 2009.

A. Funahashi, Y. Matsuoka, A. Jouraku, M. Morohashi, N. Kikuchi, and H. Kitano, “CellDesigner 3.5: a versatile modeling tool for biochemical networks,” Proceedings of the IEEE, vol. 96, no. 8, pp. 1254–1265, 2008.

M. S. Cline, M. Smoot, E. Cerami et al., “Integration of biological networks and gene expression data using Cytoscape,” Nature Protocols, vol. 2, no. 10, pp. 2366–2382, 2007.

D.-H. Le and Y.-K. Kwon, “NetDS: a Cytoscape plugin to analyze the robustness of dynamics and feedforward/feedback loop structures of biological networks,” Bioinformatics, vol. 27, no. 19, pp. 2767–2768, 2011.

E. O. Voit, Computational Analysis of Biochemical Systems: A Practical Guide for Biochemists and Molecular Biologists, Cambridge University Press, Cambridge, UK, 2000.

S. Wu, S. Huang, J. Ding et al., “Multiple microRNAs modulate p21Cip1/Waf1 expression by directly targeting its 3′ untranslated region,” Oncogene, vol. 29, no. 15, pp. 2302–2308, 2010.

P. Sætrom, B. S. E. Heale, O. Snøve Jr., L. Aagaard, J. Alluin, and J. J. Rossi, “Distance constraints between microRNA target sites dictate efficacy and cooperativity,” Nucleic Acids Research, vol. 35, no. 7, pp. 2333–2342, 2007.

Y.-S. Jung, Y. Qian, and X. Chen, “Examination of the expanding pathways for the regulation of p21 expression and activity,” Cellular Signalling, vol. 22, no. 7, pp. 1003–1012, 2010.

W. Wang, H. Furneaux, H. Cheng et al., “HuR regulates p21 mRNA stabilization by UV light,” Molecular and Cellular Biology, vol. 20, no. 3, pp. 760–769, 2000.

U. Wittig, R. Kania, M. Golebiewski et al., “SABIO-RK—database for biochemical reaction kinetics,” Nucleic Acids Research, vol. 40, database issue, pp. D790–D796, 2012.

R. Milo, P. Jorgensen, U. Moran, G. Weber, and M. Springer, “BioNumbers—the database of key numbers in molecular and cell biology,” Nucleic Acids Research, vol. 38, supplement 1, pp. D750–D753, 2010.

I.-C. Chou and E. O. Voit, “Recent developments in parameter estimation and structure identification of biochemical and genomic systems,” Mathematical Biosciences, vol. 219, no. 2, pp. 57–83, 2009.

A. Saltelli, K. Chan, and E. M. Scott, Sensitivity Analysis, John Wiley & Sons, New York, NY, USA, 1st edition, 2000.

S. Strogatz, Nonlinear Dynamics and Chaos: Applications to Physics, Biology, Chemistry, and Engineering: With Applications to Physics, Biology, Chemistry and Engineering, Westview Press, Boulder, Colo, USA, 2000.

P. Zhou, S. Cai, Z. Liu, and R. Wang, “Mechanisms generating bistability and oscillations in microRNA-mediated motifs,” Physical Review E, vol. 85, no. 4, apart 1, Article ID 041916, 9 pages, 2012.

S. Marino, I. B. Hogue, C. J. Ray, and D. E. Kirschner, “A methodology for performing global uncertainty and sensitivity analysis in systems biology,” Journal of Theoretical Biology, vol. 254, no. 1, pp. 178–196, 2008.

L. P. Lim, N. C. Lau, P. Garrett-Engele et al., “Microarray analysis shows that some microRNAs downregulate large numbers of-target mRNAs,” Nature, vol. 433, no. 7027, pp. 769–773, 2005.

J. Krützfeldt, N. Rajewsky, R. Braich et al., “Silencing of microRNAs in vivo with ‘antagomirs’,” Nature, vol. 438, no. 7068, pp. 685–689, 2005.

M. W. Pfaffl, G. W. Horgan, and L. Dempfle, “Relative expression software tool (REST) for group-wise comparison and statistical analysis of relative expression results in real-time PCR,” Nucleic Acids Research, vol. 30, no. 9, p. e36, 2002.

D. P. Bartel, “MicroRNAs: target recognition and regulatory functions,” Cell, vol. 136, no. 2, pp. 215–233, 2009.

T. Abbas and A. Dutta, “p21 in cancer: intricate networks and multiple activities,” Nature Reviews Cancer, vol. 9, no. 6, pp. 400–414, 2009.

J. Wang, M. Lu, C. Qiu, and Q. Cui, “TransmiR: a transcription factor microRNA regulation database,” Nucleic Acids Research, vol. 38, supplement 1, pp. D119–D122, 2010.

S. Bandyopadhyay and M. Bhattacharyya, “PuTmiR: a database for extracting neighboring transcription factors of human microRNAs,” BMC Bioinformatics, vol. 11, article 190, 2010.

A. Le Béchec, E. Portales-Casamar, G. Vetter et al., “MIR@NT@N: a framework integrating transcription factors, microRNAs and their targets to identify sub-network motifs in a meta-regulation network model,” BMC Bioinformatics, vol. 12, article 67, 2011.

M. Ashburner, C. A. Ball, J. A. Blake et al., “Gene ontology: tool for the unification of biology. The Gene Ontology Consortium,” Nature Genetics, vol. 25, no. 1, pp. 25–29, 2000.

A. I. F. Vaz and L. N. Vicente, “A particle swarm pattern search method for bound constrained global optimization,” Journal of Global Optimization, vol. 39, no. 2, pp. 197–219, 2007.

W. H. Press, Numerical Recipes: The Art of Scientific Computing, Cambridge University Press, Cambridge, UK, 2007.

J. G. Doench and P. A. Sharp, “Specificity of microRNA target selection in translational repression,” Genes and Development, vol. 18, no. 5, pp. 504–511, 2004.

B. Kholodenko, M. B. Yaffe, and W. Kolch, “Computational approaches for analyzing information flow in biological networks,” Science Signaling, vol. 5, no. 220, 2012.

F. M. Khan, U. Schmitz, S. Nikolov et al., “Hybrid modeling of the crosstalk between signaling and transcriptional networks using ordinary differential equations and multi-valued logic,” Biochimica et Biophysica Acta, 2013.