Computational analysis of target hub gene repression regulated by multiple and cooperative miRNAs.

MicroRNA (miRNA) target hubs are genes that can be simultaneously targeted by a comparatively large number of miRNAs, a class of non-coding RNAs that mediate post-transcriptional gene repression. Although the details of target hub regulation remain poorly understood, recent experiments suggest that pairs of miRNAs can cooperate if their binding sites reside in close proximity. To test this and other hypotheses, we established a novel approach to investigate mechanisms of collective miRNA repression. The approach presented here combines miRNA target prediction and transcription factor prediction with data from the literature and databases to generate a regulatory map for a chosen target hub. We then show how a kinetic model can be derived from the regulatory map. To validate our approach, we present a case study for p21, one of the first experimentally proved miRNA target hubs. Our analysis indicates that distinctive expression patterns for miRNAs, some of which interact cooperatively, fine-tune the features of transient and long-term regulation of target genes. With respect to p21, our model successfully predicts its protein levels for nine different cellular functions. In addition, we find that high abundance of miRNAs, in combination with cooperativity, can enhance noise buffering for the transcription of target hubs.

INTRODUCTION

The regulation of basic cellular functions is controlled by complex networks involving interacting genes, proteins and small molecules (1). In the past decade, an additional level of regulation was discovered after the identification of a class of short non-coding RNAs called microRNAs (miRNAs). miRNAs are ∼22 nt in length and their main function is to regulate the activity and stability of specific mRNAs at the post-transcriptional level (1,2). Individual miRNAs can regulate a large number of mRNA targets and target genes can be regulated by multiple miRNAs (3,4).

The term target hub was introduced by Borneman et al. (5), who constructed a putative transcriptional control network for yeast cell differentiation. Genes, combining several transcription factors (TFs) in their promoter regions, were referred to as target hubs. Shortly thereafter, Shalgi et al. (3) defined the notion of miRNA target hub for genes that are regulated by at least 15 different miRNAs (Figure 1).

Figure 1.

Illustration of miRNA target hub regulatory network. Once an miRNA target hub gene is transcribed, its mRNAs are post-transcriptionally regulated by many miRNAs whose expression can be activated or inhibited by numerous TFs. In animals, the complementarity between an miRNA and the target hub mRNA usually results in reduced protein expression through translational repression or mRNA deadenylation followed by degradation. Target hubs are highly interconnected via protein–protein interactions, so that its regulatory network may contain interlocked network motifs, such as feedback and feedforward loops.

At least two groups (6,7) experimentally proved that pairs of miRNAs, with binding sites that reside in close proximity or partially overlap, show cooperative behaviour, adding to the complexity of the miRNA target regulation. Saetrom et al. (7) proposed a range of 13–35 nt for the distance of the binding start sites for which the phenomenon appears. According to these experimental evidences, we assumed this kind of interaction as plausible and applied it in a kinetic model of p21 regulation by multiple miRNAs to investigate the dynamic and regulatory consequences of this feature. A well-known miRNA target hub is the cell cycle regulator and tumour suppressor p21 (a.k.a. CDKN1A or Cip1/Waf1) (8). Specifically, p21 is involved in the G1 phase cell cycle arrest in response to stress signals provoked by DNA damage. The expression of p21 is dependent on environmental conditions and is transcriptionally regulated through p53-dependent and -independent mechanisms (9). Recently, numerous miRNAs were identified to regulate the expression of p21. From 266 predicted miRNA regulators of p21, Wu et al. (10) validated the ability of a subset of 28 miRNAs to target p21 using a high-throughput luciferase reporter screen. Thus, p21, to our knowledge, is one of the first experimentally verified miRNA target hubs. Moreover, the regulation induced by miRNAs can provoke phenotypic changes. Wu et al. (10) showed that overexpression or knockout of miR-372, miR-519-3p and miR-520a-3p in JAR cells leads to either increased cell proliferation or accumulation of cells in G1 phase.

In this work, we investigated the structure and dynamics of the regulatory network of the miRNA target hub p21 and several of its repressor miRNAs by using a quantitative kinetic model. By integrating information and data from the literature and databases, we constructed a regulatory map of p21 including its targeting miRNAs, TFs and interacting proteins. Using this map, we found potential ‘housekeeping’ miRNAs that are promoted by a few TFs and are associated with many distinct cellular functions. Based on this map, we derived a kinetic model to test hypotheses about mechanisms of collective miRNA repression. We computationally detected several couples of miRNAs that may cooperate in the repression of p21 and experimentally validated the cooperativity between two of them, miR-572 and miR-93. We further used these two miRNAs as an example to investigate three hypothetical target regulation mechanisms by multiple miRNAs. We found that one miRNA is sufficient to completely repress a target when highly expressed. However, the same effect can be achieved by two cooperating miRNAs, both of which are only expressed on a medium level. The model was further used to predict miRNA and p21 expressions for nine cellular functions, and it predicted high levels of p21 expression during DNA damage, DNA repair, senescence and migration, whereas low expressions of p21 for cell proliferation and apoptosis.

MATERIALS AND METHODS

General overview of the methodology

The aim of our analysis is to unravel the complex mechanisms by which gene and signalling networks are regulated by multiple miRNAs. We iteratively integrate data from the biomedical literature, high-throughput experiments and biological databases into a kinetic model of miRNA target hub regulation. The kinetic model is then used to formulate and test hypotheses about mechanisms of target regulation and cellular function-related variability. In short, the approach includes three blocks (Figure 2): (i) data retrieval, (ii) construction and analysis of the regulatory map and (iii) kinetic modelling and simulation. These blocks are discussed in detail below.

Figure 2.

