Longitudinal Liquid Biopsy and Mathematical Modeling of Clonal Evolution Forecast Time to Treatment Failure in the PROSPECT-C Phase II Colorectal Cancer Clinical Trial

Colorectal cancer is one of the most common cancers worldwide and accounts for 8% to 9% of cancer-related mortality (1), with poor 5-year survival for stage 4 disease (2, 3). Genetic alterations in the RAS pathway are responsible for primary and acquired resistance to anti-EGFR monoclonal antibodies (4–12). However, despite tailored patient selection based on genetic screening for somatic RAS mutations, 65% to 70% of patients progress within 3 to 12 months after starting therapy.

In many patients, RAS mutations occur early during colorectal carcinogenesis, manifesting as a clonal (truncal) alteration in primary and metastatic lesions (13, 14); these patients derive no benefit from anti-EGFR therapies. In seminal studies, plasma cell-free DNA (cfDNA) analysis of patients with no detectable RAS mutations at baseline has demonstrated that under the selective pressure of therapy, small and undetectable RAS-mutant subpopulations at baseline undergo clonal expansion, ultimately leading to acquired treatment resistance (15–17). Hence, although resistance to anti-EGFR therapies can be polyclonal, it often converges on biochemical signaling pathways downstream of EGFR (18). Because treatment resistance is driven by intratumor heterogeneity (ITH) that generates the substrate of variation necessary for evolution, multiple biopsies in time and space are critical to understanding the complexity of an evolving malignancy (19, 20).

Although cfDNA holds immense promise for patient management, the predictive power of liquid biopsies in metastatic colorectal cancer (mCRC) has not been demonstrated within a prospective clinical trial, and data on concordance between liquid (cfDNA) and solid (tissue) biopsies remain sparse at this stage. Most of the previously reported cohorts of anti-EGFR–resistant mCRC include patients treated with anti-EGFR therapy in combination with chemotherapy, representing a potential confounding factor in understanding genomic patterns. Moreover, although seminal studies in other cancer types such as breast cancer (21) and non–small cell lung cancer (NSCLC; ref. 22) have shown prospectively that cfDNA positivity precedes clinical progression, a major challenge to the use of liquid biopsies in the clinic is the extreme variability among patients, with times to progression ranging from 1 to 400 days in breast cancer, 10 to 346 days in NSCLC, and unknown in mCRC. This variability prevents patient-specific clinical predictions, representing an obstacle for personalized medicine.

We have recently shown that combining functional genomics, molecular pathology, and radiology in the context of well-annotated clinical trials can help in defining novel biomarkers to optimize patient selection (23, 24). In this prospective phase II trial of patients with RAS wild-type (WT) mCRC treated with single-agent anti-EGFR monoclonal antibodies, we aimed to (i) test the value of profiling subclonal mutations in the RAS pathway in cfDNA to predict response to anti-EGFR therapy; (ii) assess the mutational concordance between liquid and solid biopsies; and (iii) estimate time to treatment failure in each individual patient using mathematical modeling of cancer evolution.

The PROSPECT-C trial (ClinicalTrials.gov identifier: NCT02994888) is a study of biomarkers of response and resistance to anti-EGFR therapies in KRAS/NRAS WT chemorefractory mCRC. The trial recruited 47 patients between November 2012 and December 2016. The objectives of the study were to track and validate the known mechanisms of resistance/response to anti-EGFR therapies and to identify novel mechanisms of response/resistance to such therapies.

Patients in the study were subjected to tissue sampling from metastatic deposits at predefined time points including pretreatment [baseline (BL)] and posttreatment [disease progression (PD)]; optional partial response (PR) biopsies were conducted where clinically and technically feasible, and archival material (primary cancer or original diagnostic biopsies) was used when available. Plasma was collected every 4 weeks until disease progression (Table 1; Fig. 1A; Supplementary Table S1). Median duration of cetuximab treatment was 10.7 months (interquartile range, 2.0–137.3); 20.0%, 24.4%, and 46.7% of the patients experienced PR, stable disease, and PD (indicative of primary resistance), respectively, by RECIST v1.1. Median progression-free survival (PFS) and overall survival (OS) were 2.6 months (95% CI, 1.9–4.2) and 8.2 months (95% CI, 4.2–12.0), respectively (Supplementary Table S2). These results were found to be consistent with the previously published literature (25).

A consort diagram describing the study is reported in Fig. 1B. Patients were tested for RAS pathway aberrations in cfDNA and tissue biopsies according to different methods as described below.

