A mathematical model of bimodal epigenetic control of miR-193a in ovarian cancer stem cells.

Accumulating data indicate that cancer stem cells contribute to tumor chemoresistance and their persistence alters clinical outcome. Our previous study has shown that ovarian cancer may be initiated by ovarian cancer initiating cells (OCIC) characterized by surface antigen CD44 and c-KIT (CD117). It has been experimentally demonstrated that a microRNA, namely miR-193a, targets c-KIT mRNA for degradation and could play a crucial role in ovarian cancer development. How miR-193a is regulated is poorly understood and the emerging picture is complex. To unravel this complexity, we propose a mathematical model to explore how estrogen-mediated up-regulation of another target of miR-193a, namely E2F6, can attenuate the function of miR-193a in two ways, one through a competition of E2F6 and c-KIT transcripts for miR-193a, and second by binding of E2F6 protein, in association with a polycomb complex, to the promoter of miR-193a to down-regulate its transcription. Our model predicts that this bimodal control increases the expression of c-KIT and that the second mode of epigenetic regulation is required to generate a switching behavior in c-KIT and E2F6 expressions. Additional analysis of the TCGA ovarian cancer dataset demonstrates that ovarian cancer patients with low expression of EZH2, a polycomb-group family protein, show positive correlation between E2F6 and c-KIT. We conjecture that a simultaneous EZH2 inhibition and anti-estrogen therapy can constitute an effective combined therapeutic strategy against ovarian cancer.

Introduction

Ovarian cancer is the most lethal gynecological malignancy and the 5th leading cause of cancer death among women [1]. As ovarian cancer has few symptoms early in its course, the majority of patients are diagnosed with late stages (III and IV) of the disease. The 5-year survival rate is generally less than 20% for patients with advanced-stage disease despite therapeutic advances, whereas survival rate for patients with stage I or II disease is greater than 80% for the same period [2]. Although current chemotherapeutics demonstrate a more than 90% complete response rate in early-stage tumors, only a 20–30% partial response rate is observed in advanced stage tumors as well as in chemo-resistant relapsed tumors [3]. A better understanding of the molecular changes of ovarian carcinogenesis may lead to better therapeutic strategies for this deadly disease.

An emerging hypothesis states that cancer arises from a small population of self-renewing cancer initiating cells (CIC) [4]. These CICs are thought to possess tumorigenic potential and enhanced drug resistance within a tumor, and are able to repopulate tumor colonies in vivo. We have previously isolated ovarian CICs (OCICs) from ovarian cancer patients [5]. These OCICs exhibit increased drug resistance towards cisplatin and taxol, and are characterized by the expression of several cell surface markers, including c-KIT (CD117, or stem cell factor receptor). However, the mechanism of how OCICs arise and how these cell surface markers are transcriptionally controlled are not fully understood.

Estrogen is an important regulator of growth and differentiation in normal ovaries [6] and is involved in the development of ovarian cancer [7]. In this regard, studies have explored the risk of ovarian cancer among women using hormone replacement therapy (HRT) in the treatment of menopausal symptoms [8], [9], [10]. For example, a study involving a million women suggested that women treated with HRT were associated with an increased risk of ovarian cancer [11]. However, the role of estrogen in ovarian carcinogenesis is not fully understood.

Epigenetic modifications in the genome play a crucial role in transcriptional control and aberrant epigenetic modification is now considered a hallmark of cancer [12]. Epigenetic modifications are responsible for controlling gene expressions that allow for specific phenotypes [13], [14]. Among them, DNA methylation is one of the most characterized epigenetic modification which takes place at the 5′ position of cytosine in CpG dinucleotides resulting in the formation of 5-methylcytosine. DNA methylation which occurs at the clusters of CpG sites (called CpG islands) in the promoter region of a tumor suppressor gene results in transcriptional repression and tumorigenesis in some cases [15]. Also considered as crucial regulators of gene expression, miRs are referred to as tumor suppressor-miRs if they regulate the expression of an oncogene [14]. Studies found that tumor suppressor-miRs are frequently down-regulated by genetic or epigenetic mechanisms leading to the up-regulation of oncogenes and acceleration of tumorigenesis [16]. For example, miR-34b/c (which resides at the promoter CpG island of a gene) was found to be epigenetically silenced by DNA methylation in metastatic cancer cells [17]. Re-expression of miR-34b/c suppressed cancer invasion in vitro and in vivo.

