Systems medicine in colorectal cancer: from a mathematical model toward a new type of clinical trial.

Current colorectal cancer (CRC) treatment guidelines are primarily based on clinical features, such as cancer stage and grade. However, outcomes may be improved using molecular treatment guidelines. Potentially useful biomarkers include driver mutations and somatically inherited alterations, signaling proteins (their expression levels and (post) translational modifications), mRNAs, micro-RNAs and long noncoding RNAs. Moving to an integrated system is potentially very relevant. To implement such an integrated system: we focus on an important region of the signaling network, immediately above the G1-S restriction point, and discuss the reconstruction of a Molecular Interaction Map and interrogating it with a dynamic mathematical model. Extensive model pretraining achieved satisfactory, validated, performance. The model helps to propose future target combination priorities, and restricts drastically the number of drugs to be finally tested at a cellular, in vivo, and clinical-trial level. Our model allows for the inclusion of the unique molecular profiles of each individual patient's tumor. While existing clinical guidelines are well established, dynamic modeling may be used for future targeted combination therapies, which may progressively become part of clinical practice within the near future. WIREs Syst Biol Med 2016, 8:314-336. doi: 10.1002/wsbm.1342 For further resources related to this article, please visit the WIREs website.

INTRODUCTION

At the molecular level: colorectal cancer (CRC), as other cancers, is a multihit disease due to the accumulation of genetic mutations/alterations (inheritable at the somatic level) in genes that control cell growth, programmed death, differentiation, and cell‐tissue architecture.1 The last decade has seen the introduction of new methods for high throughput DNA and RNA analysis of sequence, copy number and expression variations, as well as somatically inheritable (passed through cancer cell generations) epigenetic alterations (particularly in methylation of promoter region CpG islands), which has allowed for the expansion of our knowledge of the alterations involved in CRC.2, 3 Significant tumor heterogeneity between patients and even within the same tumor, has been identified with cancer biology, but its significance for improved clinical management has only been more recently understood.4, 5, 6

The study of all phenomena causing the onset and dynamic behavior of the disease, through a systems biology approach, has led us to an understanding of cancer as a disease of signaling networks of the cell. We can subdivide a network into interconnected ‘pathways.’ Different individual cancers tend to be more similar in terms of altered signaling pathways, than in terms of the individual mutations possibly present in a given pathway.3 Within a given pathway mutations can be, at least partially, interchangeable and mutually exclusive.7

Therapeutic agents targeting individual specific signaling‐proteins (and de facto their pathways), is one of the major applications of personalized medicine in oncology. Targeted therapies have been (and are being) developed, and some have already been shown to confer benefits.8 More will become available during the next 5–10 years, as not only individual targeted drugs, but also (more interestingly) as combinations of distinct targeted drugs attacking distinctly altered pathways. Clinical application of personalized medicine already depends on the identification of relevant biomarkers for early cancer detection, prognostic markers for risk stratification, and predictive markers associated with response to a particular therapy. Biomarkers already play a significant role in the management of patients, and intense research efforts aim at identifying new promising molecular markers (DNA, RNA, or protein biomarkers). When identifying ideal therapies, considering both experimental and validated targeted drugs (especially ones targeting signaling‐proteins affected by excess of function), the reconstruction and mathematical dynamic modeling of the signaling‐networks (or regions of them) can be an especially useful and superior method for biomarker integration. Reasoned suggestions of targeted‐drug combinations are the natural output of an integrated approach of this kind. We also discuss the need for new innovative clinical trial designs, to verify the efficacy of combinations of targeted drugs, a goal difficult to achieve applying conventional study designs.

ROLE OF BIOMARKERS IN CRC, RESEARCH, AND CLINICAL PRACTICE

Advances in molecular biology in the last decades have helped to elucidate some of the genetic mechanisms leading to CRC. Biomarkers have begun to play an increasingly important role in the management of patients, and intense research efforts aim at identifying new promising molecular markers (at the DNA, RNA, or protein levels).9, 10, 11, 12

Biomarkers are measurable characteristics or factors that ‘could be indicators of normal biologic processes, pathological processes, or pharmacological responses to therapeutic interventions.’13 In cancer, a biomarker is ‘a molecular, cellular, tissue, or process‐based alteration that provides indication of current, or more importantly, future behavior of the cancer.’14

Clinics already utilize biomarkers, including driver mutations (both at DNA and protein levels), mRNAs, micro‐RNAs, and other somatically inheritable alterations. Oncologists use biomarkers for screening, risk stratification, prevention, diagnosis, treatment, and surveillance.15

We will focus our attention mostly on a limited group of bio‐medically relevant biomarkers related to diagnosis, prognosis, and treatment.

Driver Mutations and Other Somatically Inheritable Alterations as Biomarkers

Important gatekeeper genes, involved in CRC, are repressor genes playing a double role. They are involved with both altered networks (see Gatekeeper and Caretaker Mutations section) and in determining genetic predisposition. Some important ones include: the Familial Adenomatous Polyposis gene (APC)16; the Juvenile Polyposis genes (SMAD4, BMPR1A that belongs to the family of TGF‐β receptors)17; the Cowden syndrome gene (PTEN)18; and the Peutz‐Jeghers syndrome gene (STK11/LKB1, required for the organization of cell polarity).19 TP53 and CHK2 (Li‐Fraumeni syndrome) may also confer increased CRC risk and be involved in an altered signaling‐network.20

Mismatch repair genes are caretaker genes whose loss or reduction of function is also associated with CRC predisposition, because they decrease DNA‐replication fidelity. They are connected with the Lynch syndrome (HNPCC—hereditary nonpolyposis CRC), an autosomal dominant genetic condition. Both immunohistochemistry and genetic testing can be performed.21 For biomarker sensitivity and specificity, and a rather extensive description of the molecular diagnostics of CRC, see Refs 22 and 23, respectively.

Evidence from clinical and retrospective studies allowed for the discovery of, and introduction into clinical practice, several predictive molecular biomarkers for the identification of patients suitable for targeted therapies, and thus prevented unnecessary side effects from a‐specific therapies.11

Guidelines routinely recommend evaluation of KRAS gene mutations in patients with metastatic CRC, to predict the response to antiepidermal growth factor receptor antibody‐based therapies, such as cetuximab and panitumumab.24, 25, 26 Recent studies have reported evidence that, in addition to KRAS, mutations in NRAS predict nonresponse to anti‐EGFR therapy, and should therefore be tested for predictive purposes.9 Several studies suggest that we can also use BRAF mutations as predictive markers for ineffective EGFR‐targeted therapy.9 This similar behavior is rationally plausible because KRAS and BRAF are ‘adjacent’ in the same pathway. We rarely find a mutated KRAS and BRAF in the same cancer specimen (except in the case of distinct subclones: intra tumor heterogeneity). The National Comprehensive Cancer Network recommends the assessment of the presence/absence of wild‐type KRAS, NRAS, and BRAF, prior to deciding which treatment strategy is optimal for individual CRC patients.27 Investigators have also evaluated mutations in PIK3CA and loss of PTEN expression, as co‐predictive markers for inefficient anti‐EGFR therapies9 (they are in a distinct pathway downstream of EGFR, see Figure 2); these alterations act via AKT phosphorylation and activation; AKT1 could also be directly mutated and activated at low frequency in CRC (E17K). This pathway flows towards mTOR activation. The KRAS (MAPKs) pathway and PIK3CA pathway behave as independent synergic pathways. We observe them altered together at a frequency greater than expected by chance.28

Figure 2

Molecular Interaction Map (MIM) including most of the G0–G1 transition in CRC. We surrounded the cartouches (like in hieroglyphics), of mutated/altered signaling‐proteins in the HCT116 line, with a light yellow oval. The oncoprotein inhibitors are without cartouches and an arrow indicates their relative targets. We enlarged this new MIM by about 50% in respect to the MIM we published in 2015.90 Numbers near the interaction lines refer to an annotation list of interactions (supplementary material of our 2015 article90). Grid coordinates (numbers at the left and letters on top of the MIM) help to locate molecular species and interactions present in the MIM.

