Modeling the dynamics of Let-7-coupled gene regulatory networks linking cell proliferation to malignant transformation

The microRNA Let-7 controls the expression of proteins that belong to two distinct gene regulatory networks, namely a cyclin-dependent kinases (Cdks) network driving the cell cycle and a cell transformation network which can undergo an epigenetic switch between a non-transformed and a malignant transformed cell state. Using mathematical modeling and transcriptomic data analysis, we here investigate how Let-7 controls the cdk-dependent cell cycle network, and how it couples the latter with the transformation network. We also determine whether the two networks can be combined into a larger entity that impacts on cancer progression. Our analysis shows that the switch from a quiescent to a cycling state depends on the relative levels of Let-7 and several cell cycle activators. Numerical simulations further indicate that the Let-7-coupled cell cycle and transformation networks control each other, and our model identifies key players for this mutual control. Transcriptomic data analysis from the The Cancer Genome Atlas (TCGA) suggest that the two networks are activated in cancer, in particular in gastrointestinal cancers, and that the activation levels vary significantly among patients affected with a same cancer type. Our mathematical model, when applied to a heterogeneous cell population, suggests that heterogeneity among tumors results from stochastic switches between a non-transformed cell state with low proliferative capability and a transformed cell state with high proliferative property. The model further predicts that Let-7 may reduce tumor heterogeneity by decreasing the occurrence of stochastic switches towards a transformed, proliferative cell state. In conclusion, we identified the key components responsible for the qualitative dynamics of two GRNs interconnected by Let-7. The two GRNs are heterogeneously involved in several cancers, thereby stressing the need to consider patient’s specific GRN characteristics to optimize therapeutic strategies.


