Adaptive therapies that alternate between drug applications and drug-free vacations can exploit competition between sensitive and resistant cells to maximize the time to progression. However, optimal dosing schedules depend on the properties of metastases, which are often not directly measurable in clinical practice. Here, we proposed a framework for estimating features of metastases through tumor response dynamics during the first adaptive therapy treatment cycle. Longitudinal prostate-specific antigen (PSA) levels in 16 patients with metastatic castration-resistant prostate cancer undergoing adaptive androgen deprivation treatment were analyzed to investigate relationships between cycle dynamics and clinical variables such as Gleason score, the change in the number of metastases over a cycle, and the total number of cycles over the course of treatment. The first cycle of adaptive therapy, which consists of a response period (applying therapy until 50% PSA reduction), and a regrowth period (removing treatment until reaching initial PSA levels), delineated several features of the computational metastatic system: larger metastases had longer cycles; a higher proportion of drug-resistant cells slowed the cycles; and a faster cell turnover rate sped up drug response time and slowed regrowth time. The number of metastases did not affect cycle times, as response dynamics were dominated by the largest tumors rather than the aggregate. In addition, systems with higher intermetastasis heterogeneity responded better to continuous therapy and correlated with dynamics from patients with high or low Gleason scores. Conversely, systems with higher intrametastasis heterogeneity responded better to adaptive therapy and correlated with dynamics from patients with intermediate Gleason scores.

Significance:

Multiscale mathematical modeling combined with biomarker dynamics during adaptive therapy helps identify underlying features of metastatic cancer to inform treatment decisions.

The majority of cancer-related morbidity and mortality is driven by late-stage metastatic disease (1). Metastasis puts an increasing metabolic strain on the body and disrupts normal tissue function, which causes an overall weakening of the patient. Systemic cytotoxic treatments can lessen the overall disease burden, but a balance is needed between maximizing tumor kill, avoiding patient toxicity and poor quality of life, and managing the emergence of treatment resistance (2). Whilst often metastatic staging may be binary (metastases present or not), some cancers have more nuanced classifications, such as the oligometastatic state, where metastases-directed treatment may be curative (3). Treatment failure occurs due to resistant cells, which can arise from all lesions, from a subset of them, or from a single metastatic site. Given the wide variation in lesion sizes, locations, and characteristics, measuring the relevant properties of metastatic disease that can inform treatment decisions is inherently challenging. There is potential heterogeneity in the molecular make-up of a single metastatic lesion as well as between lesions. Whilst large tumors can be routinely monitored by imaging, newly colonized micrometastases may be too small to be detected. In some cancers, overall burden and disease state may be estimated with blood biomarkers, but at best this provides a coarse-grained view of systemic burden. Thus, new ways in which to measure features of metastases and interlesion and intralesion heterogeneity of metastatic cancers are urgently required.

In this work, we focus on metastatic prostate cancer. The clinical data that is routinely collected consists of blood biomarkers, imaging, and biopsies. Prostate-specific antigen (PSA) is a blood biomarker collected from peripheral blood, which is routinely used for prostate cancer screening. PSA is produced by most prostate cancer cells and is the gold-standard indicator for systemic tumor burden. Diagnostic imaging with either a bone scan, local MRI, or CT can be used to visualize the extent of metastases that most often seed to the bone, and less often to the lymph, lung, and liver. Biopsies of the primary lesion yield a Gleason score, which grades the level of similarity between benign and tumor tissue. Immunohistochemistry and other assays may be used to identify various markers in the primary lesion and potentially some larger metastatic lesions. However, each of these sources of data has its drawbacks. Biopsies can only capture a single time point for a small section of a visible and accessible lesion. Imaging, while able to capture whole lesions over multiple time points, is more spatially coarse-grained and has a threshold of detection that misses small lesions. Furthermore, bone scans for metastatic prostate cancer often capture metabolically active hotspots that can indicate the presence of a lesion, but quantifying its size is obscured by inflammation. Finally, blood biomarkers such as PSA, circulating tumor cells, and circulating tumor DNA can provide a finer temporal resolution, but only represent the total systemic burden, if not just a subset of the total burden.

PSA is the most widely used molecular marker for the diagnosis and management of a malignancy, however, absolute PSA alone is not always useful (4). Prior to therapy, it is not clear whether PSA is prognostic, however, PSA ≤ 4 ng/mL after 7 months of androgen deprivation therapy strongly predicts improved survival (5). While most studies have focused on using static metrics to predict outcome, an emerging approach that holds promise is using dynamic biomarkers instead. These markers base prognosis not on the absolute value of measurement at a single time point, but on the relative change over time. For example, PSA doubling time may stratify patients with prostate cancer likely to respond better to chemotherapy, progress to metastatic disease, or die from prostate cancer (4). In this paper we explore how PSA response dynamics can allow us to infer additional information about the characteristics of the disease that standard analysis of biomarkers cannot reveal.

To do so, we leverage data from a new treatment approach that has been gaining ground in the treatment of advanced cancers. This is an evolutionary-informed strategy, termed “adaptive therapy,” in which dosing is adjusted in response to biomarker dynamics. The aim is to prolong control over the disease by maintaining a certain tumor burden in exchange for delaying, or even preventing, the emergence of treatment resistance. Treatment breaks are given adaptively in a patient-specific manner, to prevent “competitive release,” which is when a small population of resistant cells begins to grow after therapy has removed the larger population of sensitive cells that was competitively holding resistant growth in check. Adaptive therapy has been successfully tested in several in vivo models (6–8), which has motivated the launch of clinical trials in prostate cancer, thyroid cancer, and melanoma (NCT02415621, NCT03511196, NCT03630120, NCT03543969). The most advanced of these clinical studies investigated an adaptive algorithm for the treatment of metastatic castration-resistant prostate cancer (mCRPC) using the Cyp 17 inhibitor abiraterone (9). In this study, tumor burden was monitored monthly using PSA, and treatment was withheld once patients had achieved a 50% drop from their pretreatment PSA level. There have been other trials that have used similar dose-minimizing algorithms, such as intermittent therapy (10–13). However, these studies started and stopped therapy based on either fixed time intervals, or absolute values of PSA as a threshold. In contrast, the adaptive therapy study normalized the PSA thresholds for each patient by their initial PSA levels; this allows adaptive therapy to be more sensitive to the nadir of cancer burden in each patient, which must remain high enough to prevent competitive release of resistant cells. In addition, this patient-specific approach may reveal additional insights into the tumor's dynamics when coupled with mechanistic mathematical models.

Mathematical modeling has played an important role in advancing the development of adaptive therapy (9, 14–26). These studies have mainly focused on how various features of sensitive and resistant cells affect the response to adaptive therapy, and in general the models applied to individual tumors and systemic biomarkers, rather than widely disseminated metastatic disease. We take a different approach and consider the treatment response dynamics to characterize the properties of the system of metastases. Relatively few mathematical models have been explicitly tailored to understanding metastases at all; those that have cover a diversity of aspects pertaining to metastatic disease, which can be represented simply as a burden (3), comprehensively by connecting primary tumor growth to metastatic burden (27–32), or be defined by the steps in the metastatic cascade (33, 34). Other models focus on specific aspects of metastasis, such as invasion and survival of metastatic clusters in the circulation, seeding timing, diversity, and initiation (35–38), metastatic spread and dissemination (39), the role of immune system (40), and secondary seeding (41).

In this work, we propose a framework for estimating individual and collective components of a metastatic system using tumor dynamics during the first adaptive therapy treatment cycle. We present an analysis of the longitudinal PSA dynamics of 16 patients with mCRPC undergoing adaptive androgen deprivation treatment (NCT02415621). We quantify on- and off-treatment times and investigate relationships with clinical variables such as Gleason score, the change in the number of metastases over a cycle, and the total number of cycles over the course of treatment. To gain a more mechanistic understanding of how the observed statistical trends may relate to characteristics of the underlying disease, we compare our clinical observations with an agent-based model (ABM) of a metastatic cancer. We investigate how the mode of drug action, total tumor burden, number of metastases, turnover rate, degree of resistance, and inter- and intrametastasis heterogeneity all have distinct effects on the cycle dynamics under adaptive therapy. Whilst we focus on prostate cancer, our theoretical insights may be applied more broadly. As such we propose that biomarkers based on the dynamics in the first cycle of adaptive therapy may allow us to better characterize an individual patient's metastatic disease and help us to better identify which patients are best suited to adaptive versus maximum tolerated dose therapies.

Clinical data

We utilize the clinical data from the 16 patients with mCRPC on an adaptive abiraterone therapy trial (NCT02415621). Patients in this trial were given abiraterone until their PSA dropped to 50% of the original value and there was imaging confirmation of no progression (9). Abiraterone was commenced again once the PSA returned to the original patient specific value, and then the cycle was repeated. The data used included Gleason score, monthly PSA measurements, and the number of lesions at several time points. If a bone met no longer shows active uptake on bone scan, it will be removed from the metastasis count, and lymph nodes that have reduced to < 1 cm will be removed from the metastasis count. Data for individual patients in this study are listed in Supplementary Table S1. Some analyses could only use a subset of the cohort (n = 14) due to missing information for some patients, so the number of patients for each result lists this value. Patients in this study were compared with a historical cohort with similar ranges for age, Gleason score, and treatment history. More information on the patient demographics in both cohorts can be found in (9, 42), and (21).

Computational model

