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.

(b) Treatment withdrawal after a 30% tumour burden reduction. This suggests that less aggressive treatment, especially if resistance is prevalent prior to treatment, will result in longer tumour control, and matches the results obtained for non-spatial models [2,3,4,5]. Note though that we are not taking into account possible risks associated with the higher tumour burden under the 30% threshold algorithm. Parameters: (n 0 , f R , c R , d T ) = (75%, 1%, 0%, 0%); n = 250 independent replicates.
Supplementary Discussion 1: Comparison of random initial seeding with seeding as a circular domain Throughout the main paper we initialise our simulations by seeding cells randomly in the domain. This was meant to recapitulate the heterogeneous structure of tumours in which nests of tumour cells are interspersed with remnants of normal tissue, tumour stroma, or areas of necrosis. In addition, it means that resistant cells are not guaranteed to be surrounded by sensitive cells, so that this appeared to be a good "worst case" scenario in which to test adaptive therapy. However, a more common way in which tumours are modelled in the literature is as discs (in 2-D) or balls of cells (in 3-D), in which the tumour grows radially outwards from its site of origin and the majority of cell division takes place on the surface of the growing tumour (e.g. [6,7]).
To examine how our conclusions are affected by such spherical tumour architecture, we repeated some of our analyses with the cells seeded as a disc in the centre of the domain, surrounded by empty space. To allow comparison with the random initial seeding in the main paper, we chose the radius of this disk to be given by: so that approximately the same total number of cells was seeded. We then seeded exactly the same number of resistant cells as in the random case within this disk (R 0 = f R n 0 10 4 ), and filled up the rest of the disk with sensitive cells. Furthermore, in order to allow cells on the surface to expand freely, we increased the domain size to l = 150px. Supplementary Figure 6a shows a comparison of the relative time gained by adaptive therapy over continuous therapy for the two types of initial conditions. We find that which of the two is predicted to benefit more from adaptive therapy depends on the initial cell density (n 0 ), the initial resistant cell fraction (f R ), and the cell turnover (d T ) in a multi-factorial fashion. Supplementary Figure 6b shows the treatment response when we assume the same parameters as in However, while seeding the cells in this way means that clones on the outside can expand more freely, it also means that cells are in stronger competition with each other at the beginning of the simulation (compare the right panel in Supplementary Figure 6b with Figure 2c). As a result, when we start with a small number of cells (small n 0 ) adaptive therapy is more effective when cells are seeded in a circular fashion than when they are placed randomly (Supplementary Figures 6a & c).
A further implication of this increased competition is that turnover has a greater beneficial effect on adaptive therapy when cells are seeded as a disk than when they are seeded randomly ( Supplementary Figures 6a & d). This is because turnover means that drug can kill cells throughout the tumour, so that not as much of the cell mass near the surface needs to be removed to achieve the 50% burden reduction required for treatment withdrawal (compare tumour sizes at t = 500d in centre panels in Supplementary Figures 6b & d). Consequently, resistant nests near the edge remain surrounded by sensitive cells for longer and can, thus, be controlled better than in the random seeding (or no turnover) case.
To sum up, these results corroborate our findings in the paper that the tumour's spatial organisation determines whether, and for how long, resistance may be controlled by adaptive therapy. This is because the architecture shapes the nature of intra-tumoral competition. We advocate that future research should study the impact of different tissue structures more explicitly, using for example the different "onco-evotypes" recently defined by Noble et al [8]. Importantly, such research should also interrogate the role of non-tumour tissue which we have neglected here. First steps in this direction have recently been presented by M A et al [9].  Figure 2a-c, but with cells seeded as a disk rather than randomly. Adaptive therapy is unable to control resistance as treatment frees resistant colonies near the edge which can subsequently expand in an unhindered way. (c) Simulations with the same parameters as (b) but with a smaller initial cell number. The smaller cell number means that the tumour radius has to be reduced by less to hit the 50% threshold at which treatment is withdrawn so that it takes longer for resistant cells to be freed (compare tumour size change between t = 0d and t = 500d in the middle panels in (b) and (c)). (d) Simulations with the same parameters as (b) but with a 30% turnover rate. Turnover means that drug kill is spread more evenly throughout the tumour so that resistant cells near the edge are contained for longer than in the absence of turnover. This corroborates our main finding that the tumour architecture impacts adaptive therapy by changing the competition dynamics of both resistant and sensitive cells. Panels (b) -(d) are based on n = 250 independent replicates.