Introduction
Let-7 microRNAs control two distinct gene regulatory networks (GRNs) that regulate cell cycling and malignant transformation of breast cancer cells (Johnson et al., 2007;. A cyclin-dependent kinases (Cdk) network controls the correct progression of the cell cycle along the G1, S, G2 and M phases (Morgan, 2007). Growth factors (GF) and E2F stimulate, while Let-7 down-regulates the expression of several components of this cdk-dependent cell cycle network (Bueno and Malumbres, 2011). Mathematical models focusing on post-translational regulations of cyclin/Cdk complexes were proposed to account for the dynamics of the Cdk network in mammals (Novak and Tyson, 2004;Gerard and Goldbeter, 2009).
However, to our knowledge, no model has been proposed to study the impact of miRNAs on this network. Let-7 is also a key component of a GRN which promotes cell transformation in response to an inflammatory stimulus ). This GRN is characterized by a positive feedback loop (PFL), where a transient inflammatory stimulus is sufficient to induce the cells to undergo a PFL-dependent epigenetic switch from a non-transformed state toward a permanently malignant transformed state. We previously proposed a model describing the dynamics of this transformation GRN . Our model suggested that a transient inflammatory signal activates an irreversible bistable switch in the expression of the GRN components, eventually leading to a stable epigenetic switch allowing cells to display increased motility and invasiveness. In this GRN, Let-7 prevents cell transformation by inhibiting the translation of Interleukin-6 (IL6) and Ras, two drivers of the inflammatory feedback loop. Let-7 being a component common to the cell cycle and transformation networks we now raise the following questions: how does Let-7 control the cdkdependent cell cycle network; does Let-7 play a coupling role between the cell cycle GRN and the transformation GRN, and can the two GRNs be combined into a larger network that impacts on cancer progression? We address these issues using experiment-based mathematical modeling of the GRNs and by analyzing transcriptomic data from The Cancer Genome Atlas (TCGA). Our approach identifies the dynamics of the two GRNs and indicates that their activation is cancer-specific and heterogeneous from patient to patient, stressing the need to consider patient's specific characteristics.
The model of the cell cycle network is composed of a set of kinetic equations describing the temporal evolution of the levels of each component of the network. It includes the mRNAs of cyclin D, E, A, and B; the active form of E2F; the various cyclin/Cdk complexes: cyclin D/Cdk4-6, cyclin E/Cdk2, cyclin A/Cdk2 and cyclin B/Cdk1; and the active form of the Anaphase-Promoting Complex, APC, which triggers degradation of cyclins A and B at the end of mitosis ( Supplementary Fig. 1).
Variables are defined in Supplementary Table 1, the kinetic equations are in   Supplementary Table 2, and the parameters are defined in Supplementary Table 3.
As a consequence of its regulatory structure the network self-organizes with sustained oscillations in the activity of the various cyclin/Cdk complexes, which correspond to successive rounds of cell cycling (Fig. 1B-E). In the absence of both Let-7 and GF, cells proliferate and sustained oscillations of the various cyclin/Cdk complexes develop (Fig. 1B). This situation bears similarity with transformed or cancer cells, which are often characterized by downregulation of Let-7 and signalindependent growth (Sotiropoulou et al., 2009;Hanahan and Weinberg, 2011).
Starting from that condition, an increase in GF maintains cell proliferation (Fig. 1C), while an increase in Let-7 suppresses cell proliferation (Johnson et al., 2007). The latter case is characterized by a stable steady state, with low levels of each cyclin/Cdk (Fig. 1D). Finally, starting from that steady state, an increase in GF permits to recover cell proliferation (Fig. 1E).
The respective impact of Let-7 and GFs can be visualized in a two-parameter plane where the dynamical behavior of the cell cycle network is represented as a function of the synthesis rate of Let-7, V SLET7 , and the level of growth factors, GF ( Fig. 2A). For large values of V SLET7 (high levels of Let-7), the Cdk network tends to a stable steady state corresponding to cell cycle arrest regardless of GF levels. This corroborates experimental results showing that Let-7 represses cell proliferation (Johnson et al., 2007;Zhu et al., 2015). In contrast, for small values of V SLET7 (low levels of Let-7), the cell cycle network is characterized by sustained oscillations. The temporal evolution of cyclin E/Cdk2 and cyclin B/Cdk1 corresponding to conditions A to F in Fig. 2A are represented in Supplementary Fig. 2A-F, respectively. The sustained oscillations in cyclin E/Cdk2 with oscillations in cyclin B/Cdk1 of very small amplitude might correspond to endoreplication (Edgar et al., 2014) where multiple rounds of DNA replication occur without entry into mitosis (see temporal evolution in Supplementary Fig. 2D, which corresponds to condition D in Fig. 2A). Previous theoretical studies already showed that the regulatory structure of the cell cycle network in mammals is capable of generating endocycles (Gerard and Goldbeter, 2009).
The dynamics of the cell cycle network are further illustrated for different levels of Let-7 in a two-parameter plane defined by the synthesis rate of the cyclin/Cdk complexes, V SMCYC , and the level of GF, (Fig. 2B-D). By increasing Let-7, the domain of sustained oscillations, corresponding to cell proliferation, is reduced and limited to higher levels of cyclin/Cdk complexes (compare Fig. 2B, C, D where V SLET7 is equal to 0, 0.25 and 3, respectively). From a stable steady state, corresponding to quiescence (condition 1 in Fig. 2C), an increase in GF or an increase in the cyclin/Cdk levels may trigger the switch to sustained oscillations (see temporal evolution of cyclin B/Cdk1 in Fig. 2E and 2F).
We concluded that progression or arrest of the cell cycle is controlled by the relative levels of Let-7 and GF, or of Let-7 and the cyclin/Cdk complexes. Thus, when designing efficient anti-cancer strategies, the model stresses the importance to consider the relative, rather than the absolute, expression levels of network components displaying opposing effects on cell cycling. Let-7 couples the cell cycling and transformation networks. Let-7 belongs both to the cdk-dependent cell cycle GRN and the malignant transformation network (Fig.   1A). Therefore, we here determine the qualitative role of Let-7 as a coupling factor between the two GRNs, by analyzing the mutual impact of the GRNs on their respective dynamics.
Starting from a non-transformed, quiescent, cell state, defined by high Let-7 and low cyclin B/Cdk1 levels, a transient Src signal induced by inflammation triggers a stable down-regulation of Let-7, eliciting the switch of the cell cycle network from a stable steady state to sustained oscillations (Fig. 3A). Thus, transient inflammatory signals can promote persistent cell proliferation. On the opposite, as observed in the experiments  We concluded that transient modifications in the expression of the components of the inflammation-dependent bistable transformation network can impact the long-term behavior of the cell cycling network when Let-7 couples the two GRNs.
To determine if the cell cycle network can modulate the dynamics of the transformation network, we simulated overexpression of all cyclins by increasing their synthesis rates, V SMCYC . Cyclin overexpression promotes uncontrolled cell proliferation of cancer cells (Gillett et al., 1994;Pok et al., 2013). Starting from a nontransformed cell state (high Let-7 levels), the model shows a stable down-regulation of Let-7 when V SMCYC increases ( Fig. 4A where V SMCYC increases at t > 100h). The corresponding temporal evolution of Let-7 and cyclin B/Cdk1 is represented for different synthesis rates of cyclins (at t > 100h, V SMCYC changes from 1.5 to 2.5 in Fig.   4B and from 1.5 to 12 in Fig. 4C). For weak overexpression of cyclins, the non- Temporal evolution of Let-7 and cyclin B/Cdk1 indicates that the switch to cell proliferation triggered by cyclin overexpression is irreversible because downregulation of cyclins to their initial levels will not restore cell cycle arrest ( Supplementary Fig. 3A, condition 1, 1200h < t < 1500h). This is the consequence of the irreversible bistable switch at the core of the transformation GRN . The model predicts that a stronger decrease in cyclin synthesis can eventually (1) the cholangiocarcinoma cohort, which displays the highest CPI levels; (2) the hepatocellular carcinoma cohort, which shows high CPI and low NSTI values, and (3) the prostate adenocarcinoma cohort characterized by low CPI and high NSTI. This analysis revealed that non-tumor (blue dots) and tumor samples (red dots) cluster separately in cholangiocarcinoma ( We concluded that the network activation varies from patient to patient and depends on tumor-type. Thus, an increase in tumor suppressors or a decrease in oncogenes reduces the probability of stochastic switches to a transformed, proliferative, cell state.

Dynamics of cell proliferation and cell transformation of a heterogeneous
However, if cells are already in a proliferative, transformed state, similar changes in tumor suppressors or oncogenes do not permit to revert back to a more "healthy" cell phenotype, which highlights an irreversible process in cancer progression.

Discussion
Tumorigenesis rests on many biological features which include sustained proliferative signaling, evading growth suppressors, resisting cell death, promoting angiogenesis, ensuring replicative immortality and eliciting invasion and metastasis (Hanahan and Weinberg, 2011). Here, we built a mathematical model to analyze the dynamical properties of a Let-7-dependent mechanism coupling cell proliferation and an epigenetic switch driving malignant transformation.
Our mathematical model illustrates qualitatively how Cdk-dependent and transformation networks may interact, and proposes a mechanism, acting through Let-7, which suggests that cyclin overexpression can promote cell proliferation, while inducing and accelerating malignant transformation. Indeed, overexpression of cyclins progressively sponges the free form of Let-7. The latter will no longer be available to repress the components of the transformation network, leading to the activation of the epigenetic switch. This effect is known as competing-endogenous, ceRNA, effect. CeRNAs regulates other RNA transcripts by competing for shared miRNA and were involved in tumorigenesis (Salmena et al., 2011;Tay et al., 2014;Chiu et al., 2017). Here Let-7 is the shared miRNA between both networks. Let-7 was shown to be involved in different ceRNA mechanisms. Indeed, Let-7e can modulate the inflammatory response in vascular endothelial cells through a ceRNA effect (Lin et al., 2017). Imprinted H19 lncRNA, which plays important roles in development, cancer and metabolism, modulates Let-7 availability by acting as a molecular sponge and causing precocious muscle differentiation (Kallen et al., 2013).
Moreover, amplification of MYCN mRNA levels in neuroblastoma can sponge Let-7 thereby rendering LIN28B dispensable for cancer progression (Powers et al., 2016).
Note however that experimental and theoretical studies indicate that a ceRNA effect between multiple RNA transcripts and the shared miRNA is effective only in the presence of adequate expression levels of the transcripts and the miRNA (Gerard and Novak, 2013;Denzler et al., 2014;Tay et al., 2014).
Our mathematical model further shows that random variations in gene expression from cell to cell can create large fluctuations in the global network dynamics leading to stochastic switches of some cells to a transformed and proliferative state. These stochastic switches in the GRN dynamics could be a source of heterogeneity in cancer cell populations (Tang, 2012;Patel et al., 2014). Stochastic switches in gene networks were identified in hematopoietic tumor stem cells (Dingli et al., 2007), in the appearance of mammary tumor in mice (Bouchard et al., 1989), and in the differentiation and maturation of T lymphocytes (Davis et al., 1993). These switches may also give rise to an equilibrium in population of cancer cells (Gupta et al., 2011). Transcriptomic analysis of TCGA data suggests that the coupling between the cell cycle and a malignant transformation networks, and the activation of these networks in tumors are cancer type-specific, with predominant activation in gastrointestinal cancers. Our PCA analysis reveals inter-patient heterogeneity in network activation in tumors, which stresses the need to consider patient-specific characteristics when optimizing therapeutic strategies by reversing network dynamics of activated GRNs (Biankin et al., 2015).
In conclusion, by means of transcriptomic data analysis and modeling-based investigations, we identified a Let-7-dependent connection between two major GRNs involved in tumorigenesis, and whose activation is cancer-and patient-specific. We anticipate that a better characterization of the dynamics resulting from the combination of other GRNs specific for each patient will help providing a global GRN activation map for personalizing and optimizing cancer treatment.    passes from 0.1 to 0 for 500h < t < 1000h), (C) NF-kB inhibition (k AA1NFKB = k AA2NFKB = k AA3NFKB = 0 for 500h < t < 1000h), or (D) PTEN over-expression (V SMPTEN passes from 0.001 to 10 for 500h < t < 1000h). Other parameter values are as in Supplementary Table 3.               Table S1, the kinetics equations are found in Table S2, while the description of the parameters and their numerical values used in the simulations are in Table S3.
In the model, we consider for simplicity a constant total concentration of NF-kB