In the initial 22 consecutive eligible patients (cohort 1), cfDNA was tested by digital droplet PCR (ddPCR) using a tiered approach based on frequency of RAS pathway aberrations previously reported as being associated with primary and acquired resistance to anti-EGFR treatment (18). Patients in this cohort also had tissue biopsies sequenced: A total of 84 cores from 22 of 22 (44 cores), 4 of 4 (8 cores), and 16 of 22 (32 cores) biopsies at BL, PR, and PD, respectively, were tested; 6 patients did not have progression biopsies (2 deceased prior to assessment, 2 were considered clinically unfit, and 2 declined biopsy). Twenty available archival samples from 8 patients were also sequenced (Supplementary Table S1; Supplementary Fig. S1).

From the remaining 23 eligible patients on the trial, PD and/or BL plasma samples from patients (n = 17, cohort 2) with primary resistance [PFS ≤ 3 months (11 patients, 12 samples)] and long-term benefit [PFS > 6 months (4 patients, 13 samples)] were subjected to targeted sequencing with the Roche Avenio cfDNA Expanded Kit, covering 77 cancer-related genes. In the 2 long-term responders for whom material was available from the first cohort (4 samples), a similar approach was adapted to validate our findings by sequencing BL and PD cfDNA (Supplementary Fig. S1).

We initially investigated cfDNA concentration (Supplementary Table S3) and RAS pathway hotspot mutations in 143 plasma samples using ddPCR (10 samples could not be tested because of hemolysis; Supplementary Fig. S1). Eleven patients (50%) had RAS pathway aberrations in their baseline cfDNA: 6 patients with KRAS/NRAS mutations, 2 patients with BRAF^V600E mutations, 1 patient with PIK3CA^E545K mutation, and 2 with ERBB2 amplification (Supplementary Fig. S2; Supplementary Table S4). Detection of RAS pathway aberrations in baseline cfDNA was significantly associated with inferior PFS (HR, 3.41; CI, 1.24–9.37; P = 0.02) and worse OS (HR, 2.78; CI, 1.09–7.11; P = 0.03), and showed a trend toward poor response rate (0% vs. 36.4%; P = 0.09) compared with WT patients (Fig. 2A–D). In order to validate these findings, we tested the prevalence of RAS pathway aberrations in all the patients who showed primary resistance to cetuximab in the second cohort. Considering that hotspot-based methods such as ddPCR allow testing of only a limited number of known potential drivers, thus perhaps underestimating the biological and clinical relevance of polyclonal resistance, we subjected all baseline samples in the second cohort to next-generation sequencing (NGS) of a broad panel of 77 cancer-related genes. Interestingly, in keeping with data from the first set, known RAS pathway aberrations were found in 6 of 11 (54.5%) cases (Supplementary Fig. S3; Supplementary Table S5). In several cases, resistance could be attributed to multiple drivers (Fig. 2E). In the same cohort, 2 patients with no common RAS pathway aberrations displayed mutations in EGFR^T739P, EGFR^C326R, FGFR2^R203C, and KRAS^A18D. Among these, the KRASA18D is a gain-of-function mutation with transforming activity that lies within the GTP binding region of the KRAS protein (26) and might impair response to anti-EGFR treatment, thus deserving further validation (Supplementary Table S5). Interestingly, in either of the two cohorts, no significant differences in cfDNA concentration were observed between patients with or without baseline aberrations (Supplementary Fig. S4).

Next, we screened for drivers of acquired resistance in the first cohort of patients. Nineteen patients (86.3%) displayed RAS pathway aberrations at progression. As expected, the majority of patients with KRAS mutations had amino acid substitutions in codons 12, 13, and 61 (Supplementary Fig. S2). Interestingly, in 75% of the patients who achieved PR, aberrations were detectable in cfDNA in advance of clinical–radiologic progression: In the first patient, MET amplification was detected in cfDNA 2 months in advance of the radiologic progression. In the second patient, a KRAS^Q61H A-T mutation was detected 1 month prior to PD. In the last patient, who experienced a remarkable PR with PFS of 10 months, two mutations were detected. A KRASQ61H A-T mutation was found as early as 8 months ahead of PD, whereas another KRAS^G12D mutation emerged 4 months prior to radiologic progression (Fig. 2F). Polyclonal resistance due to multiple independent RAS-mutant subclones was observed in primary and secondary resistance cases and was supported by the slope of the APC clonal mutation that was used as a reference to infer subclonal composition (Fig. 2G; Supplementary Fig. S5). Specifically, polyclonal resistance could be determined because (i) the sum of the clone frequencies of distinct subclones cannot be larger than the whole, indicated by the frequency of a clonal mutation like APC (pigeonhole principle; ref. 27), and (ii) the trajectories over time of the two mutant subclones do not travel together and one can even overtake another, showing they are independent (Fig. 2G). The occurrence of polyclonal resistance in patients with acquired resistance (PFS > 6 months) was further confirmed by NGS of progression, best response, and baseline cfDNA of 4 patients with long-term clinical benefit in the second cohort and 2 patients in the first cohort for whom no mutations were initially detected by ddPCR (Fig. 2H; Supplementary Table S5). Multiple drivers of resistance emerged over treatment in 83.3% of cases and were confirmed to be subclonal by comparing their variant allele frequency to one of the truncal mutations in genes such as APC and or TP53 (Supplementary Fig. S5).