Recent studies highlighted that circulating DNA, from blood samples, could be convenient for analysis of tumor mutations in real time. Liquid biopsies can identify potentially clinically relevant mutations that were not detected in the primary tumor tissue at the time of biopsy.29, 30 Liquid biopsies can be very useful for the early detection of a relapse, and can be used in monitoring the molecular evolution of tumors in response to targeted therapy.29 Potentially, this could also be used for early detection in patients known to be at high risk for certain cancers.

As mentioned above, mutations in genes involved in the DNA mismatch repair system (MLH1, MSH2, MSH6, and PMS2) result in alterations in highly repeated DNA sequences (microsatellites). Actually, MSI is a diagnostic marker for Lynch Syndrome (HNPCC, hereditary nonpolyposis CRC), an inherited CRC syndrome. Clinical trials and retrospective studies have also reported that patients with MSI tumors do not benefit from 5‐fluorouracil (5‐FU) adjuvant systemic chemotherapy, while patients with microsatellite stable tumors do. In addition, some studies have also found that MSI lines were more sensitive to SN‐38 (irinotecan, an inhibitor of topoisomerase‐1) than MSS lines. The value of using MSI as a predictive marker for these chemo‐sensitivities remains controversial and is still under evaluation.9

Proposed Panels of mRNAs

Mesenteric lymph node (ln) metastasis is the single most important prognostic characteristic in CRC. The ln status, used for staging, is a decisive selection criterion for postoperative adjuvant therapy. In addition to histopathology, a study of 14 biomarker mRNAs (using RT‐PCR on mesenteric lns material collected surgically) showed that high CEA, low MUC2 and a high KLK6 mRNA levels together formed a strong ‘trio’ for staging and prediction of outcome in CRC patients.31

Marshall et al. tested RNA extracted from peripheral blood cells and created a test based on a seven‐gene biomarker panel (ANXA3, CLEC4D, LMNB1, PRRG4, TNFAIP6, VNN1, and IL2RB). The authors derived the panel from a 196‐gene expression profile considering 112 CRC patients (including those with stage  I, II, III, and IV disease) and 120 controls. The panel was to some extent validated with material from an additional 202 CRC patients (at all disease stages), and 208 controls. Marshall et al. reported 72% sensitivity and 70% specificity for this initial study. They also developed a commercial blood test: the ColonSentry® seven‐gene mRNA biomarker panel. The test is supposed to determine the risk of developing cancer; with approximately 70% sensitivity and specificity it is not ideal for clinical purposes32; however, sensitivity is similar to other well‐established screening tools. Arguably, it could be used as an indicator for further screening.

Oncotype DX® Colon is a diagnostic test for newly diagnosed patients with stage II and III colon cancer. The test works by examining the tumor tissue at a molecular level, in order to provide information about the individual biology of each patient's tumor. To do this, the test evaluates 12 genes (expression levels measured through RT‐PCR) within the colon tumor tissue, to determine the likelihood that the cancer cells will spread, or metastasize, within 3 years after diagnosis. The seven cancer‐related genes included in the test are: BGN, INHBA, FAP, genes associated with activated stroma; MK167, MYBL2, MYC, which are associated with cell cycle activation; and GADD45B (which is related to genotoxic stress). ATP5E, GPX1, PGK1, UBB, and VDAC2 are used as reference control genes.33, 34 Oncotype DX® Colon combines the measurement of these 12 genes into an individualized result called the Recurrence Score®. The Recurrence Score® is a number between 1 and 100 that correlates with the likelihood that an individual patient's colon cancer will recur. This information can help healthcare providers, and their patients, to make informed decisions about whether or not to implement additional chemotherapy treatment following surgical removal of the primary tumor.

Oncotype DX® Breast35 is a 21 gene panel (also using mRNA expression) for breast cancer. This panel is at a more advanced stage of clinical validation and approved for reimbursement by the FDA in the United States.

Notice that these panels try to predict the risk of recurrence. Identifying the most effective combination of targeted drugs against pathways altered in individual cases is not addressed with these existing panels.

At a more basic level, some investigators have suggested that mRNAs of genes belonging to the same pathway could have relatively similar average levels of expression.36 Consequently, ranking thousands of normalized mRNAs expression levels, we should observe a prevalent statistical clustering (‘Enrichment Score’) of (one to two dozen) mRNAs (from the same pathway) at the top, bottom or middle of our general ranking, in a nonrandom fashion. This is the key principle of Gene Set Enrichment Analysis (GSEA). Microarray data sets referred to independent experiments, on similar biological material, display a more homogeneous behavior when analyzed using GSEA.36 This approach has improved the interpretation of mRNA expression levels.

Noncoding RNA Transcripts: microRNAs and Long Noncoding RNAs

MicroRNAs are short RNA sequences, about 18–24 bases in length. MicroRNAs can modulate the function of either tumor suppressors or dominant oncogenes, modulating the transcription or translation rates, or the half‐life, of important CRC mRNAs and/or proteins.37

They can themselves have deregulated levels of expression or base changes (especially SNPs— single nucleotide polymorphisms). SNPs can also potentially change the affinity between microRNAs and microRNAs binding sites.38, 39 Even for SNPs related to microRNAs, we can have ‘passenger irrelevant changes’ and ‘driver alterations.’

Several microRNAs appear consistently altered in CRC, of which roughly two third are upregulated, and one third are downregulated.40 In fact, miR‐143 is classified as a tumor suppressor, and numerous studies have shown that downregulation of some microRNAs (including miR‐342, miR‐143, and miR‐145) plays a role in colon tumorigenesis, while miR‐21 appeared to be elevated in CRC in at least seven studies.40, 41, 42 These microRNAs can also be coded for in DNA regions where copy‐number gain or loss is detected, and microRNA promoter methylation has been observed in neoplastic tissue and is believed to play a significant role in CRC. Interestingly, miR‐342 has been shown to play a role in altered epigenetic patterns, correlating to hypermethylation, as observed in aggressive colon carcinomas.42, 43 High methylation levels observed clinically in patients with poor prognosis have been shown to be causative in at least animal models.43 Functional miR‐342 downregulates DNMT1 (the DNA methyltransferase isoform responsible for copying DNA methylation patterns throughout cell replication), while loss of miR‐342 allows for the increased DNMT1 expression (and thus copying of the denovo methylation mediated by DNMT3b in aggressive CRC), and is correlated to hypermethylation of tumor suppressor gene promoters (and their resultant silencing).42, 43

MiR‐135b can repress the expression of APC,44 while MiR‐143 represses the translation of KRAS thereby attenuating the effect of mutated KRAS.45 This could be significant for anti‐EGFR therapies.

MiR‐34a appears to inhibit cell cycle progression by targeting the FRA1/FOSL1 (FOS genes family).46 FOS inhibition means inhibition of FOS:JUN dimerization (AP1 complex, see also Figures 2 and 3).

Figure 3

