Spatial structure impacts adaptive therapy by shaping intra-tumoral competition.

Background: Adaptive therapy aims to tackle cancer drug resistance by leveraging resource competition between drug-sensitive and resistant cells. Here, we present a theoretical study of intra-tumoral competition during adaptive therapy, to investigate under which circumstances it will be superior to aggressive treatment.
Methods: We develop and analyse a simple, 2-D, on-lattice, agent-based tumour model in which cells are classified as fully drug-sensitive or resistant. Subsequently, we compare this model to its corresponding non-spatial ordinary differential equation model, and fit it to longitudinal prostate-specific antigen data from 65 prostate cancer patients undergoing intermittent androgen deprivation therapy following biochemical recurrence.
Results: Leveraging the individual-based nature of our model, we explicitly demonstrate competitive suppression of resistance during adaptive therapy, and examine how different factors, such as the initial resistance fraction or resistance costs, alter competition. This not only corroborates our theoretical understanding of adaptive therapy, but also reveals that competition of resistant cells with each other may play a more important role in adaptive therapy in solid tumours than was previously thought. To conclude, we present two case studies, which demonstrate the implications of our work for: (i) mathematical modelling of adaptive therapy, and (ii) the intra-tumoral dynamics in prostate cancer patients during intermittent androgen deprivation treatment, a precursor of adaptive therapy.
Conclusion: Our work shows that the tumour's spatial architecture is an important factor in adaptive therapy and provides insights into how adaptive therapy leverages both inter- and intra-specific competition to control resistance.

Introduction

“Can insects become resistant to sprays?”—this is the question entomologist Axel Melander raised in an article of the same title in 1914[1]. At a site in Clarkston, WA, Melander had observed that over 90% of an insect pest called the “San Jose scale” was surviving despite being sprayed with sulphur-lime insecticide[1]. If we were to ask the same question in cancer treatment today we would be met with an equally resounding “yes”: while for most cancers it is possible to achieve an initial, possibly significant, burden reduction, many patients recur with drug-resistant disease, or even progress while still under treatment. Drug resistance can develop in a number of ways, including genetic mutations that alter drug binding, changes in gene expression, which activate alternative signalling pathways, or environmentally mediated resistance[2-4]. In the clinic, the main strategy for managing cancer drug resistance is to switch treatment with the aim of finding an agent to which the tumour is still susceptible[3,4]. Similarly, Melander suggested it might be possible to tackle sulphur-lime resistance by switching to oil-based sprays[1]. However, he also foresaw the possibility of, and the challenges arising from, multi-drug resistance[1]. As an alternative, Melander proposed that it might be possible to maintain insecticide sensitivity through less aggressive spraying as this would promote inter-breeding of sensitive and resistant populations and, thus, dilute the resistant genotype[1]. Allocation of spray-free “refuge” patches in the neighbourhood of plots in which an insecticide is used is one modality of modern pest management and is even required by law for the use of certain agents in the US (e.g., Bt-crops[5]). Similarly, research into antibiotic resistance has investigated strategic modulation and combination of treatments to suppress, and ideally reverse, resistance evolution (see Baym et al.[6] for a comprehensive review). For example, Abel zur Wiesch et al.[7] found in a meta-analysis that “adjusted cycling” of drugs, where treatments were switched when resistance was detected, reduced the evolution of antibiotic resistance in hospital wards. Moreover, Hansen et al.[8] have shown that by maintaining drug-sensitive bacteria they can slow the emergence of resistant cells in a bioreactor.