S(t) R(t) N(t) under AT Drug on TTPCT TTPAT
Supplementary Figure 7: Comparison of the treatment predictions of the ABM and ODE model (Equations (2)-(4)) for the same initial density (n 0 = 75%) but different initial resistance fractions. We assume no cost or turnover. For the ABM the mean and standard deviation of n = 250 independent replicates are shown. When the initial resistance fraction is small, the ABM predicts significantly longer TTP than the ODE model. When the resistance fraction is large, the converse is true. This indicates that because of the impact of space, different initial numbers of resistant cells (and thus independent colonies) result in distinct progression dynamics in the ABM (see also Figure  5).

S(t) R(t) N(t) under AT
. We previously showed that K Eff is an important factor in determining the benefit of adaptive therapy [5]. As such, the differences in K Eff explain the discrepancy in the predicted benefit of adaptive therapy. (e) The reason why the effective carrying capacities differ is the fact that in the ABM a cell has four potential neighbouring sites into which it can divide. This allows for more division to take place than predicted by the ODE model. To show this, we compare the effective carrying capacity in the ABM for different values of cost and turnover to an expression derived by Kimmel et al [10] which accounts for the neighbourhood size:  [11]; negative differences indicate that continuously treated colonies are rounder than adaptively treated ones) at 250d ((n 0 , f R , c R , d T ) = (25%, 0.1%, 0%, 0%); n = 1000 independent replicates.). The observed correlation (Pearson's correlation coefficient, r = 0.41, p-value < 0.01) indicates that instances in which nests under adaptive therapy are more branched than nests under continuous therapy (smaller convexity values) are associated with poorer performance of AT. To compute a nest's convexity we used first the connectedComponentsWithStats() function in openCV [12] to identify individual nests, followed by the convexHull() function to find their convex hull. Supplementary Figure 12: Fitting the patient data assuming a disk-like tumour growth model, where the tumour begins as a disk in the centre of the domain and expands radially outwards. Tumours were initialised as described in Supplementary Discussion 1 (except that resistant cell numbers were not matched, so that f R here represents the true initial resistant fraction) and fitted using the same methods as described in Section 2.4 of the main text, allowing n 0 , f R , c R , and d T to vary on a patientspecific basis. Shown are the results for the subset of 11 patients whose dynamics are displayed in Supplementary Video 4. (a) Plot of the inferred cost (c R ) and turnover (d T ) values, revealing a correlation between the two values also in this case (Pearson's correlation coefficient: r = −0.87, p < 0.01). (b) Examples of treatment dynamics predicted by the model for three patients (n = 250 independent replicates in each case). This shows that, as in the case of random initial seeding, faster cycling patients are associated with more diffuse ("carpet-like") intra-tumoral architecture, whereas patients who are cycling slowly are fitted by simulations in which the dynamics is dominated by few, large sensitive colonies ("patch-like"). Supplementary Figure 13: Comparison of the patient fits for the ODE and the ABM model (varying n 0 , f R , c R , d T ). ODE fits taken from our prior work [5]. ABM simulations show mean and standard deviation from n = 25 independent replicates. (a) Comparison of the r 2 values of both models, showing that the goodness-of-fit of both models is generally comparable, although the ODE achieves slightly better fits on average (each point is one patient; n = 65). This is not unexpected given the extra complexities involved in fitting the ABM, such as stochasticity, longer run times, and gradientbased vs heuristic optimisation. Nevertheless, we believe that our new ABM results provide valuable insights, as the ODE model allowed us to formulate hypotheses about why the patients cycle at different rates only in terms of cell kinetic parameters (e.g. cost and turnover), which are difficult to measure in practice. In contrast, the ABM revealed that differences may additionally manifest themselves in the spatial architectures of the tumours, a hypothesis that is more easily testable, for example, by analysing histology. (b) Comparison of the treatment trajectory predicted by both models for two patients for whom the ABM fits better (P39 and P62), and three patients for whom the ODE achieved better fits (P60, P101, and P99). We can see that the ODE is better able to capture the peak PSA values in some fast cycling patients. (c) ODE and ABM fits for a representative subset of the patients that were excluded from the cycling analysis, because the ABM was unable to capture the cycling dynamics. In all cases the ODE fails to recapitulate the dynamics of these patients, too and, in fact, these patients were also excluded from the analysis in our prior paper [5]. This suggests that the observed dynamics in these patients is not explicable by our simple 2population model, and that other factors, such as a third tumour cell population or non-tumour cells (immune cells or stroma) may be at play. (d) Comparison of the parameter values inferred by the two models (each point is one patient; n = 65). Inset text shows the Pearson's correlation coefficient. There is generally little correlation, likely because of differences in how these parameters act on the growth dynamics in the two models. The exception to this is turnover, indicating an important role in both models. Overall, the ABM results presented here corroborate and complement our prior ODE analysis by providing an interpretation of how the inferred differences in parameter values may manifest in the spatial architecture of the tumours.