Simplified scheme of a model promoter region regulating MYC and CCND1 transcription. The model promoter region consists of a number of important Transcription Factor Binding Sites, depicted in arbitrary order. Each TF considered to bind to its promoter is a component of the MIM. Schematic representation of the promoter region, with activators (a) and repressors (b) assumed to bind to each TFBS in the model. Arrows indicate the potential for binding at the TFBS. We also show the DNA binding regions of activators (white) and repressors (black). [Reprinted with permission from Ref

MiR‐365 can target and contrast the repressor activity of RB, favoring the crossing of the restriction point between G1 and S cell cycle phases.47, 48

Altered expression levels of some microRNAs in liquid biopsies (for instance miR‐92 and miR‐141) can be of diagnostic interest.49, 50

Measuring microRNAs from stool can also be of diagnostic interest; this approach has been validated for elevated miR‐21 expression and the hyper‐methylation of miR‐34b/c promoter regions.51, 52, 53

High miR‐21 expression appears to be associated with poor prognosis and the worst therapeutic outcomes, and this appears to be true for several cancer types.54, 55, 56, 57, 58, 59, 60, 61, 62

In the case of excess of microRNA function, we can think of the possibility of microRNA inhibition using antisense nucleotides. In the case of microRNA loss of function, we can envisage microRNAs replacements.37, 63 The field of microRNAs is less developed than that of driver mutations/somatically inheritable alterations, but its integration into a general more comprehensive picture could be of outmost relevance.

Investigators define ‘long noncoding RNAs’ as non‐mRNAs longer than 200 nucleotides; they are in principle deprived of an open reading frame (ORF). These longer noncoding RNA transcripts can exert regulatory functions interacting with other RNAs, with proteins or specific DNA sequences.

They can sometimes interact with different targets in a modular way. There is room for multiple specific interactions and specific modulations of function.64 The mammalian genome encodes many thousands of noncoding transcripts including both short (<200 nucleotides in length) and long (>200 nucleotides) transcripts. For some of them, low levels of expression and low levels of evolutionary conservation introduce the suspicion that we may sometimes be dealing with transcriptional noise, but one must also wonder why it has been maintained evolutionarily given the obvious energy expenditure. It may be that maintenance of this ‘noise’ is outweighed by the benefit of variation, or that they indeed hold a yet to be discovered function.

A chromatin signature, consisting of a short stretch of histone protein H3 tri‐methylation at the lysine in position 4 (H3K4me3), which corresponds to promoter regions, usually accompanies significant transcription. A longer stretch of tri‐methylation of histone H3 at the lysine in position 36 (H3K36me3) covers the entire transcribed region.65, 66 This chromatin signature can be present not only for mRNAs, but also for noncoding RNAs.67 ncRNAs co‐expressed with the mRNAs of a given pathway can sometimes be given a guilt‐by association role, at least as a working hypothesis. RNA interference can help to give a role to long intervening noncoding RNAs (lincRNAs).68

We have mentioned large ncRNAs particularly to make the reader aware of the multiple new levels of complexity that will undoubtedly emerge, even for CRC, during the next few years.

Clinical Practice Guidelines

In addition to basic, translational, and clinical research levels, we have present day clinical practice guidelines.

Management of CRC patients has slowly improved during the last 30 years; minor variations in cancer management and outcomes still exist amongst industrialized countries. National and international organizations provide clinical practice guidelines (CPGs) to support clinicians.

CPGs are defined as ‘systematically developed statements to assist clinicians and patients in decisions about appropriate health care for a specific clinical circumstance.’69 The guidelines include recommendations for treatment and management decisions to improve the quality, effectiveness, and efficiency of patient care. In the oncology field, the heterogeneity of cancer diseases and the complexity of therapeutic decisions (for instance in the case of comorbidities in elderly cancer patients), have led to the development of guidelines to ensure conformity of practice, promote conformity in the delivery of medical care and to develop minimal required standards of good clinical practice by quality assurance.

Guidelines base their recommendations on the systematic review of relevant medical literature and expert input at the time of their generation. Updating and revision is a continuous process, reflecting incoming data and new clinical information that may affect clinical practice standards. The guideline process not only involves development but also validation, dissemination, and assessment of impact.

Collaborative groups of multidisciplinary experts are required to define consensus guidelines for clinical practice and several initiatives are available from the Oncology Community, with the goal that these international efforts will eventually substitute preexisting national and local guidelines. A large number of clinical guidelines are available, according to tumor types, aimed at cancer prevention, diagnosis, assessment of treatment outcomes, as well as for supportive palliative care.

The guidelines of the National Comprehensive Cancer Network (NCCN) represent an alliance of 26 of the world's leading cancer centers in the United States. A study published in 2013 validated 2012 NCCN colon cancer practice guidelines.70 Statistical analysis documented a significant survival benefit for patients who received treatments adhering to NCCN guidelines.26

European Society for Medical Oncology (ESMO),25 The American Society of Clinical Oncology (ASCO),24 and The European Registration of Cancer Care (EURECCA) have the same goals as the NCCN.71 A variety of scientific societies also curate and regularly update evidence‐based guidelines for the use of biomarkers in CRC.32, 72, 73

The American Society for Clinical Pathology, the College of American Pathologists, the Association for Molecular Pathology, and the American Society of Clinical Oncology produced New Guidelines on Colorectal Cancer Molecular Testing in 2015, and included KRAS, extended RAS, BRAF, and dMMR/MSI as recommended molecular markers.74

It is important to distinguish the stage of an established good clinical practice guideline from the stage of innovative clinical trials, and goals projected toward in coming years. Technological and conceptual evolution is roughly 5–10 years ahead of routine clinical practice. In these developmental phases, costs are an additional issue.

Artificial Intelligence decision support systems, including ones developed by US Cancer Centers with IBM (under the Watson eponym), have obvious relationships with more classical clinical guidelines. They could overtake them in the long run, at least in some ways.75

An authoritative advocate for lighter regulations, allowing for more rapid introduction of discoveries into clinical practice, is Vincent DeVita. In his last book (‘The death of cancer’),76 he discusses this important issue. The most important innovations are linked to the introduction of more targeted drugs correcting more altered pathways, and something must be done to accelerate this process.

REASONED COMBINATIONS OF TARGETED INHIBITORS OF ONCO‐PROTEINS AFFECTED BY EXCESS OF FUNCTION: A DYNAMIC SIGNALING‐NETWORK MODEL

We focus on future goals aiming towards personalized cancer combination therapy, using drugs that are targeted inhibitors of onco‐proteins affected by excess of function in the context of a signaling‐network of (bio) chemical interactions. Selective inhibitor drugs are chemicals that respond to the general laws governing chemical interactions, of which biochemical interactions are a direct subset. One can easily introduce them into a dynamic model.

We present in Figure 1 a flow chart of the different phases of construction and implementation of our dynamic model of biochemical interactions/reactions among signaling proteins (G0–G1–S cell cycle phase).

Figure 1

Flow chart of the different phases of implementation of our dynamic model.

To avoid the explosion of a variety of specificities related to different cancers, we will refer, as an exemplary model especially for solid carcinomas, to CRC. We are aware of the fact that, especially in solid carcinoma tumors, the architecture of signaling‐network subregions and their pathways come back repeatedly.1 At least some of the driver and gatekeeper mutations tend to be present (albeit with frequency variations) in a multiplicity of them.77, 78, 79

An additional restriction we introduced was the concentration of our attention to the G0–G1 cell cycle phase, when the cell commits to the irreversible decision of entering the S phase, with decisive subsequent consequences on cell replication. The abundance of driver/gatekeeper mutations in the signaling‐network subregion involved in this cell cycle phase suggests it is of crucial relevance.48, 80, 81

Gatekeeper and Caretaker Mutations

To understand the mechanisms of preneoplastic lesions (benign polyps in the colon) and subsequent tumor emergence and evolution, we need to identify the genes that drive tumorigenesis. All genomes contain thousands of somatic mutations (often called passenger mutations), but only a few of them ‘drive’ a normal cell to evolve into a cancer cell, by altering genes which confer selective growth advantages to tumor cells.1, 82, 83 The Vogelstein group in Baltimore refers to a basic fundamental distinction between gatekeeper and caretaker genes.1, 84, 85 The co‐presence of a discrete number of somatically inheritable alterations sustains the actual malignancy of a cancerous cell: each gatekeeper alteration affects a ‘gate’ in a biochemical pathway of the signaling network. This discrete number of ‘hits’ (let say 5–10 hits) is required to make the cancer cell malignant when we clinically detect the tumor, independently from its past evolutionary history.

How long a cell cluster takes to evolve to a malignant tumor depends largely on the ‘caretaker mutations’ it acquires. These driver mutations reduce the fidelity of cell replication, both in terms of DNA bases and chromosomal fragmentation/rearrangement. Defects in caretaker's ability to preserve proper DNA sequence can increase the base‐level mutational probability of driver/gatekeeper mutations, while chromosomal instability/fragmentation may result in copy number decreases in gatekeeper repressor genes, or increases in gatekeeper dominant oncogenes and chimeric oncogenes. Caretaker level mutations will facilitate/accelerate the appearance of errors at the gatekeeper level.1, 84, 85 We can see an altered caretaker level as having a sort of relativistic effect of time‐scale compression on cancer evolution.

Recently, the importance of multiple inheritable chromatin alterations86, 87, 88 has emerged, they can affect the transcription of large and a‐priori nonspecific sets of genes. Inside these sets, we can have crucial up or downregulation of expression of driver/gatekeeper genes.

Somatically inheritable epigenetic alterations can also affect driver/gatekeeper genes, as is seen with CpG methylations in their promoter regions.89 MicroRNAs37 and other nontranscribed RNAs64 can also affect the levels of specific genes.

G0–G1 Cell Cycle Phase and the G1–S Restriction Point

There are three main waves of transcription in a mammalian cell cycle: the G1 to S transition, (the most studied cell cycle phase and the focus of this section), the G2 to M transition and the M to G1 transition.48

During the G1 to S transition, phosphorylation of transcriptional inhibitors (RB + proteins of the pocket family) by cyclin‐dependent kinases (for instance CyclinD:Cdk4 and CyclinE:Cdk2) releases them from transcription factors (E2F1‐3/DP1‐2 for instance). The consequence is activation of G1–S gene transcription (as depicted in our molecular interaction map (MIM; see Figure 2) and supplementary material of Tortolina's 2015 Oncotarget paper90). In cooperation with other transcription factors (Figure 3), E2F1‐3:DP1‐2 act to transcribe dozens of genes, in different cellular contexts and differentiation conditions.91 E2F1‐3:DP1‐2 targeted genes can be positive cell cycle gene regulators, and include CCND1, CCND3, JUN, MYC, MYCN, CCNE1, CCNE2, CDC25A, CDK2, E2F1‐3, NPAT, MYB, MYBL2, and TFDP1. Among the negative cell cycle gene regulators targeted are: CDKN1C, CDKN2C, CDKN2D, E2F4‐8, RB1, RBL1, and TP53. This wave of transcription initially reinforces a positive feedback loop, further activating G1–S transcription. CyclinE:Cdk2 further phosphorylates the pocket proteins, favoring the transcription of more Cyclin D/E, generating an initial positive feedback loop.

Negative feedback loops (including E2F4‐8:DP2 activation) terminate the first wave of gene expression at the transition from G1 to S phase.48 An example of a negative feedback loop is the transcription of E2F6‐8 downstream of E2F1‐3:DP1 transcription factors. This happens shortly before the phosphorylation (P‐) of E2F1‐3 by CyclinE:Cdk2. P‐E2F1‐3 detach from their promoters. We can consider the combination of the two effects (phosphorylation and detachment) as a negative feedback loop, because E2F6‐8 now can repress transcription at the same promoters previously activated by E2F1‐3. Another target of E2F1‐3 transcription is the SKP2 protein, involved in an ubiquitin ligase pathway degrading E2F1‐3,48 an additional component of the negative feedback loop.

Transcription of multiple G1‐S cell cycle genes is irreversibly committed to at the G0→G1 transition point, just before the earliest newly transcribed genes required not just by quiescent‐cell maintenance, but by new cell replication.

Mammalian cells (and yeast) have this irreversible ‘restriction point,’ after which the cell is committed to the S and subsequent phases, independently of signals from the environment.

De‐repression of G1‐S transcription allows progression into S phase in an unrestrained fashion. When this happens because of a discrete number of mutations/hits in the signaling‐network upstream, one can consider this multihits event a hallmark of cancer. Mutations/alterations in G0‐G1 signaling‐proteins (see our MIM pathways upstream of MYC/CCND1 transcription, Figure 2) can precisely disrupt (through a discrete number of hits) this 'irreversible restriction point.' The gate(s) to S phase and further cell cycle progression will remain open. This aspect is crucially relevant for cancer.

A Dynamic Signaling‐Network Model

We modeled a network subregion related to the G0–G1–S cell cycle transition, downstream of TGFβ, WNT, HGF, and EGF‐family receptors. We also included Extracellular Matrix Integrin receptors and the ALK receptor (see Figure 2). The new MIM is about 50% larger than the 2015 MIM we presented in Oncotarget.90 The signals downstream of the membrane receptors propagate through a complex network, involving cross talk amongst interacting pathways, and strong feedback loops on different levels.

Our MIM is especially relevant because it represents the network (multiple pathways) just preceding a crucial restriction point concerning irreversible decisions about cell replication. Driver gatekeeper mutations affecting the pathways reconstructed in our MIM are quite frequent (but not exclusively) in CRC: they open closed gates precisely during this phase of the cell cycle. In perspective, in our MIM reconstruction, we move in the direction of an extensive representation of the crucial G0–G1–S transition, and therefore we should find in it most of the cancer gatekeeper mutations.

A MIM is a diagram convention that is capable of unambiguous representation of networks containing multiprotein complexes, protein modifications, and enzymes that are substrates of other enzymes. This graphical representation makes it possible to view all of the many interactions in which a given molecule may be involved, and it can portray competing interactions, which are common in bio‐regulatory networks.92, 93, 94, 95, 96 Alternative syntaxes have also been proposed.97, 98

Our dynamic model suggests that the duration of the initial G0‐G1 cell cycle phase is short (30–60 min). It is made mostly of biochemical interactions and posttranslational modifications, e.g., phosphorylation→de‐phosphorylation, as our MIM also shows. It is probably a preparatory time where the cell focuses only on future transcription/translation of new proteins. At an experimental level, FBS (fetal bovine serum) stimulation of a quiescent mouse epidermal cell line resulted in AKT and GSK‐3β phosphorylation after 15–30 min.99 Similar times for MAPK phosphorylation have also been reported.100

The last 4 h of the G1 phase involve the transcription and translation of new proteins required for entering the S phase (we have mentioned some of them above, in G0–G1 Cell Cycle Phase and the G1–S Restriction Point section). Because of a multiphosphorylated RB1, in cooperation with other pro‐transcriptional G0–G1 pathways, additional proteins will appear in the early S phase, including proteins required for DNA duplication and finally the duplication of the components of the entire cell mass.48, 91

In the G0–G1 cell cycle phase that we have reconstructed, we do not have a transit through multiple waves of subsequent new transcriptions/translations and accompanying protein degradations, thus the reduced size of the MIM makes data handling simpler (while simultaneously sufficiently detailed) to interrogate interactions at a biochemical interactions level.

We moved ahead from the stage of a static signaling‐proteins interaction descriptive map (MIM), to dynamic mathematical modeling. Our reactions’ dynamic simulations used ODEs. As of the end of 2014, they involved 660 reactants (basic and modified species, complexes and inhibitors), 348 reversible reactions and 174 catalytic reactions, for a total of 870 reactions processed with one tool (348 × 2 + 174 = 870 reactions),90 while by June of 2015 we had increased the number of reactions modeled by almost 50%.

Notice that, from a basic perspective, chemical interactions/reactions, are at the very basis of the life of any cell, and therefore of life itself. For dynamic mathematical simulations, efficient software is available for achieving numerical solutions of ordinary differential equations (ODEs). We can achieve a satisfactory representation of the temporal evolution of large networks of biochemical interactions, and used Matlab's Simbiology software to do so.101 After the introduction of a mutated function, for instance through an altered/mutated signaling‐protein, Simbiology is able to converge rapidly to a new equilibrium (offering a numerical solution).

We have submitted our model as SBML to BioModels,102, 103 submitted by Dr. Vijayalakshmi Chelliah on behalf of Dr. Lorenzo Tortolina. MODEL1601250000.

A model can be beneficial if it facilitates handing of data regarding important aspects of a given phenomenon. As is typical of experimental science, it is not necessary to initially have 100% of the picture to make relevant cognitive progress. This general consideration also applies to the performance of our progressively developing dynamic model.

Extensive pretraining of our model, inputting information from >100 pertinent articles, is the predominant strategy we have followed. We achieved semi‐quantitative predictivity (statistically very significant), and noticed that our parameters can be discretized in eight log‐scale intervals (ranges) keeping the same good correlation with experimental results. A discretization in four log‐scale intervals (ranges) is suboptimal, but still correlated (unpublished data), thus our model is performing quite well.

We aim to further improve our model with more extensive acquisition of parameters (proteins, P‐protein concentrations and improved reaction rate data). Recent observations suggest that variations in protein concentration, in different tumor types, could explain (at least in part) different responses to altered pathways inhibitors.104

A direct parameterization of rates requires known or derived structures of the involved molecules. Sophisticated computational pipelines are available to predict protein structure, and protein–protein interaction sites. We need extensive docking simulations to derive all possible interactions among the molecules, and can use molecular dynamics simulations to derive Gibbs free energies and to predict equilibrium concentrations. Introducing a general approximation at the level of association rates, we can finally derive from the Gibbs free energies both association and dissociation rates.105

Starting from an initial ‘physiological condition,’ the model can be adapted to simulate individual cancer pathologies, implementing (sequentially) discrete numbers of alterations/mutations in relevant onco‐proteins. We verified some salient model predictions using the mutated CRC lines HCT116 and HT29. We also validated the behavior of our pretrained model against subsequent external results.

During the training phase of our software, we also experimentally revealed that MEK inhibition by a MEK inhibitor drug (CI1040) decreases CCND1 mRNA stability via the transcriptional regulation of an intermediary gene.90

Mutations/alterations present in the HCT‐116 line are: ErbB2 1X, PTEN (60%), KRAS, PI3K, β‐Catenin, TGFβ receptor II, and E‐Cadherin.

Mutations/alterations present in the HT‐29 line are: ErbB2 2X, PTEN (100%), BRAF, PI3K, APC, and SMAD4.

The following oncoprotein inhibitors were used: Azakenpaullone (GSK3β inhibitor), CI‐1040 (MEK1/2 inhibitor), Perifosine (AKT inhibitor), PI103 (PI3K inhibitor), and XAV939 (promotes β‐catenin degradation).

The syntactic rules for drawing a MIM were according to the papers90, 92, 93, 94, 95, 96 and Table 1.

Table 1

Summary of our Biochemical Interactions/Reactions Pathways (Fuzzy Logic Definitions)

1. Pathway [ErbB‐family receptors – PI3K – PTEN – AKT – ramifications a, b, c, d];

2. Pathway [ErbB‐family receptors – Shc – Grb2 – SOS– GAP– KRAS – BRAF – MEK – ERK – AP1 – TFBSAP1, transcription agonist];

3. Pathway [ErbB‐family receptors – E‐Cadherin (Cadherin/Catenin adhesive complex)];

4. Pathway [ErbB‐family receptors – PLCγ – PIP2 – PKC – BRAF – MEK – ERK – AP1 – TFBSAP1, transcription agonist]; the terminal parts of pathway 2 and 4 are the same];

5. Pathway [WNT – Frz/LRP5/6 – Dvl – AXIN – APC – GSK3β – β‐catenin – TCF7L2 – TFBSTCF7L2, transcription agonist];

6. Pathway [TGFβ‐receptors – SMAD2/3 – SMAD4 – TFBSSMAD, transcription antagonist];

7. Pathway [TGFβ‐receptors – TAK‐1 – TAB2 – NLK – TCF7L2 – TFBSTCF7L2 (TCF7L2 binding site), transcription agonist], converging with 8];