We are interested in miR-193a which is located at chromosome 17q11.2 and embedded in a CpG island. Studies have found that miR-193a was epigenetically silenced by promoter hypermethylation in acute myeloid leukemia (AML) and lung cancer [18], [19]. Importantly, miR-193a targets the stem cell marker, c-KIT, for repression in AML [18], [20]. Another experimentally confirmed target of miR-193a is E2F6 [21]. Up-regulation of E2F6 has been observed in breast cancer cells upon treatment with estrogen [22]. Within the E2F family of transcription factors which are involved in cell cycle control through transcriptional activation or repression [23], E2F6 was found to be a transcriptional repressor capable of associating with the histone-lysine N-methyltransferase EZH2, which is an important constituent of the polycomb complex [24], [25], and with DNA methyl-transferase DNMT3b, allowing for its recruitment to the repressive complex [26].

As the role of miR-193a in ovarian cancer is currently unknown, we hypothesize that E2F6 protein may suppress the expression of miR-193a, both directly as a DNA-binding transcriptional repressor and indirectly through promotion of the polycomb complex assembly. Indeed, the existence of the binding sites for E2F6 protein in the promoter of miR-193a as demonstrated by ENCODE ChIP-Seq data [27] seems to suggest that long-term suppression of miR-193a by E2F6 may lead to epigenetic silencing of miR-193a. These events may lead to the up-regulation of c-KIT and trigger the onset of the OCIC phenotype, thus promoting ovarian carcinogenesis. In this study, we use a mathematical model to explore the intricacies of this epigenetic control of miR-193a which could prove to be significant as it predicts a switching behavior in c-KIT that could be critical for ovarian carcinogenesis.

Materials and Methods

Cell culture

Ovarian cancer cell line A2780 was propagated in RPMI1640 (Invitrogen, Carlsbad, CA) supplemented with 10% FBS, and 50 units/ml of penicillin/streptomycin. The cell was incubated at 37°Cwith 5% CO2.

Cloning of miR-193a expressing vector

The sequences which transcribe miR-193a pre-miRNA was amplified by PCR using specific primers F: 5′-CGCGGATCCAGTTTCTCGGCGCATAACTC-3′ and R: 5′-CCCAAGCTTCGCTATTTCTCCAGCGAAGTG-3′, using genomic DNA of IOSE cells which express miR-193a. The PCR products were ligated into yT&A Cloning Vector (Yeastern Biotech, Taiwan) for sequencing confirmation. miR-193a-yT&A was digested with BamH I and Hind III, and inserted into multiple cloning site of pSilencer 4.1-CMV puro siRNA expression vector, which was predigested with BamHI and HindIII.

Extraction of RNA and Quantitative RT-PCR

RNA extraction was performed using TRIzol reagent (Invitrogen Carlsbad, CA) according to the manufacturer’s protocol. To remove potential contaminating DNA from the complementary DNA, 1 µg of total RNA was treated with DNase I (Amplification Grade, Invitrogen) prior to reverse transcription. RNA was then reverse-transcribed using random primers (for c-KIT or E2F6 expression) or TaqMan reverse transcription kit (Applied Biosystems, Foster City, CA, USA). Quantitative real-time RT-PCR was then performed using ABI StepOne real-time PCR system (Applied Biosystems). Primers are available upon request. Relative expression of miR-193a, c-KIT and E2F6 was calculated using the comparative Ct method.

Results

miR-193a targets c-KIT and E2F6