Schematic representation of the approach. Biomedical knowledge and quantitative data are integrated using mathematical modelling to describe the regulatory role of miRNAs in signalling pathways and gene regulatory networks. It considers data retrieval, iterative cycles of model construction and calibration and computational simulations. The mathematical model obtained is suitable to test hypotheses and simulate complex biological scenarios associated with cellular function-related variability.

Data retrieval

Information about protein interactions was extracted from the Human Protein Reference Database (release 9.0) (11) and the database STRING (release 9.0) (12). A list of experimentally verified TFs for the considered target hub was generated from the literature. This list was complemented by putative TFs that are associated with conserved TF-binding sites (human, mouse and rat) residing in the 5-kb upstream region of the target gene. This information was extracted from the table of TFs with conserved binding sites (hg18) in the USCS genome browser (13). Information about miRNA–target interactions can be extracted either from databases of validated interactions, e.g. miRecords (14), Tarbase (15) and miRTarBase (16), or from predicted ones, which can be found in databases like miRWalk (17) and miRGen 2.0 (18). In the case of p21, we used the information from Wu et al. (10) where a list of predicted p21-regulating miRNAs was subjected to experimental validation. There exist other miRNAs than those included in our network that can regulate p21 expression (10,19). However, in order to enable the construction of a kinetic model, we considered only those miRNAs for which there are consistent quantitative data describing their repression abilities on p21. A list of TFs controlling the expression of the miRNAs was constructed using information of experimentally proved TFs of miRNAs contained in TransmiR (release 1.0) (20). In addition, we generated a list of putative TFs of miRNAs with binding sites in the 10-kb upstream region of the miRNA genes from the databases PuTmiR (release 1.0) and MIR@NT@N (version 1.2.1) and from the table of TFs with conserved binding sites (tfbsConsFactor) in USCS hg18 (13,21,22). We selected this region following the method in Shalgi et al. (3), where the 10-kb upstream region is indicated as ‘putative regulatory region of miRNAs’. For further investigation, we considered only the miRNA TFs that are experimentally verified or predicted in two of three resources.

Construction and analysis of the regulatory map

We first converted the information described in the previous section into a regulatory network map in CellDesigner (23). For assessing the reliability of the network structure, we computed a confidence score for each interaction included in the p21 regulatory map, similar to the procedure used in several interaction databases like STRING, iRefWeb or IntAct. The computed scores ranged between 0 (totally uncertain interaction) and 1 (most reliable), and they were established by integrating weighted scores for publications reporting an interaction, experimental method(s) used, interaction type and computational predictions. The scoring system was inspired by the one used in IntAct (24). A detailed description of how the confidence scores were calculated is provided in Supplementary Materials Section 1. All the interactions (p21–protein, TF–p21, TF–miRNA and miRNA–p21) displayed in the map and their corresponding scores can be found in Supplementary Excel file.

In our map, the network motifs including a TF that regulates p21 and a p21-targeting miRNA were considered as feedforward loops. Those in which p21 is repressed by a TF and a targeting miRNA were classified as coherent feedforward loops. Those in which p21 is activated by a TF and repressed by an miRNA were classified as incoherent feedforward loops. Meanwhile, the network motifs including a TF that promotes a targeting miRNA and interacts with p21 were classified as feedback loops. Based on Gene Ontology (GO) terms, we derived tailored TF lists of p21 and its targeting miRNAs for different cellular functions (cell proliferation, apoptosis, immune response, inflammation response, cell cycle control, DNA damage response, cell senescence, DNA repair and cell motility and migration; Supplementary Excel file).

Kinetic modelling

Based on the regulatory map, a kinetic model was constructed using ordinary differential equations (ODEs) (25). Precisely, the kinetic model accounts for the evolution in time of the mRNA (mp21) and protein (p21) expression levels of the miRNA target hub p21, the p21-targeting miRNAs considered (miR, i = 1 … 15) and the complexes formed by targeted mRNA and miRNA, [mp21|miR]. In all, the model is constituted by 32 time-dependent variables and 64 parameters:

For the p21 mRNA (mp21), processes considered in the model are basal synthesis (k) mediated by its TF (F(TF)), basal degradation (k) and association with an miRNA (k). For each p21-targeting miRNA (miR), processes considered in the model are basal synthesis (k) mediated by its TF (F(TF)), basal degradation (k) and association with the p21 mRNA target (k). For each complex ([mp21|miR]), processes considered in the model are association of miR and p21 mRNA to form a complex (k) and degradation (k). For p21 protein (p21), processes considered in the model are mRNA-mediated synthesis of protein (k) and degradation (k). Additional algebraic equations accounting for the total measurable amounts of p21 mRNA (mp21) and each p21-targeting miRNA (miR) were also included

The model variables were normalized to the non-repressed levels of p21 mRNA and protein (mp21 and p21 levels were assumed to be 1 when no miRNA repression occurs). Values of the model parameters characterizing the reaction rates were assigned either by using fixed values taken from relevant publications (e.g. protein, mRNA and miRNA half-lives) or by applying data fitting techniques (26,27). In the second case, we processed and used quantitative data, describing the repression exerted by the different miRNAs on p21 mRNA and protein levels (10), to estimate the values of the parameters k and k. These parameters were numerically calculated in biologically meaningful intervals using an iterative estimation method, which combines global and local optimization algorithms. The parameter estimation was carried out using SBtoolbox2 (28). Computational simulations were performed using MATLAB (Mathworks, MA, USA). The construction and calibration of the kinetic model is described in detail in Supplementary Materials Section 2.

Experiment methods

For validating the predictive ability of our kinetic model, we performed additional experiments. Precisely, SK-Mel-147 melanoma cells were transfected with mature miRNA mimics of two miRNAs targeting p21 (miR-572 and miR-93), either individually at a concentration of 100 nM or in combination at 50 nM concentration. Cells were pulse-treated with 250 nM of the genotoxic substance doxorubicin and protein lysates were collected at 0, 2, 4, 6, 8, and 24 h after treatment. Expression levels of p21 were measured using immunoblotting as explained in Supplementary Materials Section 4.