Recently, the concept that treatment de-escalation can delay the emergence of resistance has found application also in oncology. Standard-of-care cancer treatment regimens aim to maximise cell kill through application of the maximum tolerated dose (MTD), in order to achieve a cure. In contrast, an emerging approach called adaptive therapy proposes to focus not on burden reduction, but on burden control in settings, such as advanced, metastatic disease, in which cures are unlikely[9-12]. Eradication strategies free surviving cells from intra-tumoral resource competition, which would otherwise inhibit resistance growth. Adaptive therapy aims to leverage this competition by maintaining drug-sensitive cells in order to avoid, or at least delay, the emergence of resistance[9,11]. A number of pre-clinical studies have demonstrated the feasibility of this approach in ovarian[10], breast[13], colorectal[14], and skin[15] cancer. While large-scale, randomised clinical trials are outstanding, a pilot trial of adaptive therapy in metastatic castration-resistant prostate cancer achieved not only an at least 10 month increase in median time to progression (TTP), but also a 53% reduction in cumulative drug usage in comparison to a contemporaneous control cohort[16]. Further clinical trials in castration-sensitive prostate cancer and melanoma are ongoing (clinicaltrials.gov identifiers NCT03511196 and NCT03543969, respectively).

In addition to testing its feasibility, there has been significant interest in characterising the underpinning eco-evolutionary principles of adaptive therapy through mathematical modelling. We identify three key results. The first insight was derived from approaches which represent the tumour as a mixture of drug-sensitive and resistant cells modelled as a system of two or more ordinary differential equations (ODEs) with competition described by the Lotka–Volterra model from ecology[10,14,15,17-20] or by a matrix game[21]. These analyses have demonstrated that less aggressive treatment allows for longer tumour control under a range of assumptions on the tumour growth law (exponential:[14,17,19,21]; logistic:[10,15,17,19]; Gompertzian:[17-19]; dynamic carrying capacity:[14,20]), and the origin of resistance (pre-existing:[10,14,17-19,21]; acquired[15,18,19,22]; cancer stem-cell-based:[20]). Furthermore, this work predicts that adaptive therapy will be most effective in cases where cures are unlikely due to pre-existing resistance and where at the same time conditions (resistance fraction, proximity to carrying capacity) are such that inter-specific competition with drug-sensitive cells is strong (see ref. [19] for a comprehensive and formal summary of these results). The second key result is that while these conclusions broadly transfer to more complex, spatially explicit tumour models, the strength of spatial constraints on resistant cell growth is important[14,23]. Bacevic et al.[14] showed in a two-dimensional (2-D), on-lattice, agent-based model (ABM) of a tumour spheroid that longer control is achieved if resistance arises in the centre of the tumour compared to when it arises on the edge. Gallaher et al.[23] corroborated this result in a 2-D, off-lattice setting with resistance modelled as a continuum, and further demonstrated that tumour control was adversely affected by high cell motility and cell plasticity. Thirdly, models focussed on metastatic prostate cancer have illustrated how these concepts may be realised in a specific disease pathology[16,20] and how we may enhance tumour control by using a multi-drug approach[20,24,25].

But, what does the competitive landscape of a resistant cell actually look like? With whom do resistant cells compete and at what rate? Even though competition is a key ingredient of adaptive therapy, to the best of our knowledge, no study to-date has explicitly quantified it. Moreover, non-spatial work typically models competition phenomenologically, using Lotka–Volterra dynamics (or a Gompertzian analogue[18,19]). However, this assumes perfect mixing of cells, and is likely an inaccurate description of the dynamics in solid tumours. We, and others[26-29], have recently shown that spatial constraints may alter the nature of tumour evolution away from what may be expected from non-spatial models, an observation which has also been made in the antibiotic resistance community[30,31]. Better understanding the impact of space on the ecology and evolution during treatment is therefore important, and may help to more accurately identify for whom and how to adapt treatment in oncology, and beyond.