Analysis of RNASeq data
RNASeq and miRNASeq data of the different patient cohorts were from TCGA database (http://firebrowse.org/). For each cohort, we converted the "scaled_estimate" in the "illuminahiseq_rnaseqv2_unc_edu_Level_3_RSEM_genes" file into TPM by multiplying by 10 6 .  Table S2. Kinetic equations of the mathematical model

Model for the inflammatory circuit driving the epigenetic switch to cell transformation
Coupling the Cdk network and the epigenetic switch through Let-7 miRNA Where: Note: We focus our study on qualitative rather than quantitative aspects in the expression of the different components of the network. Indeed, our goal is to analyze the dynamic implications of the regulatory structure of the epigenetic switch linking inflammation to cell transformation and of the Cdk network driving the cell cycle, i.e. the wiring diagram, or the topology, which is crucial for the network behavior (see (Wagner, 2005)). Moreover, many parameters have not been determined experimentally. The numerical values have been selected to yield a non-transformed state without the inflammatory signal, Src; or a transformed state within less than 100h with transient inflammatory signal, which corresponds to experimental observations . For the Cdk network, parameter values have been chosen to reach cell cycle periods ranging between 20h and 50h, which correspond to the doubling time in most mammalian cells (Alexiades and Cepko, 1996). The time units for the parameters in the model are expressed in hours, while the concentrations are expressed in µM. In addition, the mechanisms controlling miRNA degradation are complex and not fully understood (Zhang et al., 2012). As a consequence, we assume that each complex between a mRNA and its regulating miRNA is rapidly targeted for degradation.   Figure S1. Detailed scheme of the model. The model for the epigenetic switch linking inflammation to cell transformation originates from a previous study , while the model for the post-transcriptional regulation of the Cdk network by Let-7 is based on a skeleton model for the Cdk network driving the mammalian cell cycle . Let-7, which is involved in multiple positive feedback loops, can repress the synthesis of each cyclin.  Table S3. Temporal evolution of Let-7 and cyclin B/Cdk1. From a non-transformed and quiescent cell state, V SMCYC changes from 1.5 to 7 (at t > 100h), and promotes the epigenetic switch and cell proliferation.
For 1200 h < t < 1500 h, V SMCYC is set back to its initial value (= 1.5), cell proliferation and the transformed cell state are maintained, which is characterized by sustained oscillations in cyclin B/Cdk1 and low levels in Let-7. In A, at t = 1500h, V SMCYC changes from 1.5 to 0.1, preventing cell proliferation without recovering a non-transformed cell state since cyclin B/Cdk1 tends to a low stable steady state level, while Let-7 levels remain low. In C, for 1500h < t < 2000h, the synthesis rate of Lin28, V SLIN28 is set to 0, allowing the recovery of a stable non-transformed state defined by high levels of Let-7 and impeding cell proliferation, characterized by a low, stable steady state, level of cyclin B/Cdk1. In C, for t > 2000h, basal conditions are used. (B, D) Temporal evolution of Lin28 and STAT3 corresponds to conditions of panels A and C, respectively. Other parameter values are as in Table S3, which correspond to basal conditions.