8. Pathway [WNT – Frz/LRP5/6 – TAK‐1 – TAB2 – NLK – TCF7L2 – TFBSTCF7L2, transcription agonist], converging with 7].

9. Pathway [Integrins – ILK – PIP3 – AKT – ramifications a, b, c, d]

10. Pathway [Integrins – FAK – AKT – ramifications a, b, c, d]

11. Pathway [Integrins – FAK – PI3K – AKT – ramifications a, b, c, d]

12. Pathway [Integrins – FAK – Grb2 –SOS– GAP– KRAS – BRAF – MEK – ERK – AP1 – TFBSAP1, transcription agonist]

13. Pathway [Integrins – FAK – SRC – ramifications e and f]

14. Pathway [ALK – SRC – ramifications e and f]

15. Pathway [ALK – PI3K – ramifications a, b, c, d]

16. Pathway [ALK – PLCγ – PIP2 – PKC – BRAF – MEK – ERK – AP1 – TFBSAP1, transcription agonist]

17. Pathway [c‐Met – Shc– Grb2 – SOS– GAP– KRAS – BRAF – MEK – ERK – AP1 – TFBSAP1, transcription agonist];

18. Pathway [c‐Met – PI3K – AKT – ramifications a, b, c, d];

19. Pathway [c‐Met – SRC – ramifications e and f];