RESULTS

The regulatory map of p21

According to the approach described in the previous section, biomedical knowledge about p21 regulation retrieved from publications and databases was integrated into a regulatory map (Figure 3). The map constructed in CellDesigner provides an up-to-date summary of information concerning transcriptional and post-transcriptional regulation of p21 and its regulating miRNAs. Using CellDesigner, the regulatory map is compliant with the Systems Biology Graphical Notation (29). The annotated SBML file of the regulatory map is freely available at www.sbi.uni-rostock.de/resources/software/target-hub. Furthermore, the map was complemented with an annotated Excel file containing accession numbers; associated GOs; PubMed identifiers and other details of the protein interaction partners, miRNAs and TFs. In addition, this file contains the confidence score for each interaction in the map. These confidence measures were established integrating weighted scores for publications reporting an interaction, experimental method(s) used, interaction type and computational predictions (see Supplementary Excel file for details).

Figure 3.

The regulatory map of p21. Visualization of the annotated SBML file integrating information about the regulation of miRNA target hub p21. miRNAs are classified according to their mechanisms of p21 regulation (deadenylation denoted by blue box and translation repression denoted by green box). TFs are divided into three groups: TFs regulating only p21 expression (red box), TFs regulating only miRNAs (yellow box) and those that regulate the expression of p21 and some of the miRNAs (light blue box). The latter accounts for putative feedforward loops. p21 protein interaction partners are clustered into groups (grey boxes). The small boxes, highlighted in green, account for putative feedback loops formed by p21, given miRNAs and the interacting partner (e.g. p21→AKT1→miR-132⊣p21). Highlighted in blue are TFs that down-regulate expression of the considered gene targets. The purple box accounts for cellular functions assigned to the considered TFs. A high resolution version of the regulatory map is available in Supplementary Data.

Overall, the map suggests a complex regulation of the miRNA target hub p21 with multiple transcriptional and post-transcriptional regulatory motifs. According to our analysis, certain TFs, including EGR1 and SP1, promote the expression of groups of miRNAs regulating p21. This supports the hypothesis that miRNA cooperativity plays a role in the repression of p21. GO analysis shows that some of the cellular functions, in particular cell proliferation, apoptosis and immune response, are prominent in the regulation of the p21-targeting miRNAs by TFs (Supplementary Materials Figure S4A). For example, miRNA-345 and miR-93 are predicted to be promoted by three and seven TFs, respectively. However, the coverage of cellular functions of miR-345 (six functions) is more extensive than miR-93 (three functions). This finding implies that miRNAs like miR-345, which are promoted by a few TFs but associated with many cellular functions, can be considered as ‘housekeeping’ miRNA. On the other hand, miRNAs like miR-93, which are promoted by many TFs but associated with a few cellular functions, can be considered cellular function-specific miRNAs (Supplementary Materials Figure S4B).

Kinetic modelling of target hub regulation by multiple miRNAs

We derived a kinetic model with a structure suitable to investigate the hypothesis of enhanced miRNA repressive ability associated with miRNA-binding sites in close proximity on the 3′-untranslated region (UTR) of the target mRNA (6,7). Enhanced repression in this context originates from a cooperative interaction between the two miRNAs whose binding sites are close. To test the consequences of this hypothesis, we first predicted the binding sites for p21-targeting miRNAs on its 3′-UTR and then we derived a kinetic model to analyse mechanisms for pairs of miRNAs repressing the same target hub. Figure 4 shows the matrix of potentially cooperating and non-cooperating miRNA pairs targeting p21 for post-transcriptional regulation (see Supplementary Materials Section 2.1 for explanation and Supplementary Excel file for binding patterns). Our results indicate that several couples of miRNAs targeting p21 display the conditions as described by the previous studies (6,7) and thus can cause enhanced target repression. The only exception is miR-639, whose binding site neither overlaps with nor is in close proximity to the other miRNA-binding sites.

Figure 4.

(A) Matrix of cooperating and non-cooperating p21 targeting miRNAs. The intersections of pairs of miRNAs denote their potential for cooperation. Grey cells indicate a non-interacting pair and blue cells denote potentially cooperating pairs based on binding site proximity range (13–35 nt) defined by Saetrom et al. (7). Red cells, however, indicate that cooperation cannot be established due to binding sites with extensive overlap or miRNA pairs that share the same binding site. The figures inside specify the fraction of binding sites interacting with the respective partner (the figures in brackets represent the miRNA on the x-axis). (B) Binding sites of regulatory miRNAs in the p21 3′-UTR.

Using the approach previously described and based on the regulatory map, we derived the kinetic model using ODEs, accounting for the evolution in time of the mRNA and protein expression levels of the target hub p21, the p21-targeting miRNAs and the complexes integrated by the mRNA and the miRNAs. To substantiate the cooperative effect associated with the proximity of miRNA-binding sites, we defined a group of new variables ([mp21|miR]), accounting for ternary complexes composed of p21 mRNA and two potential cooperating miRNAs (miR and miR). For these new variables, we considered two processes: (i) the association of p21 mRNA with miR and miR into the complex (k) and (ii) the degradation of the complex (k). Their dynamics are represented by the following differential equations, which complement the already established set of equations for the particular case of miRNA cooperativity:

Model validation