Extending our previously published model (43), our computational framework simulates a system of metastases in which the metastatic lesions are modeled using multiple off-lattice ABMs. Each instance of an ABM, representing an individual lesion, is set in its own domain, but is subject to the same systemic therapy protocol. To limit the complexity, we do not explicitly model the primary tumor but focus just on the metastases and their diversity to investigate how the heterogeneity within and amongst the metastases is affected by the systemic treatment.

Basic model setup

Figure 1 shows an overview of our modeling framework. During the initial growth period, metastases are seeded with various phenotypic compositions. The example in Fig. 1A shows a system of three metastases seeded at different time points. The metastases are allowed to grow until their total burden reaches a preset value; this ends the seeding/growth period and marks the point at which metastases are discovered and the treatment phase begins. The treatment schedule follows the adaptive therapy algorithm used in the clinical trial by Zhang and colleagues (9), where PSA is used to guide treatment decisions, such that treatment is applied until the PSA is 50% of the original value and restarted on returning to the original value. In the computational model, PSA is calculated as the total burden (sum of all cells in all metastases), under the simplifying assumption that all cells equally produce PSA. The burden is monitored at 1-week intervals; this is on a faster time scale than the clinical trial to account for the smaller spatial scale used in the models. After the total tumor burden is reduced to the 50% threshold, the drug is discontinued. The metastases are allowed to grow unimpeded until reaching the original burden again, at which point, the drug is reapplied. Note that if the tumor burden does not decline below 50% due to increasing resistance, the drug will stay on for the duration. The treatment cycles repeat until progression occurs, which we define as the time at which the tumor burden exceeds 120% of the original burden, or the time at which the simulation was terminated upon reaching an imposed time limit. The treatment burden thresholds and time limits vary, and they are specified in the results.

Figure 1.

Mathematical model overview and design. A, Simulated treatment trajectory. During the seeding/growth phase, metastases are seeded and grown until reaching a specific total burden at diagnosis (Dx). The treatment period follows, during which, therapy is applied until the tumor burden reaches half of the original value, and therapy is withdrawn until the tumor burden reaches the original value again, where the treatment cycle begins anew. B, For each metastasis, individual cells are seeded in a separate off-lattice domain. Total burden used for treatment decisions is defined as the sum of all cells in the whole system of metastases. C, ABM decision flowchart for each individual cell over a time step. D, After reaching the designated burden, the treatment period begins. E, Cells are considered druggable if they lie near the outside edge of the mass (see Supplementary Methods for details), and they can only be killed by the drug if they are druggable. F, Metastases are either subjected to a cell cycle–dependent or cell cycle–independent drug, acting at the point in the decision flowchart shown.

Figure 1.

Mathematical model overview and design. A, Simulated treatment trajectory. During the seeding/growth phase, metastases are seeded and grown until reaching a specific total burden at diagnosis (Dx). The treatment period follows, during which, therapy is applied until the tumor burden reaches half of the original value, and therapy is withdrawn until the tumor burden reaches the original value again, where the treatment cycle begins anew. B, For each metastasis, individual cells are seeded in a separate off-lattice domain. Total burden used for treatment decisions is defined as the sum of all cells in the whole system of metastases. C, ABM decision flowchart for each individual cell over a time step. D, After reaching the designated burden, the treatment period begins. E, Cells are considered druggable if they lie near the outside edge of the mass (see Supplementary Methods for details), and they can only be killed by the drug if they are druggable. F, Metastases are either subjected to a cell cycle–dependent or cell cycle–independent drug, acting at the point in the decision flowchart shown.

Close modal

Metastasis seeding

Each metastasis is seeded initially with 30 cells. This allows a small initial seed with a large enough distribution of phenotypes to be well represented when we have intratumor heterogeneity. The metastases are seeded during the seeding/growth period, and we end the growth period at a specific total tumor burden NDx. The metastases may be seeded all at once or staggered. If the seeds are staggered, we determined that the simplest rate of seeding is at regular intervals Tseed; this was chosen to ensure that all metastases will be seeded during the growth period before reaching the preset total burden. We define the time interval between each metastasis seeding as |${T}_{seed}\ = \ \frac{{{N}_{Dx}}}{{k\mathop \sum \nolimits_{i\ = \ 1}^m i}}$|⁠, where m is the number of metastases, and k = 150 cells/d was calibrated within the relevant parameter space. Whilst we realize this fixed interval seeding approach is an oversimplification of how metastases actually seed, our purpose here is more to contrast scenarios where seeds might be disseminated early with those where dissemination is more continuous prior to diagnosis. Subsequent work could investigate more complex seeding scenarios.

Cell interactions and decisions

In the off-lattice ABM, which is derived from our earlier work (43), each metastasis is composed of individual tumor cells. These cells vary in proliferation rate, modeled as a rate of progression through the cell cycle. However, proliferation is also subject to each cell's interactions with its spatial neighborhood (Fig. 1B). The life cycle of a single cell is illustrated in the flow diagram shown in Fig. 1C. At each time point, a cell is subjected to stochastic death (i.e., cellular turnover). If it survives, and if the cell has enough space to proliferate (i.e., there is space to fit a cell adjacent to the dividing cell without any overlap), then it will continue progressing through the cell cycle. If there is not enough space, the cell enters a quiescent state and does not progress through the cell cycle during that time step. If the cell is mature (i.e., it has finished the cell cycle), the cell will divide.

Cell death

Cells may die from either random cell turnover (mentioned above) or from exposure to the drug after diagnosis (Fig. 1D). Importantly, while random death can occur anywhere within the tumor (see blue population in Fig. 1C and E), the death by drug exposure is assumed to only penetrate to the cells on the periphery (see black population in Fig. 1E) as in our previous model (43). For solid tumors, the high spatial density induces blood vessels to form quickly but poorly, which means they are often highly permeable and lack functional lymphatics. This microenvironment can lead to increased pressure and vessel occlusion in the tumor interior (44). To abstract and simplify these blood vessel dynamics, we assume that the tumor periphery is where drug delivery mainly occurs and therefore where treatment is applied and withdrawn. The tumor periphery is determined using the occupation of neighboring cells and grid spaces, as explained in detail in the Supplementary Methods. In prostate adaptive cancer therapy, abiraterone acetate (an antiandrogenic CYP17A inhibitor) is administered adaptively in combination with continuous administration of leuprolide (androgen deprivation therapy). However, to extend the modeling framework to consider applications of adaptive therapy beyond prostate cancer, it is important to consider a range of treatment mechanisms and each mechanism's effect on adaptive dynamics. In particular, we explore two theoretical drug mechanisms: (i) cell cycle–independent, which occurs at a specific rate at every time step and (ii) cell cycle–dependent, which acts only upon proliferating cells during cell division (see Fig. 1F). The probability of cell death per time step for the cell cycle–independent drug is:
formula
where Pd,max is the maximum probability of cell death per day, and r is the level of resistance. Here, we assume that resistant cells have a fitness cost, in that their rate of division is slower than that of a fully sensitive cell. Specifically, the level of resistance follows a linear trade-off with the proliferation rate (see range in Supplementary Table S2) and is calculated as r = 1 – p, where p is normalized from 0–1. The probability of cell death per division by the cell cycle–dependent drug is:
formula

In addition, we assume stochastic death (turnover) can occur throughout the entire tumor. In both cases, dead cells are removed after 18 hours (43).

Simulation environment

The model was programmed in the Java language (ver. 11.0.2), using the Abstract Window Toolkit for visualization. The platform for the program was an Intel Core i7–8559U processor as well as the high-performance computing cluster at H. Lee Moffitt Cancer Center and Research Institute (Tampa, FL).

Data availability

Data for individual patients in this study are listed in Supplementary Table S1. The source code for the computational model is available on Github at: https://github.com/MathOnco/multiMets. All other raw data are available upon request from the corresponding author.

The aim of this paper is to investigate whether the dynamics during the initial cycle of adaptive therapy can be used to learn about the broad structural characteristics of the underlying disease. To begin with, we analyze these cycling dynamics in adaptive prostate cancer therapy patients, focusing on timing instead of the particular shape of burden dynamics, because we are interested in the connection to metastatic characteristics rather than the biological mechanisms of drug response in prostate cancer specifically.

Regrowth times are typically longer than treatment response times in the first drug cycle

We analyze the mCRPC trial patients to observe the variability in the dynamics of the first treatment cycle, which is the period of first applying therapy, removing it upon 50% PSA reduction, and regrowth until initial PSA levels are reached. Supplementary Table S1 lists individual patient data relevant to this study. We find that there is a range of cycle times (Fig. 2A), with differences in response times (when the treatment is on), and in the regrowth times (when the treatment is off). Further, we can visualize this initial cycle in “cycle space,” by plotting the response time on the x-axis and the regrowth time on the y-axis for each patient (Fig 2B). While there is variation between patients, we note that on average, the time that the drug is applied is shorter than the time that it takes the PSA to return to the original value after the drug is withdrawn (i.e., the points mostly lie above the diagonal). To illustrate this for individual patients, we show the PSA levels over time for three patients distributed across this cycle space (Fig. 2C).

Figure 2.