To investigate the relationship between miR-193a and c-KIT/E2F6 in ovarian cancer (Fig. 1A), we over-expressed miR-193a in A2780 ovarian cancer cells which have low levels of miR-193a expression. Expression of miR-193a resulted in a significant down-regulation of c-KIT mRNA and E2F6 mRNA, thus suggesting that c-KIT mRNA and E2F6 mRNA are targets of miR-193a in ovarian cancer cells (Fig. 1B).

Figure 1

miR-193a repress E2F6 and c-KIT.

(A) Schematic diagram of c-KIT and E2F6 mRNA showing predicted miR-193a binding site at their 3′UTR. The detailed binding sequences (MREs) of these sites are also shown along with miRNA sequence alignments. (B) qPCR analysis of miR-193a, c-KIT and E2F6 in A2780 ovarian cancer cells. Overexpression of miR-193a in A2780 cells leads to down-regulation of c-KIT and E2F6 mRNAs.

Model equations

Interactions between the components of the miR-193a–E2F6–c-KIT regulatory module are shown in Fig. 2A. Tight coupling of network components and potential non-linear kinetics can give rise to complex system dynamics, which can be understood with the help of a mathematical model. In our kinetic model (Fig. 2B) there are six components, which represent: miR-193a (R m), E2F6 mRNA (R e), c-KIT mRNA (R c), the complex of miR-193a with E2F6 mRNA (R em), complex of miR-193a with c-KIT mRNA (R mc), and E2F6 protein (P). The chemical mass-action equations for the respective concentration of these molecules or complexes, R m, R e, R c, R em, R mc, and P, are stated as follows:

Figure 2

The proposed miR-193a regulatory network.

(A) Activation of estrogen receptor (ER) signaling may increase the increase the transcription of E2F6 which may subsequently increase the expression of c-KIT through miR-193a mediated ceRNA mechanism. On the other hand, in addition to transcriptional repressor, E2F6 may lead to epigenetic silencing of the miR-193a promoter resulting in the up-regulation of c-KIT. The description for each number in the pathway is listed in S1 Table. (B) The biochemical interaction for this proposed miR-193a regulatory network. Model variables names for molecule species are also given.

Here, the association rate of miR-193a–E2F6 mRNA is k 1 and association rate of miR-193a–c-KIT mRNA is k 2; the dissociation rate of complexes miR-193a–E2F6 mRNA is k m1 and miR-193a–c-KIT mRNA is k m2; the translation rate of E2F6 protein is k 3; the transcription rate of miR-193a is k 4, of E2F6 mRNA is k 5 and of c-KIT mRNA is k 6; the degradation rate of miR-193a is δ m, of E2F6 mRNA is δ e, of c-KIT mRNA is δ c, of miR-193a–E2F6 mRNA is δ em, and of miR-193a–c-KIT mRNA is δ mc and the degradation rate of E2F6 protein is δ p; and the strength of the miR-193a transcription inhibition by E2F6 protein is denoted as K 4.

The first term on the right-hand side of Eq. (1) corresponds to arrow 9 in Fig. 2A expressing miR-193a transcription inhibition by, primarily, E2F6 protein, but also by polycomb complex proteins activity. The value of the parameter K 4 is thus intended to account for both the direct transcription repression by E2F6 and for the long-term gene silencing performed by and dependent of the levels of EZH2 and DNMT3b. (The concentration of E2F6 is a variable of the dynamical model while levels of essential polycomb proteins, which are not explicitly present in model equations, are not subjects to regulation by the module and their levels are assumed constant at least at the module equilibration time scale.) The second (respectively, third) term on the right-hand side of Eq. (1) is the rate of miR-193a release from the complex with E2F6 mRNA (respectively, the complex with c-KIT mRNA). Interactions corresponding to arrows 6 and 7 in Fig. 2A, expressing hybrid mRNA complex formation, are given, respectively, by the fourth and the fifth term on the right-hand side of Eq. (1). The last term in Eq. (1) is a first-order miR-193a degradation term.