For validating the structure of the kinetic model, we carried out overexpression experiments for miR-572 and miR-93, two of the p21-targeting miRNAs that may cooperate to regulate p21 according to our analysis. Cells were transfected with those miRNAs either individually (100 nM) or in combination (50 nM each) and treated with doxorubicin, a genotoxic stress-inducing agent. Expression levels of p21 were measured by immunoblotting at different times after doxorubicin treatment. Thus, we obtained the data characterizing the dynamics of p21 expression after genotoxic stress in four scenarios: (i) endogenous miRNA expression, (ii) only miR-572 overexpressed (100 nM), (iii) only miR-93 overexpressed (100 nM) and (iv) both miRNAs partially overexpressed (50 nM each). We further simulated the model by configuring it to the designed experimental scenarios and compared the model predictions with the experimental data (see Supplementary Materials Section 3.1 for details). As shown in Figure 5, the model predictions are in agreement with the experimental observations for p21 response after genotoxic stress. Overall, overexpression of miR-572 or miR-93 is able to reduce the p21 up-regulation after genotoxic stress. In the scenarios 2 and 3 the total amount of the transfected miRNA is equal. However, the differences in p21 response are due to the different repression efficiencies of the two miRNAs. Furthermore, the data in scenario 4 indicate that the combined partial overexpression of both miRNAs, with binding sites in close proximity and therefore potentially cooperating, induces stronger down-regulation of p21 response. These results validate the hypothesis of cooperativity for this couple of miRNAs and also the ability of our method to detect and characterize miRNA cooperativity.

Figure 5.

Regulation of p21 expression by miR-572 and miR-93. The experimental data of p21 protein expression in response to genotoxic stress (Exp) are compared with the model predictions (Sim) in four biological scenarios: (i) both miRNAs are normally expressed (NTC), (ii) miR-572 is overexpressed (miR-572), (iii) miR-93 is overexpressed (miR-93) and (iv) both miRNAs are overexpressed (miR-572 + 93). n.u.: normalized unit.

To further test the predictive ability of the model, we used it to predict miRNA-regulated p21 expression patterns in different biological contexts. First, we extracted and normalized the tissue-specific miRNA expression levels from the database miRNAmap (release 2.0) (30). Next, for each miRNA expression profile obtained, we computed the steady-state levels of p21 by changing accordingly the initial concentrations of the corresponding miRNAs in our model. Then, we retrieved qualitative information about p21 expression levels in the same tissues from the database ArrayExpress (version as of January 2012) (31). For making the simulation results and experimental data comparable, the tissue-specific p21 expression levels were discretized in two categories (low and high). The process is explained in detail in Supplementary Materials Section 3.3. Assuming that epigenetic regulation exerted by miRNAs plays an important role for tissue-specific protein expressions, our model is able to correctly predict the relative p21 expression levels in 9 of 12 tissues (75%; Supplementary Materials Figure S6).

Computational analysis of miRNA cooperativity

Using the calibrated and validated model, we hypothesized three mechanisms for the target regulation conducted by pairs of miRNAs: (i) independent target regulation, in this case the target site proximity does not produce an effective enhanced repression (represented by setting k = 0); (ii) interdependent target regulation, in which the combined binding of both miRNAs is required to repress the target (k = 0, k = 0 and k > 0) and (iii) synergistic target regulation, which is a combination of the previous two target regulation mechanisms (k > 0, k > 0, k > 0). We investigated the consequences of the three mechanisms for target repression by using the miRNA pair miR-93 (miR) and miR-572 (miR) as an example (see Supplementary Materials Section 3.1 for details). Figure 6 shows the course of the p21 protein expression levels for the three repression mechanisms upon different initial concentrations of both miRNAs, miR(0),miR(0) [10−1, 102]. To make our result comparable with the experimental results generated by Wu et al. (10), we computed the expression level of p21 protein at 48 h.

Figure 6.

The three regulatory mechanisms proposed, particularized for the miRNA pair miR-93 and miR-572. For each depicted target regulation mechanism, the p21 expression level is computed for different combinations of the initial concentrations of miR-93 and miR-572. The black symbols stand for four different miRNA initial concentration scenarios: both miRNAs are normally expressed ([1, 1], filled inverted triangle), miR-572 is overexpressed ([1, 10], filled diamond), miR-93 is overexpressed ([10, 1], filled circle) and both miRNAs are overexpressed ([10, 10], filled triangle). The legend bar represents the p21 expression ranging from basal levels (1) to full repression (0). n.u.: normalized unit.

In case of independent target regulation (Figure 6A), sufficient up-regulation of one of the two miRNAs is able to induce the complete repression of p21 protein synthesis. In combination with the up-regulation of a second miRNA, a reduced overall level of the miRNAs is required to silence p21. In contrast, effective p21 silencing is possible only when both miRNAs are sufficiently up-regulated. In other words, the independent overexpression of a single miRNA is not able to induce p21 silencing, rather results in a low miRNA-specific gene down-regulation (Figure 6B). This mechanism may explain the relatively poor repression ability of some miRNAs when they are experimentally overexpressed alone (4). Finally, p21 silencing can be achieved through synergistic target regulation in two ways (Figure 6C): (i) independent overexpression of any of the two miRNAs, which requires high levels of miRNA expression (at least 102-fold upregulation) and (ii) combined modest upregulation of the two cooperative miRNAs, which reduces the amount of miRNA required for gene silencing by at least one order of magnitude. Moreover, through this mechanism, the system becomes more sensitive to miRNA expression changes both for individual and combined miRNA down-/up-regulation.