In keeping with previous evidence, RAS-mutant clones emerged during cetuximab treatment and faded away once treatment was interrupted (17). Intriguingly, despite a period of 12 months without any treatment and radiologic evidence of PD, these mutations remained undetectable when the patient entered fourth-line treatment, supporting the rationale for potential rechallenge with anti-EGFR treatment (Supplementary Fig. S6; ref. 28).

In 2 patients who showed a sustained benefit from cetuximab, lasting 33 and 16 months, respectively, we tested whether the durable response observed was due to persistent EGFR pathway inhibition due to a phenomenon of oncogenic addiction. Consistent with this hypothesis, an EGFR amplification was observed in the BL liver biopsy of 1 of these 2 patients (Fig. 3A). Treatment was halted in this patient after 16 months due to symptomatic progression of a nontarget lesion assessed by RECIST V.1.1 criteria. Interestingly, after cetuximab was withdrawn, rapid progression of the EGFR-amplified metastatic liver lesion that was biopsied at BL along with development of new liver deposits was observed within 6 weeks. In contrast, the soft-tissue pelvic mass biopsied at time of PD did not show any significant change in volume. Consistent with this, no EGFR amplification was found in the latter lesion, suggesting that this metastasis was not dependent on EGFR signaling (Fig. 3A). A rise in APC-mutant clones was observed synchronously with the increase in size of the nontarget metastasis that led to treatment discontinuation. Carcinoembryonic antigen (CEA) lagged behind, and no RAS pathway mutant clones were detected at any time point (Fig. 3B). Although these data confirm the correlation between tumor burden and cfDNA and highlight the higher sensitivity of cfDNA compared with CEA, a critical review of this case makes us wonder whether, in absence of circulating RAS pathway aberrations, this patient should have received radiotherapy to the single progressing metastatic deposit while continuing cetuximab treatment.

Based on the data gathered from our cfDNA analysis, we performed ddPCR validation and amplicon-based ultra-deep sequencing of sequential tissue biopsies (Supplementary Figs. S1 and S7) collected at clinically relevant time points to dissect the structure of RAS pathway aberrations and to test if the evolutionary patterns observed in cfDNA were represented in tissues. We started our analysis by comparing mutations detected in cfDNA with mutational data obtained by ddPCR in tissues. Twenty-four mutations were detected in cfDNA of 14 patients for whom at least one tissue sample was available for cfDNA/tissue comparison; 79% of cfDNA mutations were present in at least one tissue biopsy of all 14 patients. In patients with paired pretreatment biopsies and archival tissues (diagnostic material obtained prior to any treatment) available (6 patients and 9 mutations), 6 mutations from 5 patients were detected in both samples, suggesting that these mutations preexisted cetuximab treatment (Fig. 4A; Supplementary Table S6). In order to validate these findings, we ran amplicon-based ultra-deep sequencing of the same tissue biopsies using a custom library covering the most frequently mutated codons in KRAS, NRAS, and BRAF, plus the APC gene. NGS analysis showed agreement among different cores from the same lesion and was in line with ddPCR analysis, although many mutations could not be detected with statistical significance due to the limits of sensitivity of NGS versus ddPCR (Fig. 4A; Supplementary Table S7). Concordance in VAF between the two methods was good (R² = 0.5; P < 0.0001; Supplementary Fig. S8). Interestingly, ultra-deep sequencing showed the presence of a number of RAS pathway mutations in the samples analyzed (Supplementary Fig. S9). An average of 6.4 mutations (range, 3–9) was observed in different RAS hotspots with proven contribution to anti-EGFR treatment, supporting the notion that cetuximab resistance might often be polyclonal. In keeping with this observation, and in line with cfDNA data, purity adjustment of variants using clonal APC mutations confirmed that most of these mutations were subclonal (Supplementary Fig. S10). More importantly, the allele frequency of these abnormalities was below the detection threshold of clinically approved methods such as COBAS (Fig. 4B) or clinically validated targeted NGS at moderate depth (Supplementary Fig. S11), explaining why these patients were initially classified as RAS WT.