The rate of constitutive E2F6 gene transcription (arrow 1 in Fig. 2A) is k 5 (Eq. (2)). This rate depends on the level of hormonal stimulation and, as our model describes simplified regulation in ovarian cells, we will directly associate k 5 with the level of estrogen. The second term in Eq. (2) represents the release of E2F6 mRNA from the complex with miR-193a, whereas the third term is the rate of the opposite binding reaction. The last term in Eq. (2) is a degradation of E2F6 mRNA. The terms on the right-hand side of Eq. (3) are analogous. The first term in Eq. (4) is the rate for the binding of E2F6 mRNA and miR-193a. The second term in Eq. (4) stands for the opposite reaction, and the third term represents the degradation of R em. The terms on the right-hand side of Eq. (5) are analogous. The last equation accounts for E2F6 translation (arrow 8 in Fig. 2A) and degradation.

Steady States

Steady states give the long-term behavior of the system. The steady states of system (1)–(6) can be obtained by equating the right-hand sides of the equations to zero. Let S = (R m, R e, R c, R em, R mc, P), be any steady state of system (1)–(6). The steady state S satisfies the following equations where

and

We note that once R m is determined by Eq. (12), then the other components of the steady state S can be determined by Eqs. (7)–(11). Hence the condition for the existence of the steady state S is that k 4 lies in h[0, ∞), the range of the function h. We are interested in threshold or switching behavior of system (1)–(6), particularly in switching of the expression of c-KIT mRNA. Thus, the conditions on the parameters for the existence of multiple steady states are relevant. Here we give a criterion for the existence of three steady states of the model as follows: where is the maximal value of on the interval [0, ∞), and is given by

The proof for the criterion (14)–(15) is given in the S1 Text.

The Steady State Bifurcation Diagram in the (k 5, R c) plane

In view of Eqs. (7)–(8), the steady states of R e and R c increase or decrease in the same direction. This model prediction is consistent with the ceRNA hypothesis [28], [29] that the expression of E2F6 positively correlates with that of c-KIT.

As mentioned before, k 5 is an experimentally controllable parameter associated with the addition of hormones (e.g., estrogen), and K 4 is a measure of the inhibition efficiency of, primarily, E2F6 protein against miR-193a. As we will see in the discussion below, the interplay between k 5 and K 4 can generate an over-expression of c-KIT mRNA. Hence, it would be interesting to see the dependence of the steady state on the parameter k 5 for different values of K 4 – this is shown in Fig. 3. This diagram is referred to as a “steady-state bifurcation diagrams”, and k 5 is referred to as a bifurcation parameter. To generate experimental curves similar to those shown in Fig. 3, one can design an experiment in which cells are cultured at different concentrations of estrogen, and then measure the long-term mRNA levels.

Figure 3

Steady-state bifurcation diagrams.

The dependence of the steady states of variable R c (RNA level of c-KIT) on parameter k 5 (transcription rate of E2F6) for difference values of the parameter K 4 (inhibition efficiency of E2F6). The value of parameters are listed in S2 Text.

With parameters satisfying Eq. (14)–(15), the model predicts that there is a range of k 5 for which the system admits three coexisting steady states. For instance, for K 4 = 0.0001 (i.e., a lower inhibition efficiency; Fig. 3, green line) and k 5 is between 0.0556 and 0.1585, the model gives three steady states (out of which the high R c and the low R c states are stable; see the curve associated with K 4 = 0.0001 in Fig. 3). The importance of the right knee in a steady-state bifurcation diagram lies in the following fact: as the control parameter k 5 is increased from the low value up through the critical value k 5 R associated with the right knee, the corresponding steady state level of R c will increase and, due to the instability (respectively, stability) of the steady states on the middle branch (respectively, the upper branch) in the steady-state bifurcation diagram, the system will jump to the steady state on the upper branch at the critical value k 5 R, and then stay on the upper branch (see Fig. 3). There is a sharp difference in c-KIT mRNA levels between the steady states of R c on the lower and upper branches in the steady-state bifurcation diagram, which would suggest cancer phenotype-promoting over-expression of c-KIT mRNA for k 5>k 5 R. On the other hand, if one is initially at the steady state on the upper branch which is associated with the over expression of c-KIT mRNA, one can decrease the k 5 value to reduce the steady state level of R c. Moreover, if one continues to decrease the k 5 value below the critical value k 5 Lassociated with the left knee, then the corresponding steady state level of R c will decrease and, due to the instability (respectively, stability) of the steady states on the middle branch (respectively, the lower branch) in the steady-state bifurcation diagram, the system will descend to the steady state on the lower branch at the critical value k 5 L (see Fig. 3). This would recover the normal c-KIT mRNA level.