Timing variability of adaptive therapy cycles in patient data and the mathematical model. A–I, Patient data (n = 16) from the mCRPC adaptive therapy trial (A–C), model results using a cell cycle–dependent drug (D–F), and model results using a cell cycle–independent drug (G–I). For each row of panels, the left column shows a timeline of the first on/off treatment cycle for all patients/simulations in the cohort. The burden is normalized at the initial value when treatment starts, and the timeline is shifted so that the zero-time point (gray line) is at the switch point from treatment on to treatment off. On-treatment cycle times (Tx ON; red) are flipped about the y-axis for visualization purposes. The middle column shows the “cycle space,” where there is a point for each patient/simulation showing the first treatment cycle, with the treatment response time on the x-axis and the regrowth time on the y-axis. The numbered boxes in the cycle space point to specific patients/simulations within the space, which are displayed in the right column as a timeline over the full course of treatment. The three specific example simulations in F and I share the same parameter sets: 1) [NDx, m, r0, dr, hi, ha] = [150,000, 1, 0.4, 0.06, 0.5, 0], 2) [NDx, m, r0, dr, hi, ha] = [150,000, 1, 0.2, 0.06, 1, 0.5], and 3) [NDx, m, r0, dr, hi, ha] = [80,000, 10, 0.4, 0.033, 1, 0.5].

Figure 2.

Timing variability of adaptive therapy cycles in patient data and the mathematical model. A–I, Patient data (n = 16) from the mCRPC adaptive therapy trial (A–C), model results using a cell cycle–dependent drug (D–F), and model results using a cell cycle–independent drug (G–I). For each row of panels, the left column shows a timeline of the first on/off treatment cycle for all patients/simulations in the cohort. The burden is normalized at the initial value when treatment starts, and the timeline is shifted so that the zero-time point (gray line) is at the switch point from treatment on to treatment off. On-treatment cycle times (Tx ON; red) are flipped about the y-axis for visualization purposes. The middle column shows the “cycle space,” where there is a point for each patient/simulation showing the first treatment cycle, with the treatment response time on the x-axis and the regrowth time on the y-axis. The numbered boxes in the cycle space point to specific patients/simulations within the space, which are displayed in the right column as a timeline over the full course of treatment. The three specific example simulations in F and I share the same parameter sets: 1) [NDx, m, r0, dr, hi, ha] = [150,000, 1, 0.4, 0.06, 0.5, 0], 2) [NDx, m, r0, dr, hi, ha] = [150,000, 1, 0.2, 0.06, 1, 0.5], and 3) [NDx, m, r0, dr, hi, ha] = [80,000, 10, 0.4, 0.033, 1, 0.5].

Close modal

A cell cycle–independent drug is needed to account for short response times and longer regrowth times

To explore the impact of different parameters on the treatment dynamics, we create a cohort of 1,000 virtual patients by randomly varying the six key model parameters shown in Table 1 and running each simulation through the growth/seeding period. The final states of each of these 1,000 sets of metastases become the initial conditions for two parallel treatment arms: a cell cycle–dependent drug and a cell cycle–independent drug. For both drugs, we find variation across the simulations (Fig. 2D-I). Under the cell cycle–dependent drug treatment (Fig. 2DF), the response and regrowth times are nearly equal (10.2 ± 6.9 and 11.1 ± 6.2 weeks, respectively, for mean and SD). This can be explained by the fact that death can only occur as fast as division when the drug is on, while regrowth occurs at that same rate of division. There are also some simulations that suggest the possibility of slow response rates followed by fast regrowth. This variability is due to stochastic effects, spatial effects, and parameter effects that are explored later. In contrast, the cell cycle–independent drug kill (Fig. 2GI) on average has shorter response times and longer regrowth times (3.5 ± 2.1 and 14.0 ± 5.9 weeks, respectively).

Table 1.

Ranges of parameters used in the mathematical model.

ParameterDescriptionRange
NDx Total number of cells at diagnosis 10 × 103 to 150 × 103 
M Number of metastases seeded 1–10 
r0 Initial mean drug resistance 0–0.4 
dr Random death rate 0 × 10−2 to 6 × 10−2/day 
hi Intermetastasis heterogeneity 0–1 
ha Intrametastasis heterogeneity 0–1 
ParameterDescriptionRange
NDx Total number of cells at diagnosis 10 × 103 to 150 × 103 
M Number of metastases seeded 1–10 
r0 Initial mean drug resistance 0–0.4 
dr Random death rate 0 × 10−2 to 6 × 10−2/day 
hi Intermetastasis heterogeneity 0–1 
ha Intrametastasis heterogeneity 0–1 

A visual comparison of the cycle space for the cell cycle–dependent drug (Fig. 2E) and cell cycle–independent drug (Fig. 2H) implies that the cell-independent drug is more representative of the skewed relationship between response and regrowth times seen in the patient dynamics (Fig. 2B). We hypothesize that this is because the cell cycle–independent drug model is a closer representation of the drug action of abiraterone, which acts to stop production of testosterone, a necessary resource for most prostate cancer cells. Example simulations are shown (Fig. 2F and I) and demonstrate the representative differences over different cycle lengths. Given its superior ability to describe the patient data, we focus on the cell cycle–independent drug application in the remainder of this paper.

Smaller metastases have shorter drug cycles

We next look at how different parameters affect the cycle times, and in doing so, keep other parameters fixed. These represent idealized scenarios that may not represent realistic tumors but allow us to explore effects of individual parameters in isolation. To gain an understanding of how the metastatic distribution affects treatment cycles, we look at how cycles are affected by the total tumor burden and the number of metastases present at the start of treatment. We create a new virtual patient cohort by randomly sampling 300 values of burden and initial number of metastases (from the ranges in Table 1), while keeping the other parameters fixed. First, we assume that drug resistance is 0, cell turnover is 0, and the tumors are completely homogeneous. We observe that the shortest response times (Fig. 3A) happen when there are many metastases and a lower total burden, while the longest times occur when there is a large burden and few metastases. The regrowth times follow a similar pattern but with slightly longer times (Fig. 3B). We isolate several examples in this space (labeled 1–4). Examples 1 and 2 both have 2 metastases and either a high or low burden, respectively, while examples 3 and 4 both have 9 metastases with either a high or low burden, respectively (Fig. 3C, left). Example 1 has the longest treatment cycles, examples 2 and 3 have similar cycles, and example 4 has the shortest treatment cycles (Fig. 3C, right). The reason for this difference lies in the spatial constraints on tumor growth. A large total burden split into a few metastases results in a few large lesions, while the same burden split into many metastases results in many smaller lesions. With a larger burden per metastasis, the perimeter-to-area ratio decreases, slowing the response and regrowth times, because the larger metastases must lose and regrow a larger number of cells to reach the 50% and 100% treatment switching thresholds.

Figure 3.

The effect of total tumor burden and number of metastases on adaptive therapy cycling. A and B, An array of 300 random parameter combinations is simulated and is shown colored in grayscale according to drug response time (A) and regrowth time (B), which both use the same array of values and have added noise for easier visualization. Isolines are shown at 5k, 10k, 20k, 40k, and 80k cells per metastasis, with increasing line weights. C, Example simulations with parameter combinations 1 through 4 (as labeled in A and B). The green tumors represent the multiple lesions pretreatment, while the plots show the dynamics of the lesions during the first cycle of adaptive therapy. The total burden is in red (drug on) and black (drug off), and the individual metastases are plotted in green. The light gray rectangle shows our chosen threshold of metastatic invisibility at 10k cells. D, The same simulations from A and B are projected into cycle space and colored according to their combination of total burden and number of metastases using the background colormap. The histograms show the frequency distribution of each axis separately. E, We plot the same points from D colored according to the average number of cells per metastasis for each simulation. [NDx, m, r0, dr, hi, ha] = [10,000–150,000, 1–10, 0, 0, 0, 0]. Animations of example simulations are shown in Supplementary Data, Movie S1.

Figure 3.

The effect of total tumor burden and number of metastases on adaptive therapy cycling. A and B, An array of 300 random parameter combinations is simulated and is shown colored in grayscale according to drug response time (A) and regrowth time (B), which both use the same array of values and have added noise for easier visualization. Isolines are shown at 5k, 10k, 20k, 40k, and 80k cells per metastasis, with increasing line weights. C, Example simulations with parameter combinations 1 through 4 (as labeled in A and B). The green tumors represent the multiple lesions pretreatment, while the plots show the dynamics of the lesions during the first cycle of adaptive therapy. The total burden is in red (drug on) and black (drug off), and the individual metastases are plotted in green. The light gray rectangle shows our chosen threshold of metastatic invisibility at 10k cells. D, The same simulations from A and B are projected into cycle space and colored according to their combination of total burden and number of metastases using the background colormap. The histograms show the frequency distribution of each axis separately. E, We plot the same points from D colored according to the average number of cells per metastasis for each simulation. [NDx, m, r0, dr, hi, ha] = [10,000–150,000, 1–10, 0, 0, 0, 0]. Animations of example simulations are shown in Supplementary Data, Movie S1.

Close modal

As such, the cycle times do not depend on the number of metastases or total burden, but the average size of the individual metastases (see contour lines of constant metastases size in Fig. 3A and B). To illustrate this, we plot the cohort in cycle space and color each parameter combination by the number of metastases and the total tumor burden. This shows that the full range of metastatic number and size can be projected along a single axis in cycle space, where certain parameter combinations give the same cycling dynamics (Fig. 3D). When we instead color by the average number of cells per metastasis (Fig. 3E) there is a clear pattern showing that the average size of the metastases determines the position in cycle-space. To sum up, larger metastases have longer cycles, and smaller metastases have shorter cycles.

Drug resistance and turnover affect treatment cycles