Based on the previous results (Figure 5), we can say that synergistic target regulation is the most likely mechanism by which p21 is repressed by the two cooperating miRNAs, namely, miR-572 and miR-93. Therefore, we further investigated this mechanism for different strengths of the miRNA cooperativity (named as K) and simulated four scenarios with different initial miRNA concentrations ([miR(0), miR(0)]): (i) both miRNAs are normally expressed ([1,1], blue triangles), (ii) miR-93 is normally expressed, miR-572 is overexpressed ([1,10], red diamonds), (iii) miR-93 is overexpressed, miR-572 is normally expressed ([10,1], green circles) and (iv) both miRNAs are overexpressed ([10,10], purple triangles). For each scenario, we defined an interval accounting for the strength of the miRNA cooperativity (K = k × [10−5 105]). We simulated the model, assigning values of K in steps of 100.25; for every round of simulation, we computed the expression levels of p21 protein at 48 h (Figure 7A). According to the value of K, we divided the emerging simulation results into three regions: Region I, in which the cooperative target repression induced by the pair of miRNAs is negligible compared with independent target repression; Region III, in which the repression of p21 is mainly driven by the miRNA cooperative interaction and Region II, in which independent and cooperative target regulation contributes equally.

Figure 7.

(A) p21 expression levels at different strengths of cooperativity (K) between miR-93 and miR-572. The dashed vertical lines represent thresholds of K, for which the miRNA-mediated p21 repression displays different behaviours. Different symbols correspond to different miRNA initial concentration scenarios: both miRNAs are normally expressed ([1, 1], filled inverted triangle), miR-572 is overexpressed ([1, 10], filled diamonds), miR-93 is overexpressed ([10, 1], filled circles) and both miRNAs are overexpressed ([10, 10], filled triangles). The small plot zooms in on the first scenario for illustrating the sigmoid shape of the expression levels of p21 at different K. (B) The role of RBPs in multiple miRNA-mediated target repression. Our analysis suggests the activity of RBPs in combination with the repression by multiple miRNAs here described can induce a tunable-like target repression. n.u.: normalized unit.

When the strength of cooperativity (K) is weak (Figure 7A, Region I), both miR-572 and miR-93 work independently to repress p21 gene expression. The strongest p21 repression is achieved when both are overexpressed (purple triangles). Overexpression of one or the other miRNA leads to different levels of target repression due to their different individual repression efficacies. For very high values of K (Figure 7A, Region III), p21 repression saturates and enters into an interdependent mode in which the cooperative interaction overcomes the repression by individual miRNAs. p21 silencing is achieved when both miRNAs are overexpressed (purple triangles). In the intermediate region that describes target repression for moderate values of K (Figure 7A, Region II), the performance of the target repression increases gradually. The most remarkable and steep decrease in p21 expression levels can be observed when both miRNAs are overexpressed (purple triangles). The results in Region II suggest that active modulation of miRNA cooperativity can play an important role in determining the efficacy of target repression. In line with this, some RNA-binding proteins (RBPs) like fragile X mental retardation protein (FMRP) or PUmilio-Fem-3-binding factor proteins (PUF) can enhance target repression by interacting with the miRNA-induced silencing complexes (miRISCs), while others (e.g. DND1 and HuR) are reported to counteract miRNA-mediated repression (32). We hypothesized that the activity of those RBPs in combination with the repression by multiple miRNAs here described induces a tunable-like target repression. According to this, the mechanism of synergistic target regulation can induce a sophisticated regulation of miRNA targets, in which regulation of RBPs activity shifts the performance of target repression from one region to another (Figure 7B).

The response of p21 to upstream stimulus signals

We further investigated whether the combination of changes in miRNA abundance and the cooperativity between miRNAs can affect the features of p21 response (p21) to transient stimulation (e.g. upstream signal from p53 after response to DNA damage). Subsequently, we defined a transient stimulus signal activating the synthesis of p21 mRNA (mp21) as a function of µ and τ, which represent the amplitude and duration of the signal, respectively (Figure 8A, left). The upstream signal function (S(µ,τ)) was included into the model.

Figure 8.

The dynamics of p21 to upstream signals. (A) The sketch of the stimulus signal computed (left) and the p21 response (right). The parameters μ and τ account for the amplitude and duration of the stimulus signal, respectively. (B) The predicted p21 response peak to different stimulus signals for different miRNA abundance scenarios. Each plot refers to a set of stimulus signals with the same amplitude range but different durations (µ ∈ [10−2 101] n.u., τ = [1, 10, 24] h). Five miRNA abundance scenarios were considered and they are as follows: (i) (solid red line): no expression of any miRNA, (ii) (black dash-dot line): normal expression of all miRNAs, (iii) (blue circles): normal expression with miRNA cooperativity, (iv) (green dotted line): overexpression of all miRNAs and (v) (magenta dashed line): overexpression with miRNA cooperativity. The grey dashed lines represent the threshold, which is equal to 10% of the basal p21 expression level. (C) The time-series plots of the p21 response when the p21 is stimulated by the low (µ = 0.1 n.u.) or high (µ = 1 n.u.) amplitude and long-lasting (τ = 10 hr) stimulus signals. n.u.: normalized unit.

In a series of simulations, we modulated µ and τ and computed the peak of the p21 response (Figure 8A right). Assuming non-basal synthesis for each miRNA (k( … ) = 0), we simulated the p21 response for five different miRNA abundance scenarios by modulating their initial concentrations: (i) : no expression of any miRNA (miR … (0) = 0, k = k = 0), (ii) : normal expression of all miRNAs (miR … (0) = 1, k = k = 0), (iii) : normal expression with miRNA cooperativity (miR … (0) = 1, k = k + k, k = k + k), (iv) : overexpression of all miRNAs (miR … (0) = 10, k = k = 0) and (v) : overexpression with miRNA cooperativity (miR … (0) = 10, k = k + k, k = k + k). Furthermore, we computed the peaks of the p21 response to transient stimulation for different settings of stimulus signals (Figure 8B; see Supplementary Materials Section 3.2 for details). Moreover, we compared the dynamics of the p21 response for long-lasting stimulus signals with low or high amplitude (Figure 8C; see Supplementary Materials Section 3.2 for details).

