## Abstract

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.

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

## Introduction

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.

## Materials and Methods

### 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.

#### 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 *N _{Dx}*. 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

*T*; 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

_{seed}*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

*P*is the maximum probability of cell death per day, and

_{d,max}*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:

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.

## Results

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).

### 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. 2D–F), 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. 2G–I) on average has shorter response times and longer regrowth times (3.5 ± 2.1 and 14.0 ± 5.9 weeks, respectively).

Parameter . | Description . | Range . |
---|---|---|

N _{Dx} | Total number of cells at diagnosis | 10 × 10^{3} to 150 × 10^{3} |

M | Number of metastases seeded | 1–10 |

r _{0} | Initial mean drug resistance | 0–0.4 |

d _{r} | Random death rate | 0 × 10^{−2} to 6 × 10^{−2}/day |

h _{i} | Intermetastasis heterogeneity | 0–1 |

h _{a} | Intrametastasis heterogeneity | 0–1 |

Parameter . | Description . | Range . |
---|---|---|

N _{Dx} | Total number of cells at diagnosis | 10 × 10^{3} to 150 × 10^{3} |

M | Number of metastases seeded | 1–10 |

r _{0} | Initial mean drug resistance | 0–0.4 |

d _{r} | Random death rate | 0 × 10^{−2} to 6 × 10^{−2}/day |

h _{i} | Intermetastasis heterogeneity | 0–1 |

h _{a} | 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.

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.

### 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.

Intriguingly, Fig. 5D–F 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. 5D–F). 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. 5G–I). 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.

### 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).

Next, we wanted to assign a Gleason score to the virtual patients in the cohort of 1,000 simulations used in Fig. 2C–F (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. 7A–C, 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. 7D–F. 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.

## Discussion

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.

## Authors' Disclosures

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

## Authors' Contributions

**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.

## Acknowledgments

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/).