The ability to study multiregion and/or sequential tissue biopsies also allowed for the dissection of spatial and temporal heterogeneity in response to anti-EGFR treatments. Multiregion sequencing of two different cores of two liver metastases (segments II and IV) and primary cancer resected after neoadjuvant chemotherapy in one of our patients (patient 1009) showed the presence of 7 mutations in KRAS and NRAS hotspots. Among the mutations previously tested and detected in cfDNA of this patient, an NRAS^G12C mutation appeared as clonal and common to all the metastases, whereas the other 3 mutations (KRAS^Q61R, KRAS^G12S, and NRAS^G12D) appeared as subclonal, emerged during cetuximab treatment, and, as stated above, faded away once treatment was halted (Fig. 4C; Supplementary Fig. S6), suggesting that truncal and private mutations in the same pathway might coexist and contribute to cetuximab resistance. Two patients had ERBB2 amplifications detected in cfDNA; the amplification was confirmed by chromogenic in situ hybridization in BL and PD biopsies (Fig. 4D), supporting the concept that ERBB2 amplifications are clonal and present in approximately 10% of colorectal cancers (29). A MET amplification was observed in the cfDNA of a patient with PR 2 months in advance of radiologic progression and was supported by a synchronous increase in APC-mutant alleles (Fig. 4E, top). Fluorescent in situ hybridization revealed patchy areas of amplification in both the PR and PD biopsies of this patient (Fig. 4E, bottom). Furthermore, NGS analysis of his sequential biopsies including two different areas of the primary resected colorectal cancer along with two different cores at BL, PR, and PD revealed the emergence of an NRAS^G12C mutation, suggesting possible polyclonal resistance to cetuximab in this patient (Fig. 4F, left). Interestingly, the emergence of (MET and NRAS) resistant clones in the PD biopsy of this patient appeared to be associated with a decay in other subclones, consistent with the presence of an evolutionary bottleneck leading to the survival of the fittest, treatment resistant, clone (Fig. 4F, cartoon).

Our longitudinal data set offered a unique opportunity to study the dynamics of treatment resistance quantitatively because of the frequent (4 weekly) sampling.

Previous seminal studies showed that time-series treatment resistance data sets could be interpreted using mathematical modeling of tumor evolution (30). Similar types of modeling were also applied to longitudinal liquid biopsies to demonstrate that resistant subclones preexist anti-EGFR treatment (15). Although it is likely that resistant subclones are present at diagnosis, their extremely low prevalence in the treatment-naïve cancer cell population prevents their detection by standard assays such as COBAS, as exemplified in Fig. 4B.

Here, we sought to use mathematical modeling to jointly analyze both CEA (tumor burden) and cfDNA from individual patients, with the aim of exploring the predictive power of evolutionary principles when applied to a prospective clinical trial cohort. Finally, we validated our predictions using RECIST v1.1 measurements from radiologic imaging data.

In our model, at baseline a tumor consists of a total of |$N$| cells. These |$N$| cells are divided into two distinct subpopulations: a population of treatment-sensitive cells of size |${n_s}$| and a population of treatment-resistant cells of size |${n_r}$|⁠, with |${n_s} + {n_r}\; = \;N$|⁠. Sensitive cells die under treatment at rate |${\lambda _s}$|⁠, whereas resistant cells continue to grow under treatment at rate |${\lambda _r}$|⁠. This leads to the following equation for the change of cancer cell population over time during treatment:

Equation (1) models the exponential decay of sensitive cells and the exponential growth of resistant cells (Fig. 5A) and predicts a typical U-shape for the dynamic response to treatment (ref. 30; Fig. 5B). As CEA is proportional to the total tumor burden |$\;N( t )$|⁠, Equation (1) can be applied to CEA values over time from a patient. Figure 5C illustrates the U-shape dynamics for the cetuximab therapy schedule in patient 1014. Model fits to the CEA dynamics under cetuximab treatment showed remarkably high goodness of fit (average R² = 0.995) in patients with a sufficient number of measurements (≥3 time points, n = 32) and allowed estimation of the response rate |${\lambda _s}$| and the progression rate |${\lambda _r}$| for each case (see Methods for details and Supplementary Fig. S12). Response rates varied between 0 (nonresponders) and 0.58 per week with a median response rate of 0.2 (excluding nonresponders; Fig. 5D). Progression rates ranged from 0.03 to 0.38 per week with a median value of 0.165 (Fig. 5E; Supplementary Fig. S13). Additionally, the model allowed estimation of the initial frequency of the resistant population at treatment initiation. Responders had an initially small or slow-growing resistant subclone, or both. On the other hand, most nonresponders had a nearly dominant (∼100%) resistant population preexisting at baseline (Fig. 5F). Interestingly, in two cases (1002 and 1045), the resistant population at baseline was not dominant (8.3% and 4.1%, respectively), although not as low as the responders (<1%), but the growth rate of the resistant subclone was extremely high, leading to progression even before the first CT scan (Fig. 5F, gray dots, bottom right of the plot). This led to the patient being labeled a progressor, although our CEA analysis predicted that there was an initial response that remained clinically undetected.

Importantly, cfDNA can independently inform on the dynamics of the resistant population, modeled by the second part of Equation (1):