AKT Ramifications

a. Pathway [AKT – GSK3β – APC – β‐catenin – TCF7L2 – TFBSTCF7L2, transcription agonist]

b. Pathway [AKT – mTOR – p70S6K]

c. Pathway [AKT – MDM2 – TP53 – TFBSTP53, transcription antagonist]

d. Pathway [AKT – P21 – Cyclin (D/E) / CDK (2/4) – pRB – E2F:DP – TFBSE2F:DP, transcription agonist]

SRC Ramification

e. Pathway [SRC – Grb2 – SOS– GAP– KRAS – BRAF – MEK – ERK‐ AP1 – TFBSAP1, transcription agonist]

f. Pathway [SRC – PI3K – AKT – ramifications a, b, c, d]

How We Modeled mRNA Levels

We implemented a thermo‐statistical derivation of a Transcription Rate Function for MYC and CCND1 (applicable to other genes transcribed at the G1–S boundary). The complex TCF7L2:β‐catenin cooperates (both through positive and negative modulations) with other transcription factors illustrated in Figure 3.

Regulation of transcription is typically multifactorial, involving a series of transcriptional activators, repressors and co‐factors that control the recruitment of the transcriptional machinery and RNA Polymerase (RNAP) to the transcription start site.

We applied a statistical thermodynamic framework to relate MYC and CCND1 transcription rates to the concentrations of their upstream transcriptional activator and repressor complexes. Notice that MYC and CCND1 are just two mRNAs, selected out of many parallel concomitant examples, for instance: CCND3, JUN, MYCN, CCNE1, CCNE2, CDC25A, CDK2, E2F1‐3, NPAT, MYB, MYBL2, and TFDP1. Among the negative cell cycle regulators are: CDKN1C, CDKN2C, CDKN2D, E2F4‐8, RB1, RBL1, and TP53.

The first step in building the thermo‐statistical model involved the identification of key Transcription Factor Binding Sites (TFBSs‐promoters) responsible for activation and repression of MYC, CCND1 and other genes, as well as the main TFs that bind to them.

Details of the computations are given in our 2015 paper90 and its supplementary material.

In 2015, we verified some salient model predictions using the mutated CRC lines HCT116 and HT29, and measured the amount of MYC mRNA, CCND1 mRNA, and AKT and ERK phosphorylated proteins, in response to treatments with different onco‐protein inhibitors, alone or in combination90 (Figure 4).

