Reciprocal interactions between tumour cell populations enhance growth and reduce radiation sensitivity in prostate cancer

Intratumoural heterogeneity (ITH) contributes to local recurrence following radiotherapy in prostate cancer. Recent studies also show that ecological interactions between heterogeneous tumour cell populations can lead to resistance in chemotherapy. Here, we evaluated whether interactions between heterogenous populations could impact growth and response to radiotherapy in prostate cancer. Using mixed 3D cultures of parental and radioresistant populations from two prostate cancer cell lines and a predator-prey mathematical model to investigate various types of ecological interactions, we show that reciprocal interactions between heterogeneous populations enhance overall growth and reduce radiation sensitivity. The type of interaction influences the time of regrowth after radiation, and, at the population level, alters the survival and cell cycle of each population without eliminating either one. These interactions can arise from oxygen constraints and from cellular cross-talk that alter the tumour microenvironment. These findings suggest that ecological-type interactions are important in radiation response and could be targeted to reduce local recurrence.

The growth of tumour spheroids typically begins with a phase of exponential growth which is followed by an intermediate linear growth phase and eventually superseded by a slower approach to an equilibrium size at which the net rates of cell death and proliferation across the tumour volume are balanced. Growth laws that have been proposed to describe this process include the Gompertz, Bertalan↵y and logistic models [5].
Here we use the latter to model the growth of homogeneous spheroids. The logistic growth law states that the rate of change of the tumour spheroid volume V (t) (mm 3 ) at time t (day) is given by where the parameters r > 0 (day 1 ) and K > 0 (mm 3 ) represent the growth rate and carrying capacity, or equilibrium size, of the tumour and V 0 > 0 (mm 3 ) denotes the spheroid volume at t = 0. The initial value problem (1) can be solved to give .

Lotka-Volterra model of heterogeneous spheroid growth
We assume that when the control and resistant cell populations are co-cultured to form heterogeneous spheroids, their growth dynamics and interactions can be viewed as a Lotka-Volterra system. If we denote by V C (t) and V R (t) the volumes of the tumour occupied by the control and resistant cells respectively at time t, then their evolution is described by the following system of ordinary di↵erential equations: In Eqs. (2), the positive parameters r C (day 1 ) and r R (day 1 ) are the intrinsic growth rates of the two populations and the positive parameters K C (mm 3 ) and K R (mm 3 ) their carrying capacities (i.e., the equilibrium volumes to which the control and resistant populations would evolve if they were cultured as homogeneous spheroids). The parameters C (day 1 mm 1 ) and R (day 1 mm 1 ) respectively describe the e↵ect that the control cells have on the resistant cells, and vice versa. In the absence of prior knowledge about the nature of the interactions between the two cell populations, we do not restrict C and R to be positive or negative. Indeed, as the signs of C and R vary we can distinguish six types of interactions: these are summarised in Table 1. We note that, if C = R = 0, then Eqs. (2) reduce to two logistic equations for V C and V R similar to Eq. (1).
In order to arrive at a well defined initial value problem, Eqs. (2) are supplemented by the following initial conditions: For a homogeneous spheroid, whose growth can be modelled using Eq. (1), there is an obvious physical interpretation of the carrying capacity parameter K -it is the maximum volume of the spheroid that can be supported by its environment. For a heterogeneous spheroid, such as that modelled with Eqs. (2), the spheroid's saturation size cannot be described by a single parameter. Indeed the di↵erent interactions that may exist between co-cultured cell populations manifest themselves in a multitude of possible growth regimes. Nonetheless, we can consider some simple cases. Let us first denote by K T the carrying capacity of a heterogeneous spheroid, when it exists. If C = R = 0 then it is straightforward to show that In the case of competition, when C > 0 and R > 0, we would expect Similarly, when the populations support each other, that is when C < 0 and R < 0 (mutualism), we For logistic growth, the carrying capacity K corresponds to a stable steady state (see Eq. (1)).
Similarly, the co-existence equilibrium solutions of Eqs.
Steady states 1, 2 and 3 exist for all choices of the interaction parameters C and R . By contract, the coexistence steady state is only physically realistic and stable for certain combinations of values of r C , r R , then the coexistence steady state exists in the shaded regions of Fig. 1 where We note also that if ⌘ R = ⌘ C = ⌘, say, then