We applied Equations (1) and (2) to CEA and cfDNA-mutant frequencies, respectively, for those patients in our cohort for which enough time points were available for both measurements (≥3 time points, n = 11). We found that the model described the dynamics extremely well (CEA mean R² = 0.996, cfDNA mean R² = 0.979, Supplementary Fig. S14), highlighting responders in which rise of cfDNA preceded CEA (Fig. 5G) versus nonresponders in which both cfDNA and CEA rose at the same time (Fig. 5H). As discussed previously, in patient 1007, multiple mutant clones rising in the plasma indicated polyclonal resistance, a phenomenon previously reported (18). In this case, our model allowed estimation of the progression rate for both subclones independently (KRAS^G12D = 0.303/week; KRAS^Q61H = 0.122/week). The radically different growth rates of these resistant subclones, as well as the fact that they crossed over each other, confirmed that these variants were indeed in two different cancer cell populations, corroborating polyclonal resistance (rather than two mutations in the same subclone or in nested subclones).

Importantly, the application of the model to both CEA and cfDNA allowed estimation of the progression rate from two independent measurements, allowing us to determine how accurately the dynamics observed in the plasma would recapitulate the dynamics observed in CEA later on. Remarkably, we found that despite the limited cohort, progression rate measured by CEA with Equation (1) and progression rate measured by cfDNA with Equation (2) were significantly correlated (Fig. 5I). This not only confirmed that the mutant subclones detected in plasma were indeed major contributors to resistance, but also indicated that cfDNA profiling allows quantitative prediction of the time to progression in those patients who initially responded. Because in responders cfDNA dynamics precede CEA by several weeks (Fig. 5G, the green line preceding the black line) and the growth rates of these two curves correlate, it becomes possible to use plasma to forecast the time when we expect to observe clinical progression by RECIST v1.1 measurements. Given that these measurements are considered proportional to tumor burden (RECIST v1.1 standards: 20% increase in lesion diameter, 72.8% increase in volume), we were able to test our model against the RECIST v1.1 data to prove if predictions can be made based on either cfDNA or CEA. Specifically, we plugged in the resistance parameters estimated for each patient using CEA and cfDNA (Fig. 5F) back into Equation (1) and asked the following question: If the population at baseline is |$N( {t\; = \;0} )$|⁠, at what time |$t$| do we expect |$N( t )$| to be 72.8% larger than |$N( {t\; = \;0} )$|? This simply requires solving Equation (1) for t. We found that measurements of progression rate from CEA and cfDNA using our model (Fig. 5F) were considerably precise in predicting when progression was identified through RECIST (Fig. 5J), particularly given the low accuracy of RECIST v1.1 measurements (a manual assessment of a CT scan every 3 months). In the case of patient 1007, both CEA and cfDNA KRAS^G12D predicted the time to progression with reasonable accuracy, whereas the secondary resistant subclone KRAS^Q61H, which was slow growing, was not accurate at all because of the dominant effects of G12D that anticipated progression of several months (Fig. 5J). It is important to note that our method allows prediction of which subclone will dominate the dynamics of progression because of the estimation of each subclone's progression rate independently. The estimated parameters can then be used to calculate which subclone would lead to progression first, as well as the combined effects of multiple subclones. In the case of patient 1007, the dynamics of progression were driven by the fastest-growing subclone, KRAS^G12D. Hence, as expected, accuracy in predicting progression from cfDNA relies on identifying the dominant resistant subclones. The predictions of the time to progression for all the other responders are reported in Fig. 5K. Percent error in these estimations with respect to the time from baseline to RECIST v1.1 progression is reported in Supplementary Fig. S15. We note that our current predictions have two limitations that are not dependent on our model: (i) they are tested against RECIST v1.1 measurements, which may not represent entirely accurate estimations of the exact time to progression because they are manual assessments of CT scans that are performed only every 3 months (hence progression could have happened up to 3 months minus 1 day before), and (ii) the predictions are based on the assumption that the detected mutants in the cfDNA are responsible for the majority of the resistance. If there are other undetected resistant subclones, or a component of nongenetic resistance, our model needs such information to perform an accurate prediction. Indeed, as expected, when resistance is polyclonal but not all subclones are detected or profiled using cfDNA, the errors in the predictions are higher, as for the case of patient 1014′s subdominant MET amplification. Interestingly, however, our model allows for quantifying a posteriori the magnitude of the “unexplained” resistance and determining the growth rate of the undetected resistant subclones, the “dark matter” of resistance, depending on the actual time observed with RECIST v1.1. Despite some limitations, in several cases our predicted time to progression was remarkably accurate. Especially considering the extensive interpatient variability of clinical response and the extraordinary underlying complexity of the disease, these dynamics can be captured by a relatively simple and applicable model with relatively high accuracy.