Figure 4

Response to different inhibitor treatments. Scatter plots of experimental (Y‐axis) versus simulated (X‐axis) values. Notice that only the Y‐axis represents experimental variability. Top graph to the left: ERKPP protein levels. Top graph to the right: AKTP protein levels. Bottom graph to the left: MYC mRNA levels. Bottom graph to the right: CCND1 mRNA levels. The (0;0) origin of the two axes makes reference to a theoretical complete inhibition, both experimental and simulated. Notice the (0;0) origin is never reached, either in experimental or simulated values, even for strong inhibitions. Similarly, the 45° dotted diagonal makes reference to an ideal complete correspondence between simulated and experimental data. We used the following inhibitors and inhibitor combinations—1, 2: Controls; 3: XAV939; 4: PI103; 5: CI1040; 6: Perifosine 20nM; 7: Perifosine 40nM; 8: XAV939 + PI103; 9: XAV939 + CI1040; 10: PI103 + CI1040; 11: Perifosine 20nM + CI1040; 12: Perifosine 40nM + CI1040; 13: XAV939 + PI103 + CI1040; 14: XAV939 + PI103 + CI1040 + Perifosine 20nM; 15: XAV939 + PI103 + CI1040 + Perifosine 40nM. Figure elaborated from results presented in Ref 90.

The correlation between simulated and experimental values was statistically significant in all cases considered, at a P‐value of at least <0.01 (two‐tailed). Statistical tests we implemented include Spearman's rho and R 2 coefficient for linear regression. Details are reported in Ref 90.

Our MIM reconstruction, parameterization, and mathematical modeling, is not a final working instrument. However, a continuously updatable working instrument appears realistically achievable in the near future. What we have already shown is that some deficiencies in input (inevitable at this stage), are definitely not obscuring the connection between model and experiments.

The very important message (reinforced by statistical analysis) is that an incremental path is now open, justifying more long‐term future cooperative efforts in the direction of building larger MIMs, supported by more input parameters and the direct inclusion of molecular concentrations and reaction rates.

Our present model is already quite encouraging, and we continue to work toward future developmental advancements in the model. Our long‐term goal is its utilization for the proposal of reasoned combinations of inhibitors of specifically altered pathways in the specific cancers of individual patients.

Notice that in terms of differential behavior between a normal network and a network carrying altered/mutated pathways, pathways that remained normal in a cancerous state tend not to contribute to the differential behavior between a normal and a cancer counterpart. These considerations apply in principle also to pathways not yet inserted into our MIM.

Looking at our extensive pretraining approach, some of the published literature could contain situations close to our experimental validations. Implicit/partial conditions of retrofitting are difficult to evaluate when your pretraining utilizes the pertinent fraction of data reported in results of >100 papers. However, a partial overlap between pretraining and validation could suggest that we are already covering the most relevant features of our G0→G1 network subregion.

Our Model Applied to Independently Published Results

Well after our model finalization and our own experimental verifications, the predictions of our model were tested against preclinical results obtained by independent investigators.106

The authors had examined the DiFi, LIM1215, HCA‐46, and OXCO‐2 CRC lines, before and after induction of panErb resistance, through a subsequent KRAS mutation. These lines were initially sensitive to the panErb inhibitors cetuximab or panitumab (both monoclonal anti EGFR antibodies), but resistance emerged through subsequent new KRAS mutations. The authors observed that, in their CRC lines (initially sensitive and then resistant to panErb inhibitors) the addition of MEK inhibitors (downstream in the same KRAS pathway) could only partially overcome resistance. The initial (pre KRAS mutation) sensitivity to the panErb inhibitors suggests an activation of probably more than one pathway starting with the tetramer [(ligand:EGFR)2]. Not only had the MAP‐kinases pathway to be involved, the panErb inhibitor would also continue to inhibit an additional pathway along PI3K→AKT. Following the KRAS mutation, administering a MEK inhibitor in combination with a panErb inhibitor led to more complete resensitization to, and thus effectiveness of, the therapy.

With our model, we examined the simulated behavior before and after the emergence of resistance to panErb inhibitors. To do so we generated ad hoc MIM modeling, in the absence or presence of KRAS alterations. We simulated the presence of a panErb inhibitor, a MEK inhibitor, or both.

The behavior of P‐protein/total protein for EGFR, ERK, and AKT (at 30 min–1 h) was compatible with the authors’ observations. Moreover, both c‐MYC and CCND1 mRNAs (both playing a crucial role for cell replication), were completely normalized (at 4–8 h), only by the combination of panErb and MEK inhibitors in the presence of a mutated KRAS.

Our simulations suggest that this behavior is due to a synergic effect of the two inhibitors, which target two different pathways:

MEK inhibitor downstream of the acquired KRAS mutated pathway.

PanErb inhibitor on the PI3K → AKT pathway.

This synergism appears in line with the one observed by the authors at a cellular level. The behavior of our model in this new context was quite encouraging, and able to reflect satisfactorily these independent experimental findings.

An epithelial organoid culture system, in which human intestinal stem cells (ISCs) indefinitely self‐renew and form crypt‐like organoid structures in Matrigel, has recently been developed.107, 108 Using CRISPR‐Cas9–Cancer lines mediated engineering of human intestinal organoids, they can be progressively transformed into neoplastic lesions.109, 110

Human colorectal tumors bear recurrent mutations in genes encoding proteins belonging to the WNT (APC), MAPK (KRAS), TGF‐β (SMAD4), TP53, and PI3K pathways, all present in our MIM.

We introduced in our dynamic simulator all possible permutations (32) for 0 (1), 1 (5), 2 (10), 3 (10), 4 (5), and 5 (1) mutations, and present the results obtained in Figure 5.

Figure 5

We introduced, in our dynamic model, (in the abscissae) all the possible permutations of the five mutations (and mutated pathways) present in the organoids of109, 110: WNT (APC), MAPK (KRAS), TGF‐β (SMAD4), PI3K, and TP53 pathways. In the ordinates, we show the increased transcription rates of the complex described in Figure 3. At 4–5 mutations and 20–30 A.U. of mRNAs levels, we appear to cross the cell‐cycle restriction point. At position −1 in the abscissae, we show the strong inhibitory effect of a deleted β‐catenin.111

As shown in Figure 5, the results of our model were in very good agreement with the experimental findings cited in the following references.109, 110 We can consider the transcription complex of Figure 3 (and related equations) as acting in parallel to the stimulus to cross the boundary of the restriction point between G1 and S cell cycle phases. We can consider an mRNA ratio 20–30/1 as roughly a transition threshold for a multimutated organoid, versus an mRNA ratio = 1 for 0 mutations (normal organoid). Notice that mice carrying a KO β‐catenin are unable to develop normal crypts.111

As shown by Figure 6, the behavior of our model is again quite reasonable, even more so if we analyze the performance of individual combinations: inhibitors of more distant /independent pathways tend to synergize more strongly.

Figure 6

We started with a dynamic model of,109, 110 carrying mutations in: WNT (APC), MAPK (KRAS), TGF‐β (SMAD4), and PI3K pathways. We gave simulated combinations of the following inhibitors: PanErb inhibitor, MEK inhibitor, PI3K inhibitor, AKT inhibitor, and MDM2 inhibitor. We show the effects of the simulated combination of the five inhibitors. We show permutations of 0 (1), 1 (5), 2 (10), 3 (10), 4 (5), and 5 (1) inhibitors.

CONSIDERATIONS AND QUESTIONS FOR FUTURE PROSPECTS AND DEVELOPMENTS

We think we are on a viable research track, not a presently impassable route, even considering the availability of only partial information. Our dynamic modeling could offer informatics based support for decisions about reasoned combinations of targeted therapies. Kinetic parameters of inhibitors are known and they can be introduced in a straightforward way in our kinetic model.