Next, we look at how treatment cycles are affected by death parameters in a new cohort of virtual patients with a range of drug resistance and cell turnover rates (see Fig. 4; Table 1). In this cohort, two homogenous metastases are grown to a combined burden of 80,000 cells prior to treatment. The results show that drug response time is shortest when resistance is low and turnover is high, while the response time is longest when the opposite is true (Fig. 4A). Generally, regrowth times are quicker when resistance and turnover are low (Fig. 4B), and slower when the opposite is true. The examples (1–4) examine the mechanisms that cause these trends (Fig. 4C). The effects observed with resistance are straightforward: highly resistant examples (3 and 4, blue lines) will respond slower to the drug (longer drug application periods), and also take longer to regrow (longer drug-free periods) because of the assumed resistance-proliferation trade-off. The effect of turnover is slightly more complicated. Examples 1 and 3 have an increased number of proliferating cells (green cells), due to increased open space from high random death. This increases drug penetration and shortens the drug application times. During the drug-free regrowth period, random cell death slows the net growth rate, increasing regrowth times. These effects contribute to how different combinations of these parameters affect a tumor's cycle times during adaptive therapy (Fig. 4D). In the previous section, we noted that the overall cycle time is mostly determined by the tumor size. Given that these tumors are essentially the same size when treatment begins, we see that the combination of drug resistance and turnover contributes to further variation around that position within the cycle-space, which is skewed slightly toward shorter response times. In addition, the combination of effects on growth rate by these parameters leads to different net growth rates, which can be observed by coloring the same points from Fig. 4D in terms of tumor age at diagnosis (Fig. 4E). The older metastases take longer to reach the same size and have longer treatment cycles while the younger metastases have shorter treatment cycles.

Figure 4.

The effect of the level of preexisting drug resistance and turnover on adaptive therapy cycle dynamics. A and B, An array of 300 random parameter combinations is simulated and is shown colored in grayscale according to drug response time (A) and regrowth time (B). C, Example simulations with parameter combinations labeled in A and B are shown colored according to their phenotypic state (P, proliferating; Q, quiescent; D, druggable; R, random death) and over the first cycle of adaptive therapy. The dynamics show the total burden in red (drug on) and black (drug off) and the individual metastases in green. The light gray rectangle shows our chosen threshold of metastatic invisibility at 10k cells. D, The same virtual cohorts in A and B are plotted in cycle space and colored according to their combination of sensitivity and random death rate parameters using the colormap from A and B. E, The same points from D are colored according to the age of the metastases prior to treatment for each simulation. [NDx, m, r0, dr, hi, ha] = [80,000, 2, 0–0.4, 0–6×10−2, 0, 0]. Animations of example simulations are shown in Supplementary Data, Movie S2.

Figure 4.

The effect of the level of preexisting drug resistance and turnover on adaptive therapy cycle dynamics. A and B, An array of 300 random parameter combinations is simulated and is shown colored in grayscale according to drug response time (A) and regrowth time (B). C, Example simulations with parameter combinations labeled in A and B are shown colored according to their phenotypic state (P, proliferating; Q, quiescent; D, druggable; R, random death) and over the first cycle of adaptive therapy. The dynamics show the total burden in red (drug on) and black (drug off) and the individual metastases in green. The light gray rectangle shows our chosen threshold of metastatic invisibility at 10k cells. D, The same virtual cohorts in A and B are plotted in cycle space and colored according to their combination of sensitivity and random death rate parameters using the colormap from A and B. E, The same points from D are colored according to the age of the metastases prior to treatment for each simulation. [NDx, m, r0, dr, hi, ha] = [80,000, 2, 0–0.4, 0–6×10−2, 0, 0]. Animations of example simulations are shown in Supplementary Data, Movie S2.

Close modal

Seeding asynchronicity and intermetastasis heterogeneity affect tumor size distributions

We previously assumed that all metastases were seeded simultaneously and that metastases were phenotypically homogeneous. This means that despite some stochasticity, all tumors within the system of metastases were roughly the same size. However, this is unlikely. No correlation between the initial PSA and the number of metastases in both the adaptive therapy and historical cohorts was found (Supplementary Fig. S1A), but we can track even invisible metastases in the model system (Supplementary Fig. S1B). Here we investigate the impact of having differently sized metastases, particularly those below the threshold of detectability. Metastases may be undetectable due to a slow growth rate or late seeding; therefore, we compare two possible scenarios for metastatic seeding: (i) homogeneous sequential seeding and (ii) heterogeneous simultaneous seeding. For (i), we seed the first tumor at time 0 and the remaining tumors are seeded at a constant rate throughout the growth period prior to treatment. All tumors are the same phenotype, delayed in time. For (ii), we seed all tumors simultaneously at time 0. In this case, each lesion is different, composed of a single phenotype of cells drawn from a range of proliferation rates (and therefore a range of treatment sensitivities). The amount of this “intermetastasis heterogeneity” can be set between 0 (all lesions are the same) and 1 (proliferation rates are selected from the full range of the values permitted for proliferation from Supplementary Table S2). All other parameters are fixed (Table 1).

To compare with the previous results, we show simulations using homogenous tumors seeded simultaneously across a range of total burden and number of metastases (Fig. 5, left column). We plot the number of visible metastases within this parameter space. As above, because all metastases are seeded simultaneously and are homogenous, they are similar in size at time of therapy. Therefore, we find that generally either all of the metastases are visible, or all of them are below the threshold and therefore not visible (Fig. 5A). Similarly, for sequential seeding we find that when there are only a few metastases with a large burden, they are all visible, and when there are many metastases making up a small burden, they are all undetectable (Fig. 5B). However, with a larger burden spread amongst more metastases of various sizes, some of the late-seeded, smaller metastases may not be observed. A similar pattern occurs with intermetastasis heterogeneity (Fig. 5C). We find here that the high burden, high number of metastases portion of the plot has expanded substantially so that only the single metastases or several very large metastases are completely accounted for, whilst most simulations have at least some undetectable metastases.

Figure 5.

Comparing effects of seeding timing and intermetastasis heterogeneity on metastasis visibility and cycle times. A–C, Metastases are either seeded simultaneously (A), sequentially (B), or with heterogeneous phenotypes simultaneously (C). For each case, an array of values for total tumor burden and number of metastases are plotted colored according to the percentage of visible metastases. DF, Spatial distributions of several example metastases are shown prior to treatment, and temporal dynamics of individual metastasis sizes are shown over one cycle of adaptive therapy. GI, The same simulations from the arrays in AC are plotted in cycle space using the same background colormap that designates the combination of number of metastases and total number of cells. [NDx, m, r0, dr, hi, ha] = [10,000–150,000, 1–10, 0, 0, 0, 0]. Threshold of visibility is set at 10,000 cells. Animations of example simulations are shown in Supplementary Data, Movie S3.

Figure 5.

Comparing effects of seeding timing and intermetastasis heterogeneity on metastasis visibility and cycle times. A–C, Metastases are either seeded simultaneously (A), sequentially (B), or with heterogeneous phenotypes simultaneously (C). For each case, an array of values for total tumor burden and number of metastases are plotted colored according to the percentage of visible metastases. DF, Spatial distributions of several example metastases are shown prior to treatment, and temporal dynamics of individual metastasis sizes are shown over one cycle of adaptive therapy. GI, The same simulations from the arrays in AC are plotted in cycle space using the same background colormap that designates the combination of number of metastases and total number of cells. [NDx, m, r0, dr, hi, ha] = [10,000–150,000, 1–10, 0, 0, 0, 0]. Threshold of visibility is set at 10,000 cells. Animations of example simulations are shown in Supplementary Data, Movie S3.

Close modal

Intriguingly, Fig. 5DF shows that while both homogeneous sequential seeding and heterogeneous simultaneous seeding results in a distribution of sizes of metastases, some of which are below the threshold of detectability, the changes of the tumor sizes over a cycle of adaptive therapy can be informative. Several examples in this parameter space are displayed spatially prior to treatment and the normalized sizes of metastases are shown over one cycle of adaptive therapy (Fig. 5DF). With simultaneous, homogenous seeding all tumors are nearly the same size prior to the treatment cycle as after the treatment cycle (Fig. 5D). In the sequentially seeded tumors (Fig. 5E), whilst there is size variation amongst the lesions prior to treatment, they also return to a similar size distribution after treatment. There are also different sized tumors within the simultaneous heterogeneous metastases (Fig. 5F). The most significant difference between sequential seeding and heterogeneous seeding is that the most resistant tumors will not respond to treatment. Therefore, while the larger more sensitive tumors dominate the dynamics at the start of the treatment cycles, smaller more resistant tumors will grow unimpeded during treatment and may contribute more to the burden over the course of the cycle (e.g., Example 4 in Fig. 5F).

Furthermore, we find that the size and phenotype heterogeneity amongst the metastases affects the position in cycle-space (Fig. 5GI). Ultimately, when the cycles for the sequentially seeded tumors are plotted in cycle-space (Fig. 5H), we see little difference from the baseline, but some of the cycle times have shifted (compared with Fig. 5G). The metastases in Examples 2 and 4 (Fig. 5H) are grown to a small burden and have very little difference between the oldest and the youngest metastasis, compared with those in Fig. 5G. However, for Examples 1 and 3 that are grown to a larger total burden, the earlier seeding results in a more significant size difference of the largest tumors, which is enough to elongate the cycle times for that set of metastases. We therefore conclude here that it is not actually the average size of the metastasis that affects the cycle time but rather the largest tumor dominates the dynamics. For the heterogeneous metastases (Fig. 5I), the largest tumor still dominates the cycle times in much the same way if there is a large difference in size of the metastases (Examples 1 and 2). However, when more resistant cells make up a significant portion of the total burden (Examples 3 and 4), cycle times are elongated for both response and regrowth, due to slow-growing drug-resistant cells.

Tumor heterogeneity affects treatment response