In Fig. 5L, we illustrate the clinical impact of predictive modeling for different dynamics of resistance in each patient and for different sensitivity of cfDNA detection. Assuming an example of initial resistant frequency of 0.03% (1 every 3,300 cells, the median estimated in our cohort) and the range of progression rates in our cohort, with the model we can calculate the expected time to progression according to RECIST v1.1 criteria for each rate. Initially, the resistant population will be undetectable in the cfDNA because of biological and technical limitations (white area). Assuming weekly blood profiling, at some point cfDNA will become positive for resistance variants (blue) depending on the accuracy of the detection method (different accuracies are reported in different panels), and after enough time points we will be able to fit the model and infer the crucial parameters (frequency and progression rate) that will predict progression (red). The sooner we can detect the mutant alleles in the blood, the earlier we can forecast progression, thus creating a larger window of opportunity (yellow area) to take clinical decisions for a specific patient, such as changing treatment or adjusting treatment dynamically. We extended this illustrative analysis to the case of using CEA alone to predict resistance (Supplementary Fig. S16), as well as for different time intervals of blood sampling (Supplementary Fig. S17). The latter indicates that, when possible, at least a biweekly blood sampling provides significantly greater predictive power with respect to 4 weeks, and that beyond the 4-week interval the predictive power is limited. This may help the design of future studies.

The ability to design optimal personalized treatments relies strongly on the possibility of predicting the course of the disease in individual patients. Cancer evolution is the fundamental paradigm to understand how tumors change over time (31), and evolutionary mathematical modeling of tumor growth (32, 33) provides the theoretical framework to construct predictive models in cancer.

For the first time within a prospective phase II study, we demonstrated that the combination of longitudinal plasma biopsies and solid-tissue biopsies can be coupled with mathematical modeling of tumor evolution to anticipate tumor progression quantitatively, thus affecting future clinical decisions. We validated previous retrospective and preclinical observations (16–18, 34, 35), supporting the notion that RAS pathway aberrations clonally expand during anti-EGFR treatment. More importantly, we showed that approximately 50% of patients with mCRC considered KRAS WT, and as such eligible for anti-EGFR treatment, in fact present RAS aberrations and do not benefit from cetuximab. Our data might, at least in part, explain the observation that even in pan-RAS WT patients response and benefit from cetuximab are limited (36).

Whereas two thirds of all the abnormalities detected on baseline bloods prior to cetuximab treatment were observed in RAS genes, one third were found in genes not routinely tested in clinical practice, such as PIK3CA and HER2, proved to be involved in primary resistance (37), suggesting that extending genomic testing beyond the RAS genes might be useful for patient selection. The discrepancy in RAS status between archival material (usually primary cancer resections or biopsies) and baseline bloods prior to anti-EGFR treatment can be easily explained by a number of causes including ITH or sampling/technical errors. Whether previous chemotherapy treatments (38) or evolutionary dynamics during metastatic progression (39, 40) might have had a role in priming RAS-mutant subclones remains unknown; however, it is important to point out that our data, as well as previous observations (41–43), suggest that these clones preexisted treatment. Although further studies should address if evolutionary bottlenecks induced by previous lines of chemotherapy might have had a role in selecting RAS-mutant clones, it is intuitive to think about blood-based RAS genotyping as a rational strategy to overcome hurdles in patient selection in keeping with similar approaches used for EGFR testing in NSCLC (44).

By comparing matched tissue and liquid biopsies, we demonstrated that cfDNA captures the overall evolutionary dynamics of the disease remarkably well. However, we note that sparse cfDNA sampling, such as sampling plasma every 2 to 3 months or even less frequently, does not provide sufficient predictive power to forecast treatment resistance in individual patients due to the inherent interpatient variability of malignant evolution. On the other hand, collecting frequent longitudinal cfDNA samples from each individual patient (every 4 weeks or even more frequently) was uniquely amenable to mathematical modeling of cancer evolutionary dynamics. Combining sequential mutant frequency data, tumor burden (CEA), and evolutionary modeling allowed us to measure the dynamics of resistance in each patient and then predict the estimated time to progression. This represents a fundamental step in the clinical translation of liquid biopsies for patient care as it allows the striking variability among patients to be overcome, and patient-specific predictions to be made.