Logistic model
The values of parameters ✓ hom = (r, K, V 0 ) in Eq. (1) were estimated by minimising the weighted sum of squared residuals between experimental measurements of spheroid volume and the solution to Eq. (1) for given values of ✓ hom . In mathematical terms, we sought✓ hom such that where V data i are the volume measurements taken at times t i (i = 6 for PC3 and i = 9 for DU145 spheroids, We remark that there is no consensus regarding how to weight the data when minimising the sum of squared residuals. In equation (4), we follow [1] and use the model to weight the data.

Lotka-Volterra model
The values of the parameters r C , r R , K C , K R , V C0 and V R0 describing the growth rates, carrying capacities and initial volumes for each cell population were estimated using the volume data collected from the homogeneous spheroids, and were subsequently fixed in Eqs. (2). Thus the only unknown parameters in Eqs. (2) were the interaction parameters ✓ het = ( C , R ). These were estimated by seeking ✓ het that minimises the following sum of squared residuals: where V data C (t j ) and V data R (t j ) respectively represent the volumes of the control and resistant populations within the heterogeneous spheroids at times t j , V model are the total volumes of the heterogeneous spheroids at times t i and V model The times t i at which total spheroid volumes were measured are shown in Fig. 1b (main text). The times t j at which the proportions of the control and resistant cells were measured are t j = (5, 10, 15, 19) (days) for PC3 and t j = (5, 10, 15) (days) for DU145 spheroids.
The minimisation problems (4) and (5) were implemented and solved in MATLAB using a non-linear least squares solver lsqnonlin which is an implementation of the trust-region-reflective iterative optimisation algorithm [3]. Since lsqnonlin is a local solver we use it in combination with MultiStart, a function that runs the local solver from a number of randomly selected starting points within a prescribed region of parameter space, thus ensuring more thorough exploration of the parameter space.

Radiation response modelling for theoretical study
To simulate growth of heterogeneous spheroids after exposure to radiation we used Eqs. (2) for a range of values of C and R (with other parameter values fixed) to grow in silico spheroids until they reached a volume of V IR = 0.9mm 3 . We then simulated radiation damage by calculating the surviving fraction of each population according to the linear-quadratic model [9]. In particular we set for irradiation, and ↵ (Gy 1 ) and (Gy 2 ) are the lethal lesions made by one-and two-track actions of radiation dose d. Following radiation we used Eqs. (2) to regrow the in silico spheroids until they reached their pre-irradiation volume of 0.9mm 3 . Values of the parameters ↵ and were estimated for each cell population using data collected from the clonogenic assays. The resulting growth curves were used to generate surface plots representing the regrowth time ( Figure 4 in the main text).

Supplementary Methods 2: Spatially resolved computational model
We

Oxygen dynamics
A reaction-di↵usion equation models the concentration of oxygen within our simulation domain. If we denote the oxygen concentration at location x and time t by c(x, t) (mol cm 3 ) then its evolution is described by where D is the (constant) oxygen di↵usion coe cient (cm 2 s 1 ) and (x, t) is the oxygen consumption rate (mol cm 3 s 1 ). In what follows, we fix where 0 <  Q <  P denote respectively the constant rates at which quiescent and proliferating cells consume oxygen (units: moles cm 3 s 1 ). Equation (8) is supplemented by the following initial and boundary where L is the domain length and c 1 is the background O 2 concentration. Equations (8) The state of a cell at location x and time t, denoted as state(x, t), is determined by the concentration of oxygen, c(x, t), at that location. If we define by c Q and c N the threshold oxygen levels below which cells of type become quiescent and necrotic, respectively, then we have: In addition to a state, every automaton element has a neighbourhood associated with it. The neighbourhood provides local information about the cell's surroundings and can impact the cell's next state. The two most commonly used neighbourhoods in CA models are the von Neumann and Moore neighbourhoods (Fig. 2). The number and location of neighbours can have a significant impact on cell-cell communication.
Since it has been previously demonstrated that the Moore neighbourhood can minimise artefacts associated with lattice anisotropies [10] we use the first order Moore neighbourhood. Thus, for a 2D model, every cell communicates with its eight nearest neighbours.

Cell cycle progression
Every automaton element occupied by a proliferating cell is assigned a counter that monitors its cell cycle progression. When a new cell is created following cell proliferation, its counter ⌧ cycle is a random number drawn from a normal distribution with mean⌧ cycle and standard deviation cycle where⌧ cycle represents the average cell cycle duration for cell population . After each discrete time step ⌧ , the cell cycle counter ⌧ cycle for a given cell is reduced by an amount that depends on its local neighbourhood (see Fig. 3). This mechanism imitates contact inhibition of proliferation, a well known feature of 2D and 3D cell aggregates.
Unlike Jagiella et al. [7] who assumed that a cell cycles only if it is located within a given distance from the spheroid boundary, here a cell makes this decision based on local information only.

Cell division
When ⌧ cycle  0 for a given cell, it divides to produce two identical cells. One cell occupies the same position as its parent and the other cell is placed in an adjacent automaton element. If more than one automaton element adjacent to the dividing cell is empty then the division process is complete (if a dividing cell has more than one free neighbour then the neighbour with the maximum number of neighbours is chosen to maintain cell-cell adhesion [8]). Alternatively, if a dividing cell does not have an empty adjacent automaton element, then we find the shortest chain of cells that connects the dividing cell to the spheroid's boundary.
We shift the cells in the chain along this path in order to create space for the new daughter cell (if multiple shortest paths exist we choose one at random). This algorithm mimics the way in which growing cells exert mechanical stress on their neighbours to generate spheroid expansion. Fig. 4 illustrates how the chain of cells 1, 2 and 3 is shifted towards the spheroid's boundary to create space for the new daughter cell D.

Cell death
A cell becomes necrotic if the oxygen concentration at its location falls below a threshold value c N . Necrotic cells are lysed at the rate p lys (h 1 ). When a necrotic cell is lysed it is removed from the computational grid and a chain of cells is shifted from the spheroid's boundary towards the location of the removed cell (see Fig. 5). In order to preserve spheroidicity we select the cell on the boundary that is located furthest away from the spheroid center (or choose randomly if multiple cells lie at the same distance from the center).
The resulting cell rearrangement occurs within one computational time step and ensures that the spheroid remains compact.

Algorithm summary
We now outline the computer algorithm used to implement the hybrid CA model.

1.
Initialisation (t = 0). We assign values to all model parameters (see Table 2), initialise the oxygen distribution on the spatial grid, and place a cluster of proliferating control and resistant cells at the centre of the computational grid (all other grid sites are assumed to contain culture medium), assign cell cycle duration times to the cells and set t = ⌧ where ⌧ is the length of our computational step.
2. Nutrient consumption. We update the oxygen concentration by solving Eqs. Hybrid CA model parameter estimates

Oxygen concentration
The background oxygen concentration was set to c 1 = 2.8 ⇥ 10 7 mol cm 3 [4] and the oxygen di↵usion constant within a growing spheroid was set to D = 1.8⇥10 5 cm 2 s 1 [6]. Since the di↵usion of O 2 in culture medium is likely to be higher than within the spheroid [7], we assumed that the concentration of O 2 within the medium was replenished on a faster timescale than within the spheroid. In practice this means that O 2 levels in the culture medium were held constant at the background O 2 concentration c 1 . Proliferating cells were assumed to consume oxygen with rate  P . Quiescent cells also consume oxygen but at a lower rate

Measurements of oxygen consumption by the PC3 Ctrl and PC3 Res cells cultured in monolayers
revealed that the Ctrl cells consume more oxygen than Res cells by a factor of 1.4 (see main text). As shown previously, cells in tumour spheroids are known to consume less oxygen than in monolayers [4]. In the absence of estimates of the O 2 consumption rates in the PC3 spheroids, we treated  P as free parameters and explored their influence on spheroid growth. Estimates of the values of c Q and c N were also di cult to obtain and so they were also treated as free parameters (see Sec. ).

Cell cycle duration
We

Parameter inference
We used a two step process to estimate the values of parameters c Q , c N ,  P and p lys for the PC3 Ctrl and PC3 Res cells. First, we used our hybrid CA model to simulate the growth of homogeneous spheroids for a range of values of the model parameters. We then compared the simulation results to growth curves from the PC3 Ctrl and PC3 Res homogeneous spheroids and experimental images of their spatial distributions.
In this way, we identified those values of c Q , c N ,  P and p lys that best fit the data. We then fixed these parameter values in the hybrid CA model when simulating the growth of heterogeneous spheroids.
In more detail, to identify the values of c Q , c N ,  P and p lys that best fit the growth curves for a given cell population 2 {PC3 Ctrl, PC3 Res} we ran the CA model for a single population for the following ranges of parameters: c Q 2 (0, 2.8 ⇥ 10 7 ), c N 2 (0, c Q ),  P 2 (1.0 ⇥ 10 9 , 1.0 ⇥ 10 8 ) and p lys 2 (0, 0.5). We compared the in silico growth curves to the experimental PC3 Ctrl and PC3 Res growth curves by calculating the sum of squared residuals between them and shortlisted ten sets of parameters that resulted in in silico growth curves with the lowest sum of squares for each population. We then compared the spatial distributions for the shortlisted in silico spheroids to the available images of homogeneous PC3 were subsequently used in the simulations of heterogeneous spheroids (see Table 2). Fig. 7 shows typical results from simulations of our hybrid CA model. Fig. 7a shows how the spatial distribution of a homogeneous spheroid changes over time for a single realisation of the model whereas

Cell cycle analysis
We set out to test the null hypothesis that the proportions of cells in G0G1, EdU+, and G2M phases do not change when resistant and control cells are cultured in homogeneous or mixed spheroids after five, ten, and 15 days of culture. We inspected cell distributions shown in Figure 1 and identified expdate as a potential confounding factor. For each combination of irradiation status (TRUE or FALSE), cell type (pc3 and du145), cell resistance status (resistant and control), and days after culturing (5, 10, and 15) we performed the following analysis. Due to overdispersion in the data, a Poisson regression model was not appropriate, and we instead used a negativebinomial model [6] to model cell counts (counts) under the alternative hypothesis using the model specification in Wilkinson-Rogers notation [7]: counts ⇠ group + phase + expdate + group:phase where group, phase, and expdate are main effects, and group:phase is an interaction effect between group and phase. group represents the type of culture (homogeneous or mixed), phase represents the cellcycle phase (G0G1, EdU+, or G2m), and expdate represents the date each cell count experiment was performed.
We generated half-normal plots of model fit residuals and simulated envelopes for these plots (data not shown) using the hnp package [3] in R [5]. The half-normal plots showed one model fit with unsatisfactory residual deviances, and this model was removed from further analysis. We tested the null hypothesis by testing if the inclusion of the interaction effect group:phase significantly reduces the deviance of the model using a Chisquare test as implemented in the anova.negbin function in the MASS package [6]. We corrected the 23 p-values obtained in this way for multiple testing with the Benjamini-Hochberg procedure [1]. We report six significant tests at a false discovery rate of 0.05 in Table 1.

Cell death rates
We next investigated whether the null hypothesis that the proportions of dead and alive cells did not vary between homogeneous and mixed spheroids. The cell counts shown in figure 2 indicate that these proportions vary at day 15 and to a lesser extent at day 10, and that experiment date may have an effect on death rates. Due to overdispersion in the data, we fit a quasibinomial model [6] with the specification counts ⇠ group + expdate using the glm function in R, where group and expdate are defined as in the previous analysis, and counts is a matrix of alive and dead cell counts. We generated half-normal plots with simulated envelopes as in the previous analysis (data not shown) and excluded two model fits due to poor fit on the basis of these plots. We tested the null hypothesis by testing if the inclusion of the group effect significantly reduces the deviance of the model using an F test as implemented in the anova.glm function in the stats package [5]. We correct the remaining 22 model fit p-values for multiple testing using the Benjamini-Hochberg procedure [1]. We report 9 significant tests at a false discovery rate of 0.05 in Table 2

Post-radiation growth experiments
Volumes for spheroids irradiated on day 4 after seeding were available for PC3 and DU145 spheroids seeded from several mixtures of radio-resistant and radio-sensitive cells irradiated between 2.5 Gy and 20 Gy (See Fig 4 in main text), as well as for unirradiated cells. All spheroid volumes were log 10 transformed. Logarithmic spheroid increase was calculated for each sampling date after irradiation by subtracting the log 10 volume on day 4 from the log 10 volume measured on that day. Baseline growth delay for each cell mixture and experiment was estimated as the time point at which unirradiated spheroids crossed a cutoff growth increase (2.5x, 3x, and 4x for data presented in Fig 4a, b, and e, respectively) by fitting a local polynomial regression [2] (exact method) and estimating the fit's intersection with the cutoff line through bisection [4]. Time until regrowth to the cutoff volume was estimated as the first sampling date on which a spheroid had a larger growth increase than the cutoff. Growth delay for irradiated spheroids was estimated as the time until regrowth minus the baseline growth delay. Irradiated spheroids with a negative growth delay were excluded from the analysis. Spheroids that did not regrow past the regrowth cutoff were censored.

Co-culture transwell experiments
We analyzed cell counts from a series of co-culture experiments of resistant and control cells in hypoxia and normoxia on six different days. Each experiment was performed in triplicate. Figure 3 shows that in hypoxia, counts in the non control/control groups (C/R, R/C, and R/R) may be elevated compared to the control/control co-culture (C/C). To establish whether cell counts were elevated in non-C/C co-culture in hypoxia, and whether cell counts did not differ between co-culture types in normoxia, we modelled cell counts (count) as depending on the day of the experiment (expdate) and the type of co-culture (group) using the model specification counts ⇠ group * condition + expdate Due to overdispersion in the data, we fit this model using a negative-binomial model. Half-normal plots showed a good fit of the model to the data (data not shown). We used the R function summary.negbin from the MASS package [6] to perform a Wald test on the group effects of the model. In normoxia, no culture type had a statistically significant difference in counts compared to C/C. In hypoxia, all three culture types differed from C/C (C/R: p=0.0038, z=3.54; R/R: p=8.7e-05, z=4.23; R/C: p=0.0036, z=2.86).  Figure 3: Cell counts of co-culture experiments. Each point represents a technical replicate. Each color represents an experiment date. In hypoxia, cell counts are somewhat higher in C/R, R/C, and R/R compared to C/C.