In the previous sections, we showed that drug resistance affects cycle times. Further, we showed that heterogeneity amongst the metastases (intermetastasis heterogeneity) affects cycle times. Here, we examine how different types of heterogeneity, both intermetastasis heterogeneity and intrametastasis heterogeneity (heterogeneity within a metastasis), affect the outcome and if a continuous dose would be a better option for each virtual patient in the cohort of 1,000 simulations used in Fig. 2C–F (where all parameters were varied). For this question, we tested an adaptive therapy schedule against continuous treatment to consider which situations would respond better with the aim of tumor control as opposed to tumor elimination. We find that most of the simulations result in an eventual recurrence, but in 35% of the cohort, continuous treatment leads to tumor elimination. The tumors that are eliminated, on average, have both low intermetastasis and intrametastasis heterogeneity compared with those that recur (Fig. 6A). Of those that do recur, we compare their time to recurrence (TTR) on each treatment schedule with respect to measured heterogeneity prior to treatment (Fig. 6B and C). We plot the TTR for adaptive therapy on the x-axis and for continuous treatment on the y-axis, so that for each pair of treatments a point on the line means that there is no difference in TTR for continuous treatment and adaptive therapy, while a point above the line favors continuous treatment and a point below the line favors adaptive therapy. The majority (82%) of simulations lie below the line, favoring an adaptive therapy schedule. With some exceptions, the cases with the most intermetastasis heterogeneity respond better to continuous treatment (Fig. 6B) while the cases with more intrametastasis heterogeneity respond better to adaptive therapy (Fig. 6C). To understand why, we consider two extreme scenarios that are initialized with maximum intermetastasis heterogeneity with no intrametastasis heterogeneity and maximum intratumor heterogeneity with no intermetastasis heterogeneity. When the tumors with intermetastasis heterogeneity (Fig. 6D) are treated continuously, the total burden and the number of metastases decreases, but ultimately, the total burden starts to increase again as the most sensitive tumors are killed off and the more resistant tumors take over the dynamics. If instead adaptive therapy is applied, partially responsive tumors continue to grow when they could have been killed by ongoing treatment. In this scenario, the benefit of having sensitive cells to suppress the growth of resistant cells is not applicable because they are not competing for the same space, being in different lesions. However, with intratumor heterogeneity (Fig. 6E), a mixture of phenotypes within each metastasis do compete for the same space, so that adaptive therapy can delay treatment failure. In contrast, when continuous treatment is applied to this patient, there is competitive release of resistant phenotypes, which results in earlier recurrence. Ultimately, these results show that a continuous treatment or adaptive therapy treatment strategy may be preferred to extend TTR for patients with high intermetastasis or intrametastasis heterogeneity, respectively. However, for patients with low heterogeneity, the challenge is that cure may be possible using continuous treatment, while the best adaptive therapy can do is continuous control. All in all, adaptive therapy offers less risk: if cure fails, then adaptive therapy is preferable to extend TTR, and the time gained by using adaptive therapy when it is optimal far exceeds the time lost when using adaptive therapy when continuous treatment is optimal.

Figure 6.

Comparing responses to continuous treatment (CT) versus adaptive therapy (AT) in the model cohort (n = 1,000). A, Histograms of simulations with different values of intermetastasis heterogeneity (top) and intratumor heterogeneity (bottom) that respond in a complete (cure) or incomplete (recurrence) response to continuous treatment. B and C, TTR under continuous treatment and adaptive therapy for each simulation, colored according to intermetastasis heterogeneity (B) and intratumor heterogeneity (C). Gray line indicates when there is no difference in the TTR between the two treatments. D and E, Examples of metastases with intermetastasis heterogeneity (D) and intratumor heterogeneity (E) comparing continuous treatment (top) and adaptive therapy (bottom). Each shows systemic burden and individual metastases over time (top), the number of visible metastases (>10k cells) over time (middle), and a histogram of the metastases size distribution, with each colored according to its average phenotype (bottom). The corresponding spatial distribution of each metastasis is shown below each bar. Parameters used for the example simulations are [NDx, m, r0, dr, hi, ha] = [100,000, 5, 0, 0, 1, 0] and [NDx, m, r0, dr, hi, ha] = [100,000, 5, 0, 0, 0, 1] for intermetastatic heterogeneity and intrametastatic heterogeneity, respectively. Animations of example simulations are shown in Supplementary Data, Movie S4.

Figure 6.

Comparing responses to continuous treatment (CT) versus adaptive therapy (AT) in the model cohort (n = 1,000). A, Histograms of simulations with different values of intermetastasis heterogeneity (top) and intratumor heterogeneity (bottom) that respond in a complete (cure) or incomplete (recurrence) response to continuous treatment. B and C, TTR under continuous treatment and adaptive therapy for each simulation, colored according to intermetastasis heterogeneity (B) and intratumor heterogeneity (C). Gray line indicates when there is no difference in the TTR between the two treatments. D and E, Examples of metastases with intermetastasis heterogeneity (D) and intratumor heterogeneity (E) comparing continuous treatment (top) and adaptive therapy (bottom). Each shows systemic burden and individual metastases over time (top), the number of visible metastases (>10k cells) over time (middle), and a histogram of the metastases size distribution, with each colored according to its average phenotype (bottom). The corresponding spatial distribution of each metastasis is shown below each bar. Parameters used for the example simulations are [NDx, m, r0, dr, hi, ha] = [100,000, 5, 0, 0, 1, 0] and [NDx, m, r0, dr, hi, ha] = [100,000, 5, 0, 0, 0, 1] for intermetastatic heterogeneity and intrametastatic heterogeneity, respectively. Animations of example simulations are shown in Supplementary Data, Movie S4.

Close modal

Changes in visible metastases indicate intermetastasis variability

Having characterized how different tumor parameters affect the adaptive therapy cycling dynamics, we apply our novel understanding to the patients in the clinical trial (9). In the above results, we have examined how tumor parameters in isolation or in pairwise combinations alter the adaptive therapy cycling dynamics. However, metastatic disease is a complex system and patients are likely to differ in multiple parameters at once. In the final part of this paper, we therefore aim to relate parameter combinations back to the patients in the trial using the Gleason score (45) of the primary lesion, which correlates with a number of clinically important features such as outcome. In this cohort, we observe three distinct categories of response, which we will refer to as low Gleason (6–7), mid Gleason (8), and high Gleason (9–10). Over the first cycle of adaptive therapy, patients with low Gleason scores have either the same or fewer metastases at the end of the cycle than at the start of the cycle (Fig. 7A). In contrast, all patients with mid Gleason scores have no change in the number of metastases over the first cycle, and patients with high Gleason scores have either the same or an increase in the number of metastases. Further, the mid Gleason patients are seen to have a higher average number of cycles over the full treatment course, while low and high Gleason patients have only a few (Fig. 7B). Finally, the patients with the mid Gleason scores tend to have the fastest cycle times (Fig. 7C).

Figure 7.

Comparing patient metrics (n = 14) grouped by Gleason score (low = 6–7, mid = 8, and high = 9–10) to simulation data. Patient data is shown in AC and G, and model data is shown in DF and H. A–C, Relationship between Gleason score and the change in the number of metastases (A), the total number of treatment cycles (B), and the patient's location in cycle space (C). Note that two patients were left out of this analysis due to one missing Gleason score and one with insufficient data on the number of metastases. DF, Model data from Fig. 2E and F classified into a Gleason category based on the change in the number of metastases, the number of cycles, and the timing of their cycles (algorithm explained in the Supplementary Methods; patients that were unable to be classified are not shown). Large diamonds represent the patients and simulations that best fit the average metrics of each Gleason group. G and H, Treatment dynamics of patients best representing the average of each Gleason group (G) and simulations predicted by the model for the best fits (H). These correlate with the large diamonds in AF.

Figure 7.

Comparing patient metrics (n = 14) grouped by Gleason score (low = 6–7, mid = 8, and high = 9–10) to simulation data. Patient data is shown in AC and G, and model data is shown in DF and H. A–C, Relationship between Gleason score and the change in the number of metastases (A), the total number of treatment cycles (B), and the patient's location in cycle space (C). Note that two patients were left out of this analysis due to one missing Gleason score and one with insufficient data on the number of metastases. DF, Model data from Fig. 2E and F classified into a Gleason category based on the change in the number of metastases, the number of cycles, and the timing of their cycles (algorithm explained in the Supplementary Methods; patients that were unable to be classified are not shown). Large diamonds represent the patients and simulations that best fit the average metrics of each Gleason group. G and H, Treatment dynamics of patients best representing the average of each Gleason group (G) and simulations predicted by the model for the best fits (H). These correlate with the large diamonds in AF.

Close modal

Next, we wanted to assign a Gleason score to the virtual patients in the cohort of 1,000 simulations used in Fig. 2CF (where all parameters were varied). While we have no equivalent to the Gleason score in the computational model, we can match simulations with a putative Gleason score based on comparative metrics to the patient data from Fig. 7AC, which include the change in number of metastases over the first cycle, the total number of cycles over the course of therapy, and the response and regrowth time of the first cycle (see Supplementary Methods for details). The distributions and individual values for each metric in the simulated cohort classified by Gleason group are shown in Fig. 7DF. We were able to classify 48.9% of the simulations with 7.3%, 18.4%, and 23.2% in the low, mid, and high Gleason groups, respectively. The decision tree for classification and distribution of simulation data that fits in each group, including those that did not fit in a group is shown in Supplementary Fig. S2. The patient treatment trajectories for best fit to the average metrics are shown in Fig. 7G. We observe that the patient with a low Gleason score has a decrease in the number of metastases over the first cycle and fails shortly after that. The patient with a mid Gleason score has no change in the number of metastases during the first cycle and continues with six more relatively quick cycles. The patient with the high Gleason score has an increase in metastases during the first cycle and was taken off the study at that point. The model-predicted treatment trajectories for the best fitting simulation for each Gleason group are shown in Fig. 7H. The simulations with both high and low Gleason classifications had longer cycling times during the first and subsequent cycles, while the simulation with the mid Gleason classification had faster cycles with no change in the number of metastases over the first cycle.