Even though our study clearly demonstrates the superiority of liquid versus tissue biopsies in providing clinically relevant information and, in keeping with other studies (45), highlights limitations of tissue biopsies in capturing spatial ITH, it also offered a unique opportunity to track temporal ITH upon cetuximab treatment concomitantly in plasma and tissues. This study also confirms that polyclonal resistance is a common feature in anti-EGFR–refractory patients and suggests the presence of a complex ecosystem ruling the emergence of “dominant” resistant clones (46). Genomic bottlenecking has been observed upon response to chemotherapy in gastroesophageal cancers (47, 48); here, we observed evolutionary bottlenecks in RAS pathway aberrations at time of disease progression, suggesting a hierarchical structure in the selection of the “fittest” resistant clones (35). We acknowledge that our study was limited to the RAS pathway and few other cancer-related genes; as such, other genetic (37) or nongenetic (49) determinants are likely to also play a role in this selection process and therefore in our predicted time to progression. We also recognize that the ability to predict time to disease progression using our mathematical model will have to be prospectively validated in future trials. In this context, ongoing trials such as CHRONOS (Rechallenge with Panitumumab Driven by RAS Dynamic of Resistance; ClinicalTrials.gov Identifier: NCT03227926] will offer the opportunity to couple theoretic modeling with therapeutic intervention in order to define the validity and clinical utility of “drug holidays” and windows of opportunity.

In conclusion, we show that combining liquid biopsies with mathematical modeling of tumor evolution allows quantitative anticipation of tumor progression, informing clinicians about timing of clinical decisions and future treatment strategies, facilitating the application of precision medicine with significant health and economic benefits for patients and health systems.

The PROSPECT-C trial (clinical trials.gov number NCT02994888) is a phase II, open-label, nonrandomized study of anti-EGFR monoclonal antibodies in patients with RAS WT, refractory mCRC. Patients who were at least 18 years old and had a World Health Organization performance status of 0 to 2 were considered eligible for this study if they fulfilled all the following criteria: (i) chemorefractory (at least two lines of chemotherapy) metastatic disease; (ii) KRAS/NRAS WT (on archival material according to hospital policy); (iii) measurable disease; and (iv) metastatic sites amenable to biopsy. Patients received cetuximab/panitumumab through the Cancer Drug Fund. Written informed consent was obtained from all patients. The study was carried out in accordance with the Declaration of Helsinki and approved by National Institutional Review Boards (National Research Ethics Service: 12/LO/0914). All participants were required to have mandatory pretreatment biopsies (2 cores), biopsies at 3 months (if PR by RECIST v1.1 criteria; 2 cores) and at the time of PD (2 cores from two suitable progressing metastatic sites). Treatment consisted of cetuximab 500 mg/m² once every 2 weeks until progression or intolerable side effects. All but one patient received the cetuximab mAb and was not entirely anti-EGFR naïve at the time of trial entry; indeed, this patient was switched to panitumumab due to a Common Toxicity Criteria for Adverse Events 3.0 Grade II allergic reaction after the first dose of cetuximab and had previously received 3 cycles of fluorouracil, oxaliplatin, and cetuximab combination with PR as neoadjuvant chemotherapy for liver resection in the context of the New-EPOC trial (50) 13 months before entering the PROSPECT-C trial.

cfDNA was extracted from EDTA anticoagulated blood within 1 hour after collection, and plasma was separated from the cells by centrifugation (1,500 × g for 15 minutes at 4°C) followed by a second centrifugation of the supernatant at 1,500 × g for 10 minutes at 4°C to remove all cell debris. If not used immediately, plasma was frozen at −80°C until further processing. cfDNA from 4 mL of plasma was isolated by the use of the Qiagen blood mini kit (Qiagen) according to the manufacturer's protocol.

The QX200 ddPCR system (Bio-Rad) was used, and all reactions were prepared using the ddPCR Supermix with no dUTTP for Probes. All PCR reactions were performed as duplex PCR using the relevant digital PCR assays for the WT and the mutation in question. Droplets were generated starting from 8 μL of cfDNA template using the QX200 droplet generator according to the manufacturer's protocols. The PCR reaction was performed in a C1000 Touch Thermo Cycler (Bio-Rad) using the following protocol: 95°C for 10 minutes followed by 40 cycles of 94°C for 30 seconds and 55°C for 1 minute, then 98°C for 10 minutes. Droplets were read in the QX200 droplet reader and analyzed using the Quantasoft software version 1.6.6.0320 (Bio-Rad). Fractional abundance (FA) was defined as follows: FA % = (Nmut/(Nmut + Nwt)) × 100), where Nmut is the number of mutant events and Nwt is the number of WT events per reaction. The number of positive and negative droplets was used to calculate the concentration of the target and reference DNA sequences and their Poisson-based 95% confidence intervals. ddPCR analysis of normal control plasma DNA (from cell lines) and no-DNA template controls were always included. Samples with very low positive events were repeated at least twice in independent experiments to validate the obtained results, as previously described (17).

The targeted panel was sequenced on a HiSeq 2500. Residual adapter sequences were trimmed with Skewer (51) and reads aligned with Burrows–Wheeler aligner. Base-level (base quality ≥ 25) coverage of amplified genomic regions was extracted and used to call variants with a Bayesian beta-binomial model (shearwater algorithm) implemented in the deepSNV package for R. Locus-specific error rates were estimated from a composite set of buffy coats. The beta-binomial model, which includes a global dispersion factor in addition to a site-specific error rate, was chosen because a simpler binomial model did not reflect the variability (overdispersion) observed in the normal samples (buffy coats). Mutations were called only when (i) coverage was at least 20,000× and (ii) the posterior probability of the mutation being a false positive was less than 5%. Overall, formalin-fixed, paraffin-embedded (FFPE)–specific background rates were low in buffy coats and FFPE samples, indicating sufficient removal/repair of FFPE-related DNA damage during sample preparation.