To summarize, the model predicts that there exists a threshold value k 5 R (respectively, k 5 L) in the transcriptional rate k 5 for cells to turn on (respectively, turn off) the expression of OCIC marker c-KIT. Hence, it is important to adjust the transcriptional rate k 5 to affect the expression of OCIC marker c-KIT.

Significance of the Parameter K 4

Although decreasing the value of k 5 (associated with low estrogen level) can reduce the expression of the OCIC marker c-KIT, some clinical evidence suggests that anti-estrogen therapy is only partially effective in the treatment of ovarian cancer [30], [31]. The model suggests that this might be due to the inhibition of miR-193a expression by E2F6 protein as reflected in the value of the parameter K 4 of the first term on the right-hand side of Eq. (1). This term is motivated by the fact that the expression of miR-193a depends on the transcriptional suppressive activity of E2F6 protein (associated with the level of EZH2 and DNMTs). In particular, sufficiently strong inhibition (sufficiently large K 4) would lead to a more long-term silencing of the expression of miR-193a by DNA methylation (see Eq. (1)), and thus would entail the over-expression of c-KIT mRNA; whereas for the extreme case K 4 = 0, the system (1)–(6) can only support at most one steady state (see S1 Text), which implies the non-existence of the right and left knees in the steady state diagram and loss of the bistable switch.

Let k 5 b denote the basal transcription rate constant of E2F6 mRNA. In the following, depending on the interplay between the k 5 b and K 4, we will describe two scenarios in which the second one would correspond to the clinical observation that anti-estrogen therapy are only partially effective in the treatment of ovarian cancer [30], [31].

In cells with low inhibition efficiency of E2F6 (i.e., a lower K 4), as the parameter k 5 is increased from k 5 b to k 5 2 (red solid line in Fig. 4A), the steady state of R c will first move from L 0 to L 1 via the lower branch, and, due to the instability of the steady states in the middle branch, switch from L 1 to U 1 in the upper branch, and finally move along the upper branch to U 2. On the other hand, as the parameter k 5 is decreased from k 5 2 to k 5 b (blue solid line in Fig. 4A), the steady state of R c will first move from U 2 to U 0 via the upper branch, and, due to the instability of the steady states in the middle branch, then switch from U 0 to T 0 in the lower branch, and finally move back to L 0. Hence for low inhibition efficiency of E2F6 (i.e., low K 4) and small basal rate constant (i.e., small k 5 b), addition of estrogen may turn on over-expression of c-KIT mRNA, whereas retraction of estrogen can turn off over-expression of c-KIT mRNA and make the cells recover the normal c-KIT mRNA level.

Figure 4

Steady-state bifurcation diagrams.

The dependence of the steady states of variable R c (RNA level of c-KIT) on parameter k 5 (transcription rate of E2F6) for different values of K 4: (A) K 4 = 0.1 (i.e., low inhibition efficiency of E2F6), (B) K 4 = 1 (i.e., high inhibition efficiency of E2F6). The red and blue arrow line indicates abrupt shifts of the steady state of R c with at different k values which indicated below. The value of parameters are listed in S2 Text.