While these patients and simulations represent the best fit to the average of the data, it is important to note that there is variation amongst key qualitative features within the group like outcome. In other words, many cycles does not mean longer TTR, but the model suggests that it may mean smaller lesions. We examined what other metastatic features our model fits could tell us about the virtual cohort (Supplementary Fig. S3). We find a noticeable difference in heterogeneity between the Gleason groups. We calculate both intermetastasis and intrametastasis heterogeneity just prior to treatment as the SD of the mean phenotype for each lesion and the mean of the SD of phenotypes for each lesion, respectively. The mid Gleason groups are found to have the lowest intermetastasis heterogeneity and the highest intratumor heterogeneity. This result illustrates how integrating the data from adaptive therapy cycling with mathematical modeling may help us to better characterize a patient's disease.

Here, we introduced a multiscale model to map tumor dynamics onto cycle space, showing that patient-specific response and regrowth times can imply disease characteristics like metastatic distribution, resistance, and optimal therapy strategies. Characterization of metastases is inherently challenging due to dispersion of metastatic seeds to many sites and the sparsity of data collected beyond just presence or absence of metastases. However, vigilant observation during a single adaptive therapy cycle could provide new knowledge. With further analysis and validation of the metrics introduced here, they may add supporting evidence to existing biomarkers and aid in personalizing treatment for patients with metastatic cancer.

To root the model to the clinical setting, we have compared simulation results with patient data from the first pilot trial of adaptive therapy in mCRPC (9, 21). But given the small sample size, we must not overinterpret the results (46). This was an exploratory examination of subjects each with a unique set of correlated data, in which we asked whether there were certain characteristics that defined a group. This data is unique in that we have correlated data for each patient containing primary tissue information (Gleason score), and longitudinal data for both burden (PSA) and the number of metastases. Similar studies with intermittent therapy have been done in prostate cancer (47, 48), but the available data did not contain metastasis information, which is essential for this study. In fact, most studies do not collect detailed information of metastases beyond presence or absence, or provide correlated longitudinal data, from which we can extract useful dynamic relationships. The general aim for this work was to understand basic structural and temporal processes leading to the observed dynamics of total tumor burden and individual metastases in solid tumors in general. In defining the parameters for the model, we built reasonable bounds around realistic biological quantities and related them to observable clinical metrics of metastatic cancer in general. However, these bounds are not necessarily known or measurable for all parameters, and some parameter combinations may not exist in prostate cancer specifically due to trade-offs or other biological limits. If we assume that the model can give reasonable mapping from input parameters to the output values of interest (change in metastases over the first cycle, number of cycles, and on/off cycle times), then when output values fell far enough outside of the assumed variability from the given patient data, we assumed that the input parameters were not representative of a group. It is possible with the small sample size that this limits the observed variability and unfairly cuts out some of the possible parameter combinations. So, the results here, in application to prostate cancer, should be interpreted as hypothesis-generating for future studies and a call for more meaningful and correlated data collection and availability.

During a single cycle of adaptive therapy consisting of a drug application period and a regrowth period, we can delineate multiple parameters of the metastatic system. We found that with the same initial burden, more metastases resulted in quicker adaptive therapy cycle times. However, it is not the number of metastases that matters, but the largest lesions that dominate the dynamics. Larger metastases have longer cycles, while smaller metastases have quicker cycles, similar to our previous findings where more spread out tumor seeds had faster cycles than larger patches of cells in (18). We also found that more drug resistance slows the treatment cycles, as expected with a slower response time, but also slows regrowth time due to the assumed trade-off (49). While there is an argument for trade-offs in certain cancers and in microenvironments with limited resources, it may not always be an accurate assumption in all cases (50, 51). However, in the prostate cancer adaptive therapy trial, patients needed to achieve at least 50% decline from their pretreatment PSA on abiraterone, so a significant portion of sensitive cells must be present initially to allow continued treatment response over multiple cycles. This argues for a proliferative advantage of sensitive cells over resistant cells during the growth period. In addition, we found that a faster cell turnover speeds up drug response time by increasing the surface area to aid in drug penetration but slows the regrowth time by reducing the net proliferation rate. The latter result was also found in a nonspatial model (19), where turnover increased total cycle times, which are generally dominated by long regrowth periods and shorter response times.

It is also important to recognize the limitations of imaging for characterizing metastases. A single snapshot in time may not capture smaller, later-seeded lesions. However, variation between metastases with regards to size and phenotype may be captured in the cycle dynamics. We found that the cycle dynamics are dominated by the largest tumor, but if enough resistant tumors exist in a large enough proportion, the cycle times will increase. Further, the change in the number of metastases over the cycle rather than just a single time point may yield more information. We found that change in metastases over the course of a cycle may indicate that there are differences amongst the metastases, which could be associated with the Gleason score of the primary tumor. Both lower Gleason scores (6–7) and higher Gleason scores (9–10) show more cases where the number of metastases decrease and increase, respectively, during the first cycle of adaptive therapy. These two groups also tend to have very few total cycles that are longer when compared with mid Gleason scores that on average have more cycles that are shorter duration. Classification of the model results into Gleason groups in this model revealed it may be due to differences in phenotypic heterogeneity, with respect to proliferation rate and treatment resistance. More intermetastasis heterogeneity was found in low and high Gleason groups, which ultimately causes treatment failure due to failure from specific lesions, and more intratumor heterogeneity was found in mid Gleason groups, which allows for local competition amongst phenotypes within a single lesion that is necessary for the success of adaptive therapy. It has otherwise been shown that mid Gleason scores are unique in that the cells may be in a transitional state (52), which may aid in treatment response. However, the increased response of mid Gleason scores to adaptive therapy does not necessarily correlate with outcome. While the standard-of-care continuous treatment does show the expected negative correlation between Gleason score and time to progression, in contrast to adaptive therapy, where the time to progression is higher overall and seems less correlated with Gleason score (Supplementary Fig. S4). In this analysis, we have used large-scale clinical features from patients with known Gleason scores to stratify systems of computational lesions with known cell-scale phenotypic characteristics into corresponding Gleason groups. However, it should be noted that the number of patients in the clinical trial is small; therefore, these results would need to be further explored for validation. Heterogeneity assessment for metastatic prostate cancer could include blood biomarker or histology measures for androgen receptor expression and mutations, neuroendocrine activity expression, Cyp17 expression, and immunologic components, as well as metabolic imaging modalities, including prostate-specific membrane antigen PET-CT (53–57).

Here, we have considered a minimal model of a system of metastases to gain a basic understanding of individual and systemic dynamics during a single cycle of adaptive therapy. Our aim was to capture the most generic structural elements of metastases during treatment while comparing to the mCRPC trial for guidance. In preparing a baseline model, even though is is quite complex, we did not include multiple important elements such as the effect of early or late seeding (58–60), metastases seeding metastases (61), differences in metastasis cluster size (62, 63), dormancy (64, 65), diversity of metastatic potential (62, 66, 67), drug-induced resistance at the metastatic site (68, 69), and stromal interactions in metastatic niches (52). However, our minimal model was calibrated to the dynamics of the patients in the mCRPC trial, which was limited by the available data, but additional elements could be investigated in future studies to gauge their impact on trial outcomes. In addition, the model time and space scales were both small compared with those in the mCRPC trial. The size scale was limited by the domain size, which was kept small for simplicity and computational efficiency, and the timescale was reduced in turn. Therefore, cycling times were smaller, and we needed to use 1-week intervals to check burden and make treatment decisions rather than the 4-week intervals in the clinical trial. These scale issues could be diminished by using a different, coarser model system, but would likely come at the expense of not having cell-scale spatial interactions that may be important mechanisms affecting the outcome, due to heterogeneity and competition within a lesion.

Critically, combining available metrics on different scales of resolution may aid in a better characterization of disease at this late stage when metastases may be invisible and challenging to quantify, as each modality has limited scope in isolation. Blood biomarkers are not perfect, yet they remain a cheap viewpoint of a systemic measure of tumor burden over time. Imaging also has resolution limits, but changes over the cycle may help inform treatment choice. Ultimately, there are only two reasons why continuous, maximum dose treatment succeeds: the metastases do not have preexisting heterogeneity, or the harsh treatment prevents the resistant phenotypes from taking over. However, when treatment fails, it can do so in different ways. Most simulations that resulted in recurrence favored an adaptive therapy schedule. However, some cases with the most intermetastasis heterogeneity responded better to continuous treatment. Perhaps a molecular biomarker signature could be developed to quantify tumor heterogeneity in a biopsy and/or blood sample that when combined with these results could be used to help stratify patients into optimal treatment groups. Currently, cancer staging provides only a binary view of metastatic presence or absence. But to make the most out of the limited data we can collect, we need to integrate information across spatial scales (different subpopulations/different metastases) as well as across time. Multiscale mathematical models, such as the one proposed here, can help disentangle the multi-dimensional nature of cancer and lead to a better understanding of the drivers of treatment success and failure. While this work is only a first step, it shows that model-informed analysis of biomarkers and visible metastases during a single cycle of adaptive therapy may be useful in identifying important features of the metastatic population to provide further support for treatment decision making.