For short stimulation (Figure 8B, τ = 1 h), the peak of p21 response is positively correlated with the amplitude of the stimulus signal, i.e. the stronger the stimulus signal the higher the peak. For the complete range of signal amplitude simulated, the p21 response peak is negatively correlated with miRNA abundance. That means the p21 response peak is highest if miRNAs are not expressed, whereas p21 peak is gradually reduced with the increase of miRNA abundance ( >  > ). In addition, miRNA cooperativity further lowers the p21 response peak ( > ,  > ). However, with an increasing duration of the stimulus signal (Figure 8B, τ = 10 and 24 h), the p21 response peak pattern is distorted. For signals with low amplitude (µ < 1), the p21 response peak maintains the same pattern ( >  >  >  > ). In contrast, for higher amplitude (µ > 1) the p21 response peak in scenarios with miRNA expression gradually converges towards a situation of non-miRNA repression. This behaviour is due to the consumption of the entire amount of available miRNAs (data not shown). Interesting enough, higher peaks of p21 response are reached for than for when the amplitude and duration of the stimulus signal are intense and long (µ > 8, τ = 10 h or µ > 2, τ = 24 h). This result seems contradictory to our expectation of stronger p21 repression caused by miRNA cooperation. However, the simulations of miRNA consumption dynamics indicate that under identical miRNA abundance, the model that considers cooperation shows quicker exhaustion of free miRNA, which in turn leads to higher peaks of p21 response (Supplementary Materials Figure S5A).

In addition, we defined a hypothetical threshold (10% of the basal p21 expression level) above which the p21 response peak is considered as significant. As shown in Figure 8B, the position in which the peak crosses the threshold accounts for the minimum signal amplitude required to obtain significant p21 response. The increasing miRNA abundances shift the cross point towards higher signal amplitude required (Figure 8B left). The inclusion of miRNA cooperativity further enhances this effect (Figure 8B middle and right). On the other hand, the increasing duration of the signal reduces the amplitude required. These results suggest that higher abundance of miRNAs as well as the phenomenon of cooperativity can enhance noise buffering by increasing the minimum TF activity level required to trigger significant target expression (33).

Moreover, for long stimulus duration increasing abundance of miRNA and cooperativity not only lowers the p21 response peak but also changes the shape of the p21 temporal pattern (Figure 8C). For low stimulus amplitude, the p21 response exhibits a rather flat activation, with values of p21 in the same interval for many hours, while for high signal amplitude the activation of p21 is delayed. More simulations with variation in signal duration and amplitude are included in Supplementary Materials Figure S5B.

Cooperative miRNA regulation of p21 expression for different cellular functions

We used the model to predict the p21 steady-state levels for different cellular functions and to analyse the role of miRNA cooperativity. Subsequently, we first extracted the list of GO terms associated with the TFs of the miRNAs. We then assumed activation of miRNA expression when at least one of the miRNA TFs is associated with the cellular functions under investigation. Figure 9 (left) displays the predicted activity profile of p21 targeting miRNAs for different cellular functions. In addition, we computed the p21 steady-state levels affected by cooperative vs non-cooperative miRNA action (Figure 9 right). The model predicts high levels of p21 expression during DNA damage, DNA repair, senescence and migration, while low levels are associated with cell proliferation, apoptosis and cell cycle. Moreover, in some of the cellular functions p21 steady-state levels are identical for both cooperative and non-cooperative miRNA-mediated regulation. This can be explained by the fact that the miRNAs expressed under these scenarios do not have their binding sites in close proximity and therefore do not effectively cooperate (Figure 9 right; DNA repair, cell migration and DNA damage). In contrast, cooperating miRNAs are expressed in other scenarios, and thus differences appear in the p21 steady-state levels under the cooperative and non-cooperative miRNA action (Figure 9 right; senescence and immune response). This suggests that cooperating miRNAs can be selectively expressed in specific biological scenarios.

Figure 9.

Cooperative miRNA regulation of p21 expression in different cellular functions. The statuses (on/off) of miRNAs for specific cellular functions are derived from the GO terms of their TFs. The horizontal bars represent the predicted p21 steady-state levels for different cellular functions when the miRNA cooperativity is considered (C) or not (no C). The equations and parameter values can be found in Supplementary Materials Section 3.3. n.u.: normalized unit.

CONCLUSIONS AND DISCUSSION

Structure of the mathematical model

So far only a few kinetic models of miRNA target regulation are proposed. The mathematical model of the miRNA-mediated gene silencing developed by Levine et al. (34) was continued by Nissan and Parker (35) and Zinovyev et al. (36), who developed detailed models describing the mechanism of miRNA-induced target repression. Compared to them, in our model we assumed a simplified description for the mechanism of miRNA-mediated target repression. However, all these models considered only regulation by individual miRNAs, while the model presented in this article addressed the question of the multiple and cooperative miRNAs repressing the same target. Khanin and Vinciotti (37) derived a data-driven model using microarray datasets of human mRNAs regulated by miR-124 a and used it to quantify the miRNA-mediated effect on target mRNA degradation. We followed a similar approach using quantitative data to calibrate our model, but the model here presented also quantifies the effect of miRNA repression in both mRNA and protein levels. Lee et al. (38) constructed a network model characterizing miR-204 as a tumour suppressor that was used to detect/validate 18 gene targets related to tumour progression. In line with this idea, we constructed the model based on the p21 regulatory network and used it to predict p21 expression for nine different cellular functions.

Regulation of p21

We here used a novel approach, combining bioinformatics algorithms and data from the literature and databases with kinetic modelling, to elucidate the mechanisms involved in the regulation of the miRNA target hub p21. We constructed a map based on TF prediction algorithms, published reports and databases, depicting the transcriptional and post-transcriptional regulation of p21 as well as that of its regulating miRNAs. We derived a kinetic model to test the hypotheses concerning mechanisms of collective miRNA repression and proposed at least three possible mechanisms. We computationally detected several couples of miRNAs that may cooperate and experimentally validated the cooperativity between miR-572 and miR-93. Our results suggest that these two miRNAs can induce cooperative effect to repress the expression of p21 via the mechanism of synergistic target regulation.