The goal of this paper is to study competitive suppression of resistant cells during adaptive therapy and how this is modulated by space. To do so, we developed a simple, 2-D, on-lattice, ABM in which the tumour is assumed to be composed of two cell types: drug-sensitive and resistant cells. The individual-based, spatially explicit nature of our model allowed us to directly quantify competition for space. Leveraging this, we show, for the first time, the competitive suppression of resistant cells during adaptive therapy, and we discuss how the initial proximity to carrying capacity, the initial resistance fraction, the presence of resistance costs, and the rate of cell turnover alter competition. This analysis not only provides useful validation of current theory, but also highlights a seldomly discussed factor: namely, that resistant cells compete not only with sensitive cells, but also with each other. We show that this observation is important in solid tumours, where mixing of cells is limited, because it implies that the spatial distribution of resistant cells strongly impacts treatment response under adaptive therapy. Subsequently, we discuss the implications of our new insights for the modelling of adaptive therapy using ODE models, by comparing our ABM to its equivalent mean-field ODE approximation, which we studied recently[32]. To conclude, we present an analysis in which we fitted our ABM to publicly available, longitudinal data from 65 prostate cancer patients undergoing intermittent androgen deprivation therapy. A striking feature of these data is that patients cycle between on- and off-treatment phases at different frequencies. Based on our model fits, we propose that the variations in cycling speed reflect differences in the spatial organisation of these tumours, with implications for the balance of intra- and inter-specific competition between sensitive and resistant cells within. Overall, our work helps to provide a more detailed understanding of spatial competition between sensitive and resistant cells during adaptive therapy and shows that the spatial architecture of the tumour can strongly affect treatment outcomes. While we focus here on cancer, we believe that because of the parallels with strategies explored in other areas, such as antibiotic resistance, our insights may be also of interest to the wider scientific community.

Methods

The mathematical model

Random geno- and phenotypic variation produces tumour cells, which show a degree of drug-resistance even prior to drug exposure. This may manifest as an increased ability to persist and adapt to adverse conditions such as drug exposure, or, though perhaps more rarely, it may take the form of fully developed resistance[3,4]. Selective expansion and further adaptation of this population is thought to be the cause of treatment failure in patients[33,34].

To study the evolutionary dynamics in response to treatment we consider a 2-D, on-lattice, ABM representative of a small region of tumour tissue or a metastatic site. For simplicity we assume that we can divide cells into drug-sensitive or fully drug-resistant subpopulations (Fig. 1a). We choose an on-lattice, agent-based representation as it allows us to explore the role of space and cell-scale stochasticity in a tractable, yet generalisable, way. Each cell occupies a single site in an l × l square lattice with no-flux boundary conditions, and behaves according to the following rules (Fig. 1):We denote the number of cells in each population at time t by S(t) and R(t), and the total number by N(t) = S(t) + R(t), respectively (Table 1).

Fig. 1

The agent-based tumour model.

a The tumour is modelled as a mixture of drug-sensitive (S) and resistant cells (R), where each cell occupies a square on a 2-D equi-spaced lattice. Cells divide and die at constant rates r and r, and d and d, respectively. Daughter cells are placed into empty squares in a cell’s von Neumann neighbourhood. Drug (D) will kill dividing sensitive cells at a probability d. b Flow diagram of the simulation algorithm (tEnd: end time of the simulation; dt: simulation time step). For parameter values, see Table 1.

Table 1

Summary of mathematical model variables and parameters.

Parameter	Description	Value	Comment
l	Grid size (Total number of sites: l × l)	100
dt	Simulation time step	1d
S(t)	Number of sensitive cells	0–l2
R(t)	Number of resistant cells	0–l2
N(t)	Total tumour cell number	0–l2
rS	Sensitive cell proliferation rate	0.027 d−1	Adopted from ref. [16].
cR	Resistance cost (resistant proliferation rate: rR = (1 − cR)rS)	0–50%	Lower limit:adopted from ref. [23]; Upper limit: assumption of no cost
dT	Cell death rate (relative to rS)	0–50%	Lower limit: assumption of no turnover;
			Upper limit: see discussion in ref. [32].