In cells with high inhibition efficiency of E2F6 (i.e., a higher K 4), as the parameter k 5 is increased from k 5 b to k 5 2 (red solid line in Fig. 4B), the steady state of R c will first move from L 0 to L 1 via the lower branch, and, due to the instability of the steady states in the middle branch, switch from L 1 to U 1 in the upper branch, and finally move along the upper branch to U 2. On the other hand, as the parameter k 5 is decreased from k 5 2 to k 5 b (blue solid line in Fig. 4B), the steady state of R c will first move from U 2 to U 0 via the upper branch, and, due to the fact that k 5 b is the basal transcription rate constant of E2F6 mRNA and k 5 b is larger than the k 5 component of left-knee point, then has to stop at U 0 which is the state associated with over-expression of c-KIT mRNA. Hence for high inhibition efficiency of E2F6 (i.e., a higher K 4) and large basal rate constant (i.e., large k 5 b), retraction of estrogen cannot attenuate the expression of c-KIT mRNA. We make one important observation: for sufficiently large K 4, the k 5 value associated with the left knee is very close to zero. Hence for sufficiently large K 4, it would fit the scenario discussed in this paragraph. This suggests that for sufficiently large K 4, retraction of estrogen cannot switch off over-expression of c-KIT mRNA, which is consistent with the aforementioned clinical observation.

Validation of the model in TCGA ovarian cancer dataset

In light of the above model, we reasoned that in ovarian cancer patients with low EZH2 (i.e., low K 4), increase of E2F6 mRNA (which results from an increase of transcription rate of E2F6, k 5) may induce the expression of c-KIT mRNA (R c), through miR-193a-mediated ceRNA mechanism. However, such correlation may not exist in patients with high EZH2 (i.e., high inhibition constant of E2F6 to miR-193a, K 4) as miR-193 is epigenetically silenced. To provide an independent support for this hypothesis, we analyzed the expression microarray data in TCGA ovarian cancer dataset [32]. Interestingly, ovarian cancer patients with low EZH2 level show a positive correlation between the expression of E2F6 and c-KIT (Fig. 5A). Analysis of another set of microarray data in ovarian cancer (GDS2785) also demonstrated a similar phenomenon (S1 Fig.). However, such relationship is not observed in patients with high EZH2 (Fig. 5B and S1 Fig.) or in the correlation between E2F6 and a non-miR-193a target (S2 Fig.). Taken together, this phenomenon is consistent with our prediction (Fig. 3) that mRNA of E2F6 and c-KIT demonstrates a ceRNA relationship at low K 4.

Figure 5

Correlation between expression level of E2F6 and c-KIT in expression microarray dataset from TCGA ovarian cancer patients.

Expression microarray data from 568 ovarian cancer patients are sorted according to the expression level of EZH2. The scatter plot demonstrates the correlation between expression level of E2F6 and c-KIT in ovarian cancer patients with (A) low EZH2 (bottom 12.5%) and (B) high EZH2 (top 12.5%). The R and P-value of the Pearson correlation are shown. Interestingly, as predicted by our model, patients with low EZH2 shows a positive correlation of E2F6 and c-KIT (A, R = 0.322, P<0.005). Such correlation is not observed in patients with high EZH2 (B, R = 0.2148. P>0.05). However, patients with high EZH2 can be further separated into 2 sub-groups based on the median value of c-KIT (high, red color; low, blue color). Neither of these sub-groups demonstrates a positive correlation between E2F6 and c-KIT.

Discussion

Our mathematical model postulates the existence of a bimodal epigenetic control of miR-193a in ovarian carcinogenesis. In cells with low inhibition efficiency of E2F6, there is a ∼linear relationship between transcription of E2F6 and c-KIT mRNA, thus suggesting a ceRNA relationship between these 2 genes. In this ceRNA hypothesis [28], [29], the 3′UTRs of several mRNAs contain the miR-binding element(s) (MRE) and therefore compete for binding with the miR. Theoretically, a miR can be epigenetically controlled by both promoter hypermethylation and its ceRNA in a bimodal fashion.