In the COSMIC release of June 2, 2014, the curators have estimated 4–10 driver mutations (or somatically inheritable alterations) as possible causes of an individual cancer. They are not the same, or in the same number, for different tumors. They have evolved against a background of more than 10,000 passenger mutations per tumor.79

At the stage of reconstruction of our MIM of March 2015,90 the most frequent driver and gate‐keeper mutations/alterations were already represented in the MIM, TP53, APC, KRAS, PTEN, SMAD4, PIK3CA, BRAF, and CDH1 (in an order of decreasing frequency of occurrence).

According to the COSMIC release of September 8, 2015, focused on colon carcinoma,79 we have the following probability of presence of a given mutation in an individual cancer (20 most frequent oncogenes): TP53 (48%), APC (42%), KRAS (35%), ATM (23%), PIK3CA (22%), PTEN (21%), SMAD4 (20%), FBXW7 (19%), PTCH1 (16%), ARID1A (15%), CREBBP (15%), KMT2D (15%), BRAF (14%), NF1 (14%), RB1 (14%), KMT2C (13%), KIT (12%), TRRAP (12%), CARD11 (12%), and BRCA2 (11%).

In the list shown above, we have both gatekeeper genes directly related to the signaling‐network, and caretaker genes with diversified complex functions (both fidelity/accuracy of replication and transcription). We will not enter an exhaustive analysis, but a first look suggests that a clear majority of gatekeeper genes are also present in our MIM and in our dynamic model.

From the perspective of the individual tumor of a specific patient, we could have had a Darwinian evolution to cancer, through a constellation of some frequent driver mutations/alterations and much less frequent driver mutations/alterations.112 This consideration could also apply to intratumor subclones.113

In a modern framework of personalized Oncology, to know more about these individual oncogene constellations is becoming increasingly relevant. An uncommon genotype/epigenotype in a patient could confer varying sensitivity or resistance to a specific inhibitor, and/or (more importantly) a specific combination of different inhibitors.

Several Cancer Data Portals28, 79, 114, 115, 116 can allow for in depth mutational analysis of individual tumors. In principle, for each tumor, both very frequent driver‐gatekeeper mutations and much less frequent mutations can be observed and considered.1 Individual cancers could depend on a combination of commonly altered pathways and more ‘private/unique,’ altered pathways (uncommon, but still crucial, for an individual cancer).

We aim to identify and characterize somatically inheritable driver and gatekeeper alterations, affecting cancer relevant signaling‐proteins in individual patient tumors, a task requiring advanced technology and extensive acquisition of information. Inputting this information into, and processing it with, our model, we hope to identify optimal treatments to improve patient outcome.

The multihit evolutionary process of malignant transformation is potentially causing some intratumor heterogeneity. To complicate matters, in addition to common mutations frequently presenting in the population, some subclone‐unique mutations may be present. Cancer‐genes’ mutations detected from liquid biopsies could be informative (with the benefit of easy collection), as could multiple biopsies on the primary tumor and metastases at different time points.113, 117

A very interesting finding was that CRCs with acquired cetuximab resistance (via KRAS mutation) often demonstrated a decay of KRAS mutant clones upon antibody withdrawal, thus conferring renewed cetuximab sensitivity! This indicates that intermittent drug schedules could be highly beneficial,29 and our modeling could be used to identify and exploit this and other tumor weaknesses. Since our MIM can simulate characterized individual clones one by one, in the future we may consider intratumor subclones in individual patients. Exploiting circulating tumor DNA (ctDNA), to genotype CRCs at different times of cancer evolution or therapy, appears potentially promising. The complexities of late clonal evolution could thus be exploited using repeated liquid biopsies, at different stages of patient care, to determine optimal times to start or interrupt different drugs. Insertion of these dynamic findings into our model can make it even more realistic and useful.

Perhaps, in a given cancer, we should investigate at least 100–200 of the most frequent driver and gate‐keeper alterations, to detect a clear majority fraction of them, as well as lower (2–3%) frequency mutations.

Driver‐gate‐keeper alterations involved in the G0→G1 transition are especially relevant because they lead into to the restriction point just prior to the S phase, however these should be distinguished from those involved in DNA replication fidelity but not actually required in the signaling‐network of biochemical interactions involved in cell replication control.

In alternate to mathematical modeling of the type illustrated above, multiple parameters are conventionally considered in an unconnected way, substantially as a richer set of classic bio‐markers (a nonsystems approach). If the set of biomarkers is relevant, we will be able to show important statistical correlations with prognosis and response to treatment (the two are partially correlated).

We can think in terms of well conducted and directed –omics investigations, at a variety of levels and combined levels (various DNA alterations, mRNAs, mi‐RNAs, nontranslated RNA genes, proteins, and metabolites of any kind).

We need sufficiently powerful and expensive analytical lab tools, correspondingly strong computing hardware and software, and to acquire a new larger set (and different classes) of biomarkers.

For the moment, systematic ‐omic approaches would be too expensive for routine clinical use. A successful large ‐omics pilot study, could identify a much smaller significant set for analysis (again a biomarker, or nonsystem nonnetwork‐integrated, approach) thus reducing clinical costs markedly while improving patient outcomes (see also Proposed Panels of mRNAs section).

A side comment: We can ask the side, not irrelevant, question, if we are developing a cancer care system available for only some technologically developed country and sufficiently wealthy people.

Toward Future Strategies

A translational application of our dynamic model approach envisages three main and interconnected steps.

In order to diagnose driver mutations/alterations in each patient, we can choose to detect, by NGS technologies, a discrete number of genes/signaling‐proteins, involved in an important phase of the cell cycle (for instance the G0→G1→S transition, as discussed above). Perhaps, in a given tumor (CRC in our case), we should investigate at least 100–200 of the most frequent driver and gate‐keeper alterations, to detect a clear majority fraction of them, including mutations present at a relatively low frequency.

High throughput technologies, intended not only for detecting mutations, but also copy number alterations and somatically inheritable epigenetic changes, allow for in depth analysis of alterations in individual tumors. These new technologies have contributed to highlight how individual cancers could depend on a combination of commonly altered pathways + ‘private/unique’ altered pathways (uncommon, but still crucial, for an individual cancer).

The second step requires the integration of mutation/alterations of each patient within connected molecular‐interaction pathways (a MIM) belonging to a network subregion, along which biochemical signals are propagated (see A Dynamic Signaling‐Network Model and How We Modeled mRNA Levels sections).

Starting from a MIM, we can derive a decision making strategy from the outcomes of a dynamic model carrying the gatekeeper mutations/alterations of an individual cancer, plus available targeted inhibitors of detected altered signaling proteins affected by excess of function, with the goal of proposing rational combination therapies. Using the web and drug‐focused databases,118, 119 we can suggest distinct independent targeted inhibitors, in principle one inhibitor for each altered pathway. New molecules available for clinical trials will be coming into play within the next 3–5–10 years, and more are in the pipeline. They are the research focus of most Pharmaceutical Companies, large, medium, and small, often in synergism with the academic world.

The input and rational suggestions coming from our dynamic model, combined with drug data, will have as a consequence that the platform of required pharmacological validations, both at molecular and cellular levels (inhibition of cell growth, stimulation of cell death, etc.), will be very restricted, and therefore more realistically feasible.

This final step represents an alternative to present time consuming procedures (necessitating long patient recruitment periods): small patients subsets sharing a multiplicity of common biomarker features, and extracted from larger, more heterogeneous sets. Intratumor heterogeneities could potentially generate even smaller patient subsets.113

Discussion of the Future Possibility, and Need for, of a New Kind of Clinical Trial

Large, poorly homogeneous trials are no longer sufficient. Thanks to the statistical power of large numbers of patients, even relatively small differences generally tend to show significance, perhaps a misleading result.