dD	Drug-induced cell kill probability of sensitive cells at D(t)=D_Max	0.75	Adopted from ref. [24]; verified to provide a good fit to prostate cancer in ref. [32].
n_0=N_0l^2	Initial cell density (as a percentage of total carrying capacity)	25–75%	Values within this range reported by ref. [77].
f_R:=R_0N_0	Initial resistant cell fraction (as a percentage of initial cell density)	0.1–10%	Values within this range reported by ref. [78].
f_S:=S_0N_0	Initial sensitive cell fraction (as a percentage of initial cell density)	90–99.9%	Determined by 1 − fR.

Initially, there are a total of N0 cancer cells in the tissue of which a fraction f is resistant. Generally, we will assume that the cells are spread randomly throughout the tissue, except in section “The spatial distribution of resistance impacts adaptive therapy by shaping intra-tumoral competition”  where we will explore a clustered and a disk-like configuration to dissect the role of the initial conditions more explicitly.

Sensitive and resistant cells attempt to divide at constant rates r and r (in units: d−1), respectively. If there is at least one empty site in the cell’s von Neumann neighbourhood (consisting of the four lattice neighbours, east, west, north, and south, of the cell), then the cell will divide and the daughter will be placed randomly in one of the empty sites in the neighbourhood.

Cells die at a constant rate δ (in units: d−1). For notational convenience, we will express this rate relative to the sensitive cell proliferation rate, so that d = δ/r. Note that this definition of turnover compares the cell death rate to the cells’ intrinsic proliferation rate, and is thus not the same as the sometimes measured “cell loss” rate (refs. [35,36]; see also ref. [32] for a further comparison of the two). In addition, we will make the simplifying assumption that that both sensitive and resistant cells die at the same rate, d = d = d.

Movement of cells is neglected.

The domain is sufficiently small so that drug concentration is assumed to be spatially homogeneous throughout the tissue, where is the MTD.

A sensitive cell which is currently undergoing mitosis—that is, it has attempted division and has space available in its neighbourhood—is killed by drug with probability dD(t), where .

Dead cells are immediately removed from the domain.

We consider two treatment schedules:

Continuous therapy at MTD: .

Adaptive therapy as implemented in the Zhang et al.[16] prostate cancer clinical trial: Treatment is withdrawn once a 50% decrease from the initial tumour size is achieved, and is reinstated if the original tumour size (N0) is reached:

This results in cycles of on- and off-treatment periods, which maintain the tumour burden at at least 50% its original level for as long as possible, and thereby seek to slow the expansion of resistance.

Progression was determined as a 20% increase from the pre-treatment baseline. However, occasionally it could happen that during adaptive therapy the tumour burden briefly exceeded this target at the end of the first or second off-cycle, due to rapid regrowth of the sensitive cells. As the tumour in these cases was immediately brought back under control upon treatment re-administration, and the focus of our study was progression driven by drug-resistant cells, we neglected these events and only considered a tumour to have progressed if at least 150 days had passed since the start of treatment. This criterion was applied to both adaptive and continuous therapy.

A flow-chart of our model is shown in Fig. 1b, and further implementation details are given in Supplementary Methods 1. We checked convergence (not shown), and performed a consistency analysis[37,38]. This showed that a sample size, n, upward of 250 provides a representative sample size for our stochastic simulation algorithm (see Supplementary Fig. 1 and Supplementary Methods 2 for details). The model is implemented in Java 1.8. using the Hybrid Automata Library[39]. Data analysis was carried out in Python 3.6, using Pandas 0.23.4, Matplotlib 2.2.3, Seaborn 0.9.0, and openCV 3.4.9. The time-evolution of the resistant cells’ neighbourhood composition was visualised using EvoFreq[40] in R 4.0.2. All code is available on GitHub at https://github.com/MathOnco/strobl2021_space_modulates_competition_AT (see also ref. [41]).

Comparison with the non-spatial model