J. Zhang reports personal fees from AstraZeneca, Bayer, Sanofi, and Dendreon outside the submitted work. No disclosures were reported by the other authors.

J. Gallaher: Conceptualization, software, investigation, visualization, methodology, writing–original draft, writing–review and editing. M. Strobl: Conceptualization, methodology, writing–review and editing. J. West: Conceptualization, methodology, writing–review and editing. R. Gatenby: Data curation. J. Zhang: Data curation, writing–review and editing. M. Robertson-Tessi: Conceptualization, methodology, writing–review and editing. A.R.A. Anderson: Conceptualization, resources, supervision, funding acquisition, methodology, writing–review and editing.

The authors gratefully acknowledge funding by the NCI via the Cancer Systems Biology Consortium (CSBC) U01CA232382 (supporting J. Gallaher, M. Strobl, M. Robertson-Tessi, R. Gatenby, A.R.A. Anderson) and support from the Moffitt Center of Excellence for Evolutionary Therapy.

The publication costs of this article were defrayed in part by the payment of publication fees. Therefore, and solely to indicate this fact, this article is hereby marked “advertisement” in accordance with 18 USC section 1734.

Note: Supplementary data for this article are available at Cancer Research Online (http://cancerres.aacrjournals.org/).

1.
Dillekås
H
,
Rogers
MS
,
Straume
O
.
Are 90% of deaths from cancer caused by metastases?
Cancer Med
2019
;
8
:
5574
6
.
2.
Park
DS
,
Robertson-Tessi
M
,
Luddy
KA
,
Maini
PK
,
Bonsall
MB
,
Gatenby
RA
, et al
.
The goldilocks window of personalized chemotherapy: getting the immune response just right
.
Cancer Res
2019
;
79
:
5302
15
.
3.
Scarborough
JA
,
Tom
MC
,
Kattan
MW
,
Scott
JG
.
Revisiting a null hypothesis: exploring the parameters of oligometastasis treatment
.
Int J Radiat Oncol Biology Phys
2021
;
110
:
371
81
.
4.
Crawford
ED
,
Bennett
CL
,
Andriole
GL
,
Garnick
MB
,
Petrylak
DP
.
The utility of prostate-specific antigen in the management of advanced prostate cancer
.
Bju Int
2013
;
112
548
60
.
5.
Hussain
M
,
Tangen
CM
,
Higano
C
,
Schelhammer
PF
,
Faulkner
J
,
Crawford
ED
, et al
.
Absolute prostate-specific antigen value after androgen deprivation is a strong independent predictor of survival in new metastatic prostate cancer: data from Southwest Oncology Group Trial 9346 (INT-0162)
.
J Clin Oncol
2006
;
24
:
3984
90
.
6.
Enriquez-Navas
PM
,
Kam
Y
,
Das
T
,
Hassan
S
,
Silva
A
,
Foroutan
P
, et al
.
Exploiting evolutionary principles to prolong tumor control in preclinical models of breast cancer
.
Sci Transl Med
2016
;
8
:
327ra24–
.
7.
Gatenby
RA
,
Silva
AS
,
Gillies
RJ
,
Frieden
BR
.,
Adaptive therapy
.
Cancer Res
2009
;
69
:
4894
903
.
8.
Smalley
I
,
Kim
E
,
Li
J
,
Spence
P
,
Wyatt
CJ
,
Eroglu
Z
, et al
.
Leveraging transcriptional dynamics to improve BRAF inhibitor responses in melanoma
.
Ebiomedicine
2019
;
48
:
178
90
.
9.
Zhang
J
,
Cunningham
JJ
,
Brown
JS
,
Gatenby
RA
.
Integrating evolutionary dynamics into treatment of metastatic castrate-resistant prostate cancer
.
Nat Commun
2017
;
8
:
1816
.
10.
Kratiras
Z
,
Konstantinidis
C
,
Skriapas
K
.
A review of continuous vs intermittent androgen deprivation therapy: redefining the gold standard in the treatment of advanced prostate cancer. Myths, facts and new data on a ‘perpetual dispute’
.
Int Braz J Urol
2013
;
40
;
3
15
.
11.
Strum
SB
,
Scholz
MC
,
McDermed
JE
.
Intermittent androgen deprivation in prostate cancer patients: factors predictive of prolonged time off therapy
.
Oncol
2000
;
5
:
45
52
.
12.
Crook
JM
,
O'Callaghan
CJ
,
Duncan
G
,
Dearnaley
DP
,
Higano
CS
,
Horwitz
EM
, et al
.
Intermittent androgen suppression for rising PSA level after radiotherapy
.
Obstet Gynecol Surv
2013
;
68
:
34
35
.
13.
Magnan
S
,
Zarychanski
R
,
Pilote
L
,
Bernier
L
,
Shemilt
M
,
Vigneault
E
, et al
.
Intermittent vs continuous androgen deprivation therapy for prostate cancer: a systematic review and meta-analysis
.
Jama Oncol
2015
;
1
:
1
10
.
14.
West
J
,
Ma
Y
,
Newton
PK
.
Capitalizing on competition: an evolutionary model of competitive release in metastatic castration-resistant prostate cancer treatment
.
J Theor Biol
2018
;
455
:
249
60
.
15.
You
L
,
Brown
JS
,
Thuijsman
F
,
Cunningham
JJ
,
Gatenby
RA
,
Zhang
J
, et al
.
Spatial vs. nonspatial eco-evolutionary dynamics in a tumor growth model
.
J Theor Biol
2017
;
435
:
78
97
.
16.
Cunningham
JJ
,
Brown
JS
,
Gatenby
RA
,
Staňková
K
.
Optimal control to develop therapeutic strategies for metastatic castrate-resistant prostate cancer
.
J Theor Biol
2018
;
459
:
67
78
.
17.
West
JB
,
Dinh
MN
,
Brown
JS
,
Zhang
J
,
Anderson
AR
,
Gatenby
RA
.
Multidrug cancer therapy in metastatic castrate-resistant prostate cancer: an evolution-based strategy
.
Clin Cancer Res
2019
;
25
:
4413
21
.
18.
Strobl
MAR
,
Gallaher
J
,
West
J
,
Robertson-Tessi
M
,
Maini
PK
,
Anderson
ARA
.
Spatial structure impacts adaptive therapy by shaping intra-tumoral competition
.
Commun Medicine
2022
;
2
:
46
.
19.
Strobl
MAR
,
West
J
,
Viossat
Y
,
Damaghi
M
,
Robertson-Tessi
M
,
Brown
JS
, et al
.
Turnover modulates the need for a cost of resistance in adaptive therapy
.
Cancer Res
2021
;
81
:
1135
47
.
20.
Hansen
E
,
Woods
RJ
,
Read
AF
.
How to use a chemotherapeutic agent when resistance to it threatens the patient
.
Plos Biol
2017
;
15
:
e2001110
.
21.
Zhang
J
,
Cunningham
J
,
Brown
J
,
Gatenby
R
.
Evolution-based mathematical models significantly prolong response to abiraterone in metastatic castrate-resistant prostate cancer and identify strategies to further improve outcomes
.
Elife
2022
;
11
:
e76284
.
22.
Kim
E
,
Brown
JS
,
Eroglu
Z
,
Anderson
ARA
.
Adaptive therapy for metastatic melanoma: predictions from patient calibrated mathematical models
.
Cancers
2021
;
13
:
823
.
23.
Viossat
Y
,
Noble
R
.
A theoretical analysis of tumor containment
.
Nat Ecol Evol
2021
;
5
:
826
35
.
24.
Brady
R
, et al
.
Prostate-specific antigen dynamics predict individual responses to intermittent androgen deprivation
.
Biorxiv
2019
;
624866
.
25.
Hansen
E
,
Read
AF
.
Modifying adaptive therapy to enhance competitive suppression
.
Cancers
2020
;
12
:
3556
.
26.
Silva
AS
,
Kam
Y
,
Khin
ZP
,
Minton
SE
,
Gillies
RJ
,
Gatenby
RA
.
Evolutionary approaches to prolong progression-free survival in breast cancer
.
Cancer Res
2012
;
72
:
6362
70
.
27.
Ya. Tyuryumina
E
,
Neznanov
AA
.
Consolidated mathematical growth model of the primary tumor and secondary distant metastases of breast cancer (CoMPaS)
.
PLoS One
2018
;
13
:
e0200148
.
28.
Gallaher
J
,
Babu
A
,
Plevritis
S
,
Anderson
ARA
.
Bridging population and tissue scale tumor dynamics: a new paradigm for understanding differences in tumor growth and metastatic disease
.
Cancer Res
2014
;
74
:
426
35
.
29.
Avanzini
S
,
Antal
T
.
Cancer recurrence times from a branching process model
.
Plos Comput Biol
2019
;
15
:
e1007423
.
30.
Retsky
MW
,
Demicheli
R
,
Swartzendruber
DE
,
Bame
PD
,
Wardwell
RH
,
Bonadonna
G
, et al
.
Computer simulation of a breast cancer metastasis model
.
Breast Cancer Res Tr
1997
;
45
:
193
202
.
31.
Iwata
K
,
Kawasaki
K
,
Shigesada
N
.
A dynamical model for the growth and size distribution of multiple metastatic tumors
.
J Theor Biol
2000
;
203
:
177
86
.
32.
Benzekry
S
,
Hahnfeldt
P
.
Maximum tolerated dose versus metronomic scheduling in the treatment of metastatic cancers
.
J Theor Biol
2013
;
335
:
235
44
.
33.
Coumans
FA
,
Siesling
S
,
Terstappen
LW
.
Detection of cancer before distant metastasis
.
BMC Cancer
2013
;
13
:
283
.
34.
Franssen
LC
,
Lorenzi
T
,
Burgess
AEF
,
Chaplain
MAJ
.
A mathematical framework for modelling the metastatic spread of cancer
.
B Math Biol
2019
;
81
:
1965
2010
.
35.
Szczurek
E
,
Krüger
T
,
Klink
B
,
Beerenwinkel
N
.
A mathematical model of the metastatic bottleneck predicts patient outcome and response to cancer treatment
.
Plos Comput Biol
2020
;
16
:
e1008056
.
36.
Heyde
A
,
Reiter
JG
,
Naxerova
K
,
Nowak
MA
.
Consecutive seeding and transfer of genetic diversity in metastasis
.
Proc Natl Acad Sci USA
2019
;
116
:
14129
37
.
37.
Liotta
LA
,
Saidel
GM
,
Kleinerman
J
.
Stochastic model of metastases formation
.
Biometrics
1976
;
32
:
535
50
.
38.
Liotta
LA
,
Delisi
C
,
Saidel
G
,
Kleinerman
J
.
Micrometastases formation: a probabilistic model
.
Cancer Lett
1977
;
3
:
203
8
.
39.
Gerlee
P
,
Johansson
M
.
Inferring rates of metastatic dissemination using stochastic network models
.
Plos Comput Biol
2019
;
15
:
e1006868
.
40.
Rhodes
A
,
Hillen
T
.
A mathematical model for the immune-mediated theory of metastasis
J Theor Biol
2019
;
482
:
109999
.
41.
Scott
JG
,
Basanta
D
,
Anderson
ARA
,
Gerlee
P
.
A mathematical model of tumor self-seeding reveals secondary metastatic deposits as drivers of primary tumor growth
.
J Roy Soc Interface
2013
;
10
:
20130011
.
42.
Zhang
J
,
Cunningham
JJ
,
Brown
JS
,
Gatenby
RA
.
Response to Mistry
.
Nat Commun
2021
;
12
:
329
.
43.
Gallaher
JA
,
Enriquez-Navas
PM
,
Luddy
KA
,
Gatenby
RA
,
Anderson
ARA
.
Spatial heterogeneity and evolutionary dynamicsmodulate time to recurrence in continuous andadaptive cancer therapies
.
Cancer Res
2018
;
78
:
2127
39
.
44.
Forster
JC
,
Harriss-Phillips
WM
,
Douglass
MJ
,
Bezak
E
.
A review of the development of tumor vasculature and its effects on the tumor microenvironment
.
Adv Exp Med Biol
2017
;
5
:
21
32
.
45.
Pierorazio
PM
,
Walsh
PC
,
Partin
AW
,
Epstein
JI
.
Prognostic Gleason grade grouping: data based on the modified Gleason scoring system
.
Bju Int
2013
;
111
:
753
60
.
46.
Shih
WJ
,
Ohman-Strickland
PA
,
Lin
Y
.
Analysis of pilot and early phase studies with small sample sizes
.
Statist. Med
2004
;
23
:
1827
42
.
47.
Bruchovsky
N
,
Klotz
L
,
Crook
J
,
Malone
S
,
Ludgate
C
,
Morris
WJ
, et al
.
Final results of the Canadian prospective phase II trial of intermittent androgen suppression for men in biochemical recurrence after radiotherapy for locally advanced prostate cancer
.
Cancer
2006
;
107
:
389
95
.
48.
Bruchovsky
N
,
Klotz
L
,
Crook
J
,
Goldenberg
SL
.
Locally advanced prostate cancer—biochemical results from a prospective phase II study of intermittent androgen suppression for men with evidence of prostate-specific antigen recurrence after radiotherapy
.
Cancer
2007
;
109
:
858
67
.
49.
Gallaher
JA
,
Brown
JS
,
Anderson
ARA
.
The impact of proliferation-migration trade-offs on phenotypic evolution in cancer
.
Sci Rep
2019
;
9
:
2425
.
50.
Lenormand
T
,
Harmand
N
,
Gallet
R
.
Cost of resistance: an unreasonably expensive concept
.
Biorxiv
2018
;
276675
.
51.
Jessup
CM
,
Bohannan
BJM
.
The shape of an ecological trade-off varies with environment
.
Ecol Lett
2008
;
11
:
947
59
.
52.
Frankenstein
Z
,
Basanta
D
,
Franco
OE
,
Gao
Y
,
Javier
RA
,
Strand
DW
, et al
.
Stromal reactivity differentially drives tumor cell evolution and prostate cancer progression
.
Nat Ecol Evol
2020
;
4
:
870
84
.
53.
Brady
L
,
Kriner
M
,
Coleman
I
,
Morrissey
C
,
Roudier
M
,
True
LD
, et al
.
Inter- and intra-tumor heterogeneity of metastatic prostate cancer determined by digital spatial gene expression profiling
.
Nat Commun
2021
;
12
:
1426
.
54.
Li
Q
,
Deng
Q
,
Chao
H-P
,
Liu
X
,
Lu
Y
,
Lin
K
, et al
.
Linking prostate cancer cell AR heterogeneity to distinct castration and enzalutamide responses
.
Nat Commun
2018
;
9
:
3600
.
55.
Morin
F
,
Beauregard
J-M
,
Bergeron
M
,
Nguile Makao
M
,
Lacombe
L
,
Fradet
V
, et al
.
Metabolic imaging of prostate cancer reveals intrapatient inter-metastasis response heterogeneity to systemic therapy
.
European Urology Focus
2017
;
3
:
639
42
.
56.
Logothetis
CJ
,
Gallick
GE
,
Maity
SN
,
Kim
J
,
Aparicio
A
,
Efstathiou
E
, et al
.
Molecular classification of prostate cancer progression: foundation for marker-driven treatment of prostate cancer
.
Cancer Discov
2013
;
3
:
849
61
.
57.
Labrecque
MP
,
Coleman
IM
,
Brown
LG
,
True
LD
,
Kollath
L
,
Lakely
B
, et al
.
Molecular profiling stratifies diverse phenotypes of treatment-refractory metastatic castration-resistant prostate cancer
.
J Clin Invest
2019
;
129
:
4492
505
.
58.
Hunter
KW
,
Amin
R
,
Deasy
S
,
Ha
N-H
,
Wakefield
L
.
Genetic insights into the morass of metastatic heterogeneity
.
Nat Rev Cancer
2018
;
18
:
211
23
.
59.
Klein
CA
.
Parallel progression of primary tumors and metastases
Nat Rev Cancer
2009
;
9
:
302
12
.
60.
Zhao
Z-M
,
Zhao
B
,
Bai
Y
,
Iamarino
A
,
Gaffney
SG
,
Schlessinger
J
, et al
.
Early and multiple origins of metastatic lineages within primary tumors
.
Proc Natl Acad Sci USA
2016
;
113
:
2140
5
.
61.
Gundem
G
,
Van Loo
P
,
Kremeyer
B
,
Alexandrov
LB
,
Tubio
JMC
,
Papaemmanuil
E
, et al
.
The evolutionary history of lethal metastatic prostate cancer
.
Nature
2015
;
520
:
353
7
.
62.
Araujo
A
,
Cook
LM
,
Lynch
CC
,
Basanta
D
.
Size matters: metastatic cluster size and stromal recruitment in the establishment of successful prostate cancer to bone metastases
.
B Math Biol
2018
;
80
:
1046
58
.
63.
Hu
Z
,
Li
Z
,
Ma
Z
,
Curtis
C
.
Multi-cancer analysis of clonality and the timing of systemic spread in paired primary tumors and metastases
.
Nat Genet
2020
;
52
:
701
8
.
64.
Klein
CA
,
Hölzel
D
.
Systemic cancer progression and tumor dormancy: mathematical models meet single cell genomics
.
Cell Cycle
2006
;
5
:
1788
98
.
65.
Willis
L
,
Alarcón
T
,
Elia
G
,
Jones
JL
,
Wright
NA
,
Tomlinson
IPM
, et al
.
Breast cancer dormancy can be maintained by small numbers of micrometastases
.
Cancer Res
2010
;
70
:
4310
7
.
66.
Midde
K
,
Sun
N
,
Rohena
C
,
Joosen
L
,
Dhillon
H
,
Ghosh
P
.
Single-cell imaging of metastatic potential of cancer cells
.
Iscience
2018
;
10
:
53
65
.
67.
Merino
D
,
Weber
TS
,
Serrano
A
,
Vaillant
F
,
Liu
K
,
Pal
B
, et al
.
Barcoding reveals complex clonal behavior in patient-derived xenografts of metastatic triple-negative breast cancer
.
Nat Commun
2019
;
10
:
766
.
68.
Pérez-Velázquez
J
,
Rejniak
KA
.
Drug-induced resistance in micrometastases: analysis of spatio-temporal cell lineages
.
Front Physiol
2020
;
11
:
319
.
69.
Sleeman
JP
,
Christofori
G
,
Fodde
R
,
Collard
JG
,
Berx
G
,
Decraene
C
, et al
.
Concepts of metastasis in flux: the stromal progression model
.
Semin Cancer Biol
2012
;
22
:
174
86
.
This open access article is distributed under the Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International (CC BY-NC-ND 4.0) license.