Through GO analysis, cellular functions were included in the map opening the possibility to investigate tissue- and development-dependent patterns of p21 regulation. Our model predicts high levels of p21 expression during DNA damage, DNA repair, senescence and migration, while low levels are associated with cell proliferation, apoptosis and cell cycle. Our predictions are supported by recent experimental publications describing p21 expression in the considered cellular functions. Interestingly, endogenous up-regulation of a number of miRNAs targeting p21 can rescue cells from p21-mediated senescence (19). These findings are in line with the model predictions that the miRNA-mediated p21 repression is poor during senescence, but becomes strong when cells re-enter cell proliferation. In addition, the model suggests that a low level of p21 expression is associated with apoptosis. It is known that one of the possible events triggering apoptosis initiation is the mitotic catastrophe in which cells fail to properly arrest cell cycle (39). Insufficient expression of cell cycle arresters like p21 and 14-3-3σ may trigger this phenomenon (40,41). Our simulations predict that an apoptosis-associated miRNA expression profile can contribute to the down-regulation of p21 during apoptosis, which could facilitate mitotic catastrophe.

Moreover, our map and kinetic model provide a useful means to predict and analyse deregulation in p21 expression associated with cancer-related phenotypes. For example, we detected a number of feedforward loops composed of p21, p21-targeting miRNAs and their common TFs (Supplementary Materials Table S1). The existence of feedforward loops has important consequences for the deregulation of p21 in cancer conditions. Interestingly, p21 has dual roles during the cell cycle process: under non-stressed conditions, p21 is expressed at low levels and promotes cell cycle progression; under stress conditions, p21 expression is increased through p53-dependent pathways and becomes a cell cycle inhibitor (42). Our results predict that the tumour suppressor p53 could regulate p21 through an incoherent feedforward loop, in which p53 directly up-regulates p21 and indirectly down-regulates p21 through the p53-promoted transcription of miR-125-5p. In non-small lung cancer cells, miR-125a-5p is down-regulated, suggesting that for some cancer-associated phenotypes the feedforward loop formed by p53, p21 and miR-125a-5p is deactivated (43). Our analysis suggests that the suppression of such regulatory loop can favour cancer progression by delaying the initiation of p21-triggered cell cycle arrest after DNA damage (Supplementary Materials Figure S7).

Design principles of networks involving miRNA target hubs

Levine et al. (34) hypothesized that the extent to which an miRNA can affect the mRNA levels of its target genes can be altered by miRNA-specific parameters, global parameters (controlling mRNA turnover) and the features of the cellular function in which repression occurs. Our results are in line with their conclusions and indicate that miRNA cooperativity or synergy in gene repression may be one of the mechanisms by which specificity is gained in the regulation of different targets by miRNAs. But it also suggests that cooperating miRNAs may induce fine-tuning in gene expression for cellular functions in which those close binding site miRNAs exhibit distinctive regulation.

Our kinetic model supports the findings of Doench and Sharp (6), confirming that repression of certain genes may require two or more endogenous miRNAs targeting binding sites that are in close proximity. Using a distance of 13–35 nt between two seed sites as criterion for allowing enhanced repression (7), our computational analysis detected many occurrences of binding sites in close proximity for p21. When integrated and analysed by mathematical modelling, the regulatory network derived here indicates that distinctive regulation, through different TFs, associated with different critical cellular functions, may be one of the mechanisms by which the cellular function determines efficiency of miRNA-mediated repression (Figure 10).

Figure 10.

Examples of function-specific regulation of p21. The TFs (rectangles) regulate p21 directly and indirectly through miRNAs (parallelograms): inflammatory response (A) and cell cycle (B).

Although we focused on pairwise interactions between miRNAs, the prediction of miRNA-binding sites for p21 uncovered larger putative motifs with three or more non-overlapping miRNA-binding sites in close proximity (e.g. miR-657, miR-298 and miR-208a; see Figure 4) or partially overlapping ones (e.g. miR-654-3p, miR-572 and miR-93). These motifs may account for the clustered miRNA-binding sites proposed by Rigoutsos et al. (44) and suggest the existence of more complex forms of nonlinear regression. In many cases, miRNAs can only induce mild repression of their targets (4). Based on our analysis, the mechanism of interdependent regulation can explain the relatively poor ability of overexpression of single miRNAs to silence genes, which is observed in vitro or in vivo (45).

Moreover, the mechanism of synergistic target regulation implies a possible sophisticated regulation of miRNA targets. In this mechanism, the strength of miRNA cooperativity is determined by the interplay between the miRISCs and RBPs (46).

Our analysis also confirms the enrichment of network motifs in the local architecture of p21 regulation, like coherent and incoherent feedforward loops containing miRNAs (3,47,48). Our map reveals 30 regulatory loops. This includes 27 feedforward loops integrated by TFs that regulate the expression of miRNAs and p21 (Supplementary Materials Table S1). This suggests that the combination of these regulatory loops with the cooperative repression by miRNAs constitutes a novel mechanism, by which transient and long-term cell function-specific responses are fine-tuned. For example, under this assumption the cellular function modulates duration and intensity of transient peaks in target protein concentration predicted for incoherent feedforward loops upon step-like TFs activation (49). This regulation can be modulated via cooperation with other miRNAs, expressed in a cell-function specific manner (Figure 11).

Figure 11.