To understand the impact of space we compared the ABM to the following ODE model, which we have studied previously in ref. [32]:where K is the carrying capacity, and the initial conditions are given by S(0) = S0, R(0) = R0, and N0 = S0 + R0, respectively. We set K = l2 and used the same parameter values as for the agent-based simulation otherwise (Table 1). The equations were solved using the RK45 (used when comparing the ABM and ODE model in section “How competition is modelled matters”) or DOP853 (used for faster computational performance when fitting the patient data in section “The cycling frequency of patients undergoing intermittent androgen deprivation therapy may reflect different spatial distributions of resistance”) explicit Runge–Kutta schemes provided in Scipy (for further details see ref. [32]).

Model parameters

We parametrised our model using values from the literature for prostate cancer (Table 1). We want to stress, however, that the aim of our work was to develop qualitative understanding, not to make quantitative predictions directed at prostate cancer. As such, our predictions should be interpreted not in a quantitative (“treatment X will achieve a TTP of Y months”), but in a qualitative fashion (“treatment X will achieve a longer TTP than treatment Y because of mechanism Z”).

Analysis of patient data

In order to examine whether our model could explain differences in the cycling speed of patients undergoing intermittent androgen deprivation therapy, we fitted it to the publicly available, longitudinal response data from the Phase II trial by Bruchovsky et al.[42]. The data were downloaded from http://www.nicholasbruchovsky.com/clinicalResearch.html in July 2020. So as to avoid potentially confounding effects from a change in the number of lesions, patients who developed a metastasis were excluded from analysis. Furthermore, we decided to exclude two further patients (Patients 2 & 104) from our analysis because their prostate-specific antigen (PSA) dynamics were inconsistent with the reported treatment schedules. These patients display oscillating PSA values indicative of treatment cycling but there are no reported changes in treatment, suggesting there may have been a mistake in the data reporting in these cases (see also Supplementary Fig. 2). This yielded data from a total of 65 patients. The model was fitted to the normalised PSA measurements by minimising the root mean-squared error between the data and the predictions (normalisation relative to PSA level at start of treatment). Given the stochastic nature of the ABM, each candidate fit was assembled from 25 independent stochastic replicates. Optimisation was carried out using the basin-hopping algorithm in Scipy[43] employing default search parameters and a maximum of either 50 (when fitting 2 parameters) or 75 optimisation steps (when fitting 4 parameters). To escape potential local minima, optimisation was repeated 10 times for each patient from different, randomly chosen initial conditions, with only the best fit according to the Akaike Information Criterion (AIC) taken forward for analysis. We fitted the model in three different ways: (i) each of n0, f, c, and d was allowed to be a patient-specific parameter (which we term “Model 1”), (ii) only n0 and f were allowed to be patient-specific with c and d fixed at the mean values obtained in (i) (c = 0.78, d = 0.14; Model 2), and (iii) c and d were allowed to be patient-specific with n0 and f fixed at the mean values obtained in (i) (n0 = 0.59, f = 0.04; Model 3). For 8 patients, Model 1 failed to recapitulate the cycling nature of the patients’ trajectories, yielding simple straight lines instead (e.g., Patient 52 in Supplementary Fig. 2). Owing to this discrepancy it was unclear that the associated parameter values would be representative of the tumour biology, and so we excluded these patients when computing the mean of the parameters taken forward to Models 2 and 3. For the same reasons, we excluded these 8 patients and one additional patient (Patient 13) when assessing the correlation between cycling speed, and the cost and turnover estimates of Model 3. This analysis was thus based on the data from 56 patients in total. Excluded patients are marked by a grey background in Supplementary Figs. 2 & 3, respectively. Classification of patients into “progressing” and “non-progressing” was taken from the annotation provided in the data, where progression is defined as a series of three sequential increases of serum PSA > 4.0 μg/L despite castrate levels of serum testosterone. Overviews of all fits of Models 1 and 3 are shown in Supplementary Figs. 2 & 3, respectively. Fitting was done using the lmfit package in Python[44] (version 1.0.1.). As this was a retrospective analysis of a study that was previously approved by the institutional review boards/independent ethics committees of each participating site, and whose data were available in the public domain, no further ethical approval or informed consent was required.