The Avenio panel was run in the clinically accredited Molecular Diagnostic Laboratory in the Centre for Molecular Pathology at the Royal Marsden Hospital. DNA-sequencing libraries were prepared using the Avenio ctDNA Analysis Kit (Roche), starting with 25 ng DNA, following the manufacturer's instructions, and hybridized to the Avenio Expanded Capture Kit to enrich for a panel of 77 target genes and regions. Libraries were quality-checked on an Agilent TapeStation. Sequencing was performed (150 bases, paired end) on an Illumina Nextseq 500 (High Output), 8 samples per run (approximately 100 million PE reads/sample). Data were analyzed with the Roche Avenio Custom App via a locally installed Roche server to generate variant allele frequency (VAF) and unique allele depth data (mean fold depth, 15,008 ± 2,955; unique depth, 4,710 ± 2,319). The Roche analysis pipeline supports VAF detection of SNVs to 0.5%, targeted indels and fusions to 1%, and copy-number variations over 2.3-fold with sensitivities of >99%. Variants are detectable to 0.1% VAF.

I. Chau is a consultant/advisory board member for Eli Lilly, Bristol-Meyers Squibb, MSD, Merck Serono, Bayer, Roche, and Five Prime Therapeutics. S. Rao has received other remuneration from Merck Serono. D. Watkins has received other remuneration from Amgen. N. Valeri has received honoraria from the speakers bureaus of Bayer, Merck, and Eli Lilly. No potential conflicts of interest were disclosed by the other authors.

Conception and design: K.H. Khan, D. Cunningham, B. Werner, D. Watkins, N. Starling, A. Sottoriva, N. Valeri

Development of methodology: K.H. Khan, B. Werner, I. Spiteri, J. Fernandez Mateos, N. Starling, N. Khan, N. Valeri

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): K.H. Khan, D. Cunningham, G. Vlachogiannis, I. Spiteri, J. Fernandez Mateos, A. Lampis, H. Lote, S. Hedayat, I. Chau, F. Trevisani, S. Rao, G. Anandappa, D. Watkins, N. Starling, J. Thomas, N. Khan, R. Begum, B. Hezelova, P. Proszek, J.C. Hahne, M. Hubank, C. Braconi, N. Valeri

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): K.H. Khan, D. Cunningham, B. Werner, T. Heide, A. Vatsiou, A. Lampis, I. Said Huntingford, I. Chau, N. Tunariu, G. Mentrasti, N. Starling, C. Peckitt, M. Rugge, T. Jones, P. Proszek, M. Fassan, J.C. Hahne, M. Hubank, A. Sottoriva, N. Valeri

Writing, review, and/or revision of the manuscript: K.H. Khan, D. Cunningham, B. Werner, I. Spiteri, A. Lampis, M. Darvish Damavandi, I. Chau, N. Tunariu, F. Trevisani, S. Rao, D. Watkins, N. Starling, M. Rugge, P. Proszek, M. Fassan, J.C. Hahne, C. Braconi, A. Sottoriva, N. Valeri

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): K.H. Khan, A. Lampis, M. Darvish Damavandi, S. Hedayat, G. Mentrasti, B. Hezelova, J.C. Hahne, N. Valeri

Study supervision: K.H. Khan, D. Cunningham, A. Sottoriva, N. Valeri

Other (trial physician on this study and involved in patient recruitment and patients' day-to-day care): K.H. Khan

Other (trial manager for administration the of study): A. Bryant

Other (mathematical modeling and analysis): B. Werner

This work was supported by Cancer Research UK (grant number CEA A18052), the National Institute for Health Research (NIHR) Biomedical Research Centre (BRC) at The Royal Marsden NHS Foundation Trust and The Institute of Cancer Research (grant numbers A62, A100, A101, and A159), and the European Union FP7 (grant number CIG 334261) to N. Valeri. A. Sottoriva is supported by the Wellcome Trust (202778/B/16/Z), Cancer Research UK (A22909), and The Chris Rokos Fellowship in Evolution and Cancer. B. Werner is supported by the Geoffrey W. Lewis Post-Doctoral Training fellowship. The authors acknowledge support from the National Institute for Health Research Biomedical Research Centre at The Royal Marsden NHS Foundation Trust and The Institute of Cancer Research. This work was also supported by Wellcome Trust funding to the Centre for Evolution and Cancer (105104/Z/14/Z).