Target hub regulation by cooperative miRNA repression and feedforward loops. The regulatory loops detected in our analysis and cooperative miRNA repression can synergize in the regulation of the target. This synergy can generate further fine-tuning in the response. As an example, we here display a module which combines an incoherent feedforward loop with transcriptional delay in miR1 expression (TF1→TgHub, TF1→miR1, miR1⊣TgHub) and cooperative target repression ([miR1 AND miR2]⊣TgHub, left). When miR2 expression is triggered by a different transcription factor TF2, it can modulate the transient peak and the long-term target expression levels. Our simulations indicate that, beyond a threshold in TF2-mediated expression of miR2, a step-like transcriptional activation of TF1 induces a transient peak in target hub expression (middle). The intensity of that transient peak is modulated by the levels of TF2, which promotes the expression of the second cooperative miRNA, miR2 (right). n.u.: normalized unit.

In our analysis, a kinetic model is used to describe the miRNA-mediated repression of targeted genes, although we are aware of previous models attempting to elucidate the details of the miRNA mechanism of repression (35,36). Our results, rather than contradicting their statements, complement them by explaining an additional level of fine-regulation in the miRNA-mediated repression via interactions between different miRNAs with common targets. Taken together, our results describe three mechanisms, based on the cooperation between miRNAs targeting the same gene, by which the miRNA target specificity could be realized. Our analysis provides a feasible explanation for the mild repression induced by single miRNAs at basal levels found in vivo (4). The counterpart of this mild repression would be the ability of miRNAs, via cooperation, to fine-tune transient or long-term responsiveness and cellular function.

Approach

The approach here discussed extends and improves our preliminary idea of integrating algorithms for target prediction and databases to develop more comprehensive computational models of miRNA regulatory networks (50). Given the existence of up to 834 target hubs in Homo sapiens (3), the approach proposed here can be used as a template for other miRNA target hub networks. Our approach can be used to span a gene regulatory network around a target hub of interest by not only implementing miRNA-mediated post-transcriptional control mechanisms but also transcriptional regulation induced by TFs of the miRNAs and the target hub. The procedure described here allows the integration of many sources of information describing several levels of regulation in those networks, including gene activation via TFs, miRNA-mediated gene silencing and signalling events. This information is organized in detailed and annotated regulatory maps. Some resources provide the information that is predicted by the use of computational algorithms (e.g. miRNA-binding sites prediction). Thus, it may be noted that the individual interactions in the p21 regulatory map cannot be taken as granted, although we tried to implement only those interactions that are likely to occur. However, we calculated for each interaction a confidence score that lies between 0 and 1 and is intuitive for experimentalist. Analysis of the confidence scores indicates that in general terms protein–protein interactions and miRNA–target interactions have higher scores (respectively, 0.72 and 0.78 on average) because for these interactions we only considered highly reliable information, essentially direct interactions that are experimentally validated. We note that there is an immense amount of published information on protein–protein interactions, while for miRNA–target interactions the purpose of our analysis limits us to consider only those whose actual repressive ability is measured. In line with this, moderately high scores were computed for TF–p21 interactions where typically at least one publication per interaction is found. In contrast, most of the TF–miRNA relationships are computationally predicted (investigation of transcriptional regulation of miRNAs is merely emerging now). Therefore, these interactions received rather low scores (0.32 average). However, to ensure certain confidence we consider only the TF–miRNA interactions predicted by at least two independent algorithms. Nevertheless, TF regulation of miRNAs requires experimental validation in most of the cases. Likewise, some of the p21 TFs are putative and demand further investigation. Thus, the analysis of our map provides the basis for further experimental investigations, especially for the non-confirmed interactions involved in the regulatory loops.

The combination of different network motifs in the visual map allows the detection of a rich landscape of positive and negative interconnected or overlapping regulatory loops. This enables the identification of potential network motifs, in particular feedback and feedforward loops. The latter ones play a noise buffering role (51). Under these conditions, mathematical modelling becomes a valuable tool to dissect the behaviour of those regulatory loops.

We show how sections of the regulatory map can be translated into ODEs and subsequently be characterized with biological information and quantitative data. The data should come from perturbation experiments, in which expression levels of the miRNAs considered are tuned (up- or down-regulated) and the effect on the target mRNA and protein levels are quantified. In the ideal case, the model should be refined through iterative cycles of parameter estimation, identifiability analysis, model structure modification and new perturbation experiments (26,27).

The deduced regulatory network can subsequently be used to investigate the dynamics that target hub expression undergoes due to the compound regulation conducted by dozen(s) of miRNAs and TFs. The proposed system of ODEs provides a solid fundament for these investigations. Full characterization in terms of model parameterization is desired but many times not achievable due to the lack of appropriate quantitative data. In this case, several analytical methodologies can be used to explore the design space of the system (52,53), a simplified version of which are used here to investigate the mechanisms of cooperative miRNA-mediated regulation in target hubs.

AVAILABILITY

The Supplementary Excel file, the high resolution version of the p21 regulatory map and an annotated SBML file of the kinetic model are provided at www.sbi.uni-rostock.de/resources/software/target-hub.

SUPPLEMENTARY DATA

Supplementary Data are available at NAR Online: Supplementary Tables 0–5, Supplementary Figures 0–7, Supplementary Excel file, Supplementary Text file and Supplementary References [3,6,10-12,21,22,30-32,42,43,47,49,51,54-70].

FUNDING

German Federal Ministry of Education and Research (BMBF) as part of the project CALSYS-FORSYS [0315264 to X.L. and J.V.]; miRSys-eBio [0316175A to X.L., J.V., O.W.]; ROSAge [0315892A to J.V., O.W., M.K.]; DFG [WO 991/4-1 to U.S., S.K.G. and O.W., KU 1320/5-1 to M.K.]. Funding for open access charge: University of Rostock.

Conflict of interest statement. None declared.