It is also noteworthy to point out that this ceRNA relationship is a reversible event, such that retraction of estrogen which results in lower expression of E2F6, may reduce the expression of c-KIT (Fig. 6, upper panel). Such phenomenon can be observed in normal ovarian epithelial cells or ovarian cancer cells in which the E2F6 inhibition efficiency is low. Indeed, components of transcriptional repressors such as EZH2 which is required for transcriptional suppression of E2F6 is found to be low in normal ovarian cells and a sub-set of ovarian cancer [33], [34]. Importantly, this is also consistent with the clinical evidence that the expression of E2F6 and c-KIT are correlated in ovarian cancer patients with low expression of EZH2 but not in cases of high EZH2 expression.

Figure 6

Proposed model for ovarian carcinogenesis.

Under the condition of low EZH2 and DNMT3b expression (upper panel, as indicated by light blue and red color), E2F6 only leads to a minimal suppression of miR-193a (i.e., low K 4) and the relationship between E2F6 and c-KIT exhibits a ceRNA relationship. As a result, the amount of c-KIT in the ovarian cancer cells depends on the amount of estrogen (this information is mediated by the level of E2F6). On the contrary, under the condition of high expression level of EZH2 and DNMT3b (lower panel, as indicated by dark blue and red color), E2F6 exerts a high inhibition efficiency (i.e., high K 4). Slight increment of E2F6 results in the epigenetic silencing of miR-193a and subsequently an over-expression of c-KIT. Without miR-193a, the cells may be “locked” in a high c-KIT state resulting in ovarian carcinogenesis.

However, in cells with high inhibition efficiency of E2F6, our model predicts a highly non-linear relationship between expression of E2F6 and c-KIT. Under the condition of low expression level of E2F6, there follows a ceRNA relationship between E2F6 and c-KIT. However, at high expression of E2F6 together with high expression of transcriptional repressor such as EZH2, E2F6 may lead to epigenetic silencing of the expression of miR-193a through trimethylated H3K27 and DNA methylation (Fig. 6, lower panel). Taken together, these mechanisms may lead to a non-reversible expression of c-KIT and activation of the cell growth signaling in ovarian cancer [35]. This prediction is also consistent with the previous findings that high expression of EZH2 is associated with poor survival in ovarian cancer patients [34]. Inhibition of the EZH2 activity together with anti-estrogen therapy may therefore be effective against ovarian cancer.

In summary, due to the nonlinearities and non-intuitive dynamics of the bimodal regulation of miR-193a, we developed a mathematical model that predicts c-KIT switching behavior. Importantly, a positive correlation of E2F6 and c-KIT is only observed in ovarian cancer patients with low EZH2 expression. Combination treatment of EZH2 inhibitor together with anti-estrogen therapy may be a novel strategy against ovarian cancer.

Correlation between expression level of E2F6 and c-KIT in expression microarray dataset from ovarian cancer patients (GDS2785). Scatter plot shows the correlation between expression level of E2F6 mRNA and c-KIT in ovarian cancer patients with (A) low EZH2 (n = 16) or (B) high EZH2 (n = 17) according to the median level of EZH2 expression. The R and P-value of the Pearson correlation are also shown. Although not statistically significant, a more positive correlation between expression of c-KIT and EZH2 were observed in patients with low EZH2 (R = 0.3516, P = 0.0766) than that of patients with high EZH2 (R = 0.2938, P = 0.1729).

(TIF)

Click here for additional data file.

Correlation between expression level of E2F6 and CDK5 in expression microarray dataset from TCGA ovarian cancer patients. To examine if the observation in S1 Fig. is a random event, we analyzed the correlation between expression level of E2F6 and CDK5, a non-miR-193a target in ovarian cancer patients with (A) low EZH2 (bottom 12.5%) and (B) high EZH2 (top 12.5%) as demonstrated in S1 Fig. Positive correlation is not observed in these 2 group of patients (p>0.05).

(TIF)

Click here for additional data file.

Description of biological functions of each interaction in .

(PDF)

Click here for additional data file.

The existence of steady states.

(PDF)

Click here for additional data file.

The values of parameters in and .

(PDF)

Click here for additional data file.