After 50 years of conceptually similar studies, and similar gains, if each study compared a therapy to the best treatment option previously available, and we gained only three‐months of survival time every year, most metastatic solid carcinoma tumors should have become chronic diseases, rather than the deadly illnesses which they often are. Current approaches are not enough! There is an urgent need for change!

To give an idea of typical modern clinical trials (and the difficulties they encounter), we can briefly introduce the reader to ‘umbrella trials’ and ‘basket trials.’

The umbrella design focuses on a single tumor type or histology.120 In contrast to the umbrella design, basket trials allow the study of multiple molecular subpopulations of different tumor or histological types, but sharing the same molecular pathology, all within one study.

The NCI‐MATCH study ‘Molecular Analysis for Therapy Choice’ is a new basket trial aimed at exploiting shared mutations, encountered across many cancer types, by matching them to targeted drugs. The drugs studied are associated with predictive biomarkers, and include both U.S. Food and Drug Administration–approved and investigational drugs.121 The trial (which started enrollment in July 2015) will screen up to 3000 patients with refractory solid tumors or lymphoma, and aims to enroll 35 patients into each of 20–25 biomarker subgroups. This novel type of study attempts to match patients sharing driver genetic abnormalities (regardless of tumor histology) with drug(s) expected to work on mutated pathway(s). Shared mutations are the elements of similarity across the arms in a basket trial, and although each patient subset includes a small numbers of patients, it is larger than feasible with the ‘umbrella trials’ added requirement of tissue communality.120

The identification of sufficient numbers of patients carrying the same individual (rare) combination of driver genetic aberrations implies long recruitment times, especially if we want a homogeneous tissue of origin,122 so if we can remove this requirement it could speed discovery.

A Tentative Suggestion for the Future: Meta‐Analysis of Multiple «One‐Patient Trials», through an ‘Invariant’ Mathematical Dynamic Modeling Procedure

We can consider our proposed approach as a different type of basket trial: the tissue is the same (like in an umbrella trial), but what we share is the mathematic modeling procedure used to suggest different ad hoc combination therapies, for tumors carrying different combinations of driver gatekeeper mutations.

In our approach, the procedure would be invariant because it would be using the same mathematical dynamic model, applied to the same solid tumor (CRC in our case), and because it is focused only on the quite manageable G0→G1→S cell cycle transition (see above). It appears logically possible to assess the efficacy of combination therapies suggested by the same mathematical dynamic model, through a meta‐analytic approach joining the equivalent of the outcomes of many ‘one patient trials.’123

These ‘one patient trials’ will have to have followed the suggestions coming from the same mathematical dynamic model, with the only variables being: the molecular pathology (mutations, somatically inheritable alterations) of each individual patient and the inhibitors suggested by the altered pathways.

Having used the same mathematical dynamic model (for suggesting a precise combination therapy) is the logical equivalent of a deeply homogeneous patient meta‐subset.

We think of a meta‐set of ‘1 patient sets.’ We think that a new innovative strategy could come from more advanced versions of the mathematical modeling we have been working on for some years.90, 124, 125, 126, 127

The ‘homogeneity of a treatment strategy’ would come from the ‘homogeneous procedure’ adopted for selecting the personalized‐treatment‐combination, downstream of efficient network reconstruction and modeling,90, 127 rather than traditional methods (of genetic/epigenetic biomarker‐homogeneity within a small subset of patients) which require long recruitment times.

A well performed meta‐analysis of a set of ‘1 patient subsets’ would in principle take care of the problem. We aim to perform a meta‐analysis, based on the assessment of Overall Survival (OS) and related clinical parameters (see RECIST 1.1 guidelines, for details128), of patients treated according to the pharmacological suggestions coming from dynamic modeling. We skip here additional clinical details, including the possible side effects/toxicities of combination therapies (which would be addressed at a later stage). Notice that umbrella trials, sharing one gatekeeper mutation and one corresponding inhibitor, frequently do not provide breakthroughs.129 It is difficult to anticipate what could be the average behavior of a meta‐analysis on the wavelength that we suggest here, but we project that it could be used in the next 3–5 years.

If we assume the potential validity of dynamic modeling, our approach, especially considering the continuous appearance of new drugs, could save precious time and lives. A phase‐1 toxicity study should obviously precede or have preceded these ‘1‐patient trials,’ according to usual dose escalation strategies.130

With patience and perseverance, we should try to make ethically and legally acceptable ‘treatment sets’ homogeneous for the procedural decision strategy adopted (consequence of the homogeneity of our dynamic modeling software), rather than for the personal genetic analysis of each individual patient per se (as an element of the set of a usual trial).

Considering the inevitably slow pace of standard patient‐recruiting procedures, it is important to devise a more efficient strategy to achieve rational drug combinations for each personal cancer case. Applying a standard approach (according to a nonsystem bio‐marker‐set mentality) entails recruiting difficulties in building homogeneous subsets, both in terms of subset size and required time for patient recruiting.

Assessing the performance of our new proposed approach could tentatively be feasible, but would require a partially revised set of ethical trial rules. We would give origin to a new type of basket trial based on a unique dynamic model combined with the genetic landscapes of individual tumors and drugs congruous with the altered pathways. In a broader perspective, in his recent book Dr. Vincent DeVita has repeatedly advocated for revised and more patient inclusive trial rules.76

In conclusion, it is worthwhile to explore the future possibility of a new type of proposed trial, to facilitate/accelerate the achievement of a rational combination of targeted drugs, tailored to the cancer of each individual patient. Our dynamic modeling could be conducive to informatics based decision support for reasoned combinations of targeted therapies.

CONCLUSIONS

We have briefly illustrated the process of carcinogenesis in CRC as a multihit event involving driver‐gatekeeper mutations and/or somatically inheritable alterations. This basic feature of the malignant transformation process is a characteristic shared in all cancers.

In this review, we primarily focus on and aim towards combination therapies with targeted drugs: inhibitors of signaling proteins affected by excess of function.

We start with the practical level of clinical practice guidelines, which are highly relevant as following them is associated with a significant survival benefit (as seen in patients following NCCN guidelines in the United States).70

Characterizing CRC tumors at the biomarker level can be relevant diagnostically, and for guiding clinicians to more appropriate and selective therapy.

Crucial biomarkers are driver mutations, somatically inheritable driver alterations/epimutations. Within the actual alteration of a signaling‐network, gatekeeper alterations are directly responsible for the cancer pathology of the network. Caretaker alterations are important in the acceleration of the carcinogenetic process; they reduce DNA replication fidelity, thus accelerating the appearance of new gatekeeper mutations.

From the central stage of gatekeeper mutations/alterations, we move to proposed sets/panels of mRNAs (and others RNAs). When satisfactorily validated and sufficiently sensitive and specific, they can be used to inform not only diagnostic and prognostic decisions, but also therapeutic choices.

MicroRNAs and long noncoding RNAs are potentially relevant in the process of malignant transformation. They are presently investigated at a more basic‐research level, but some of them could become part of clinical practice (both for diagnostic and therapeutic applications) within a few years.

The focus of the central part of our review considers how to propose reasoned combinations of targeted inhibitors.

The proposed strategy is the implementation of a dynamic signaling‐network model. We illustrate an example of this model that appears to behave quite promisingly in four independent validation instances.

A fascinating application of this model is the possibility of implementing a homogeneous meta‐analysis of multiple ‘one patient trials,’ through a mathematically invariant utilization of the model, where only driver gatekeeper mutations/alterations concerning a given individual patient and targeted drugs suggested by the altered pathways, are the variants introduced in the model. In the future: a clinical trial, meta‐analytically analyzing one‐patient trials, could assess the statistical validity of proposed drug combinations much more rapidly than traditional small‐targeted trials. This will produce some nonbreakthrough combinations,129 but has the potential to accelerate the discovery of the way to save individual lives.