A Spatially Resolved Mechanistic Growth Law for Cancer Drug Development Predicting Tumor Growing Fractions

Preclinical evaluation of drug efficacy plays a fundamental role in the development of oncological treatments, with the aim being to predict pharmacologically active drug concentrations and guide dose exploration in the clinic. Data for these studies come from longitudinal measurements of tumor volume in animal models with specific targets investigated by the use of both cell-derived xenograft (CDX) and patient-derived xenografts (PDX); ref. 1. Central to these preclinical studies are mathematical models used to describe the tumor dynamics, fit the experimental data, and evaluate the antitumor effect. These models usually take the form of simple growth laws for tumor volume. There are many such models which can satisfactorily capture the dynamics (2–7) so that for any given dataset it is in practice very difficult to distinguish between them (4–6, 8–10). This is compounded by the fact that typically these growth functions are purely empirical descriptions of the data, not founded in a mechanistic description of the physiologic process.

The extension of growth law modeling in preclinical trials into mechanistically based modeling would have significant advantages. It could aid discrimination between models as models could be chosen by mechanism rather than only by fit. It would also enable key physiologic processes to be better integrated in trial design and encourage the use of a broader range of data beyond tumor volume. Indeed, current preclinical growth laws typically do not even account for the variations in growth fraction dynamics and hypoxia that have been shown to be key to determining treatment response (11–14). Mathematically, there have been significant advances in mechanistic models of tumor physiology, with simulations now possible even at the level of entire vasculatures (15). However, from a practical perspective, these more complex models raise significant additional challenges, which inhibits their incorporation into current preclinical workflows. They typically are based on a detailed knowledge of, for example, tumor geometry, which is not typically acquired in preclinical trials, and require significant computational simulations. In addition, they will have more parameters, which makes it challenging to parameterize them. This raises concerns of parameter identifiability, although this may not necessarily inhibit their predictive application (16).

Here, we present a mechanistic and spatial mathematical model for tumor growth based on a well-accepted mathematical framework coupling spatial diffusive processes to tumor growth (17). Significantly, its spatial underpinning means that it predicts the size of the growing fraction of the tumor over time. Under the assumption that the tumor is spherical we show that this model can be expressed as a growth law of a similar type to that currently used in the pharmaceutical industry. It thus can be fitted and analyzed using current industry-standard methods. We demonstrate the ease of use by comparison with a current commonly adopted growth law and validate both the models using standard methods on CDX and PDX data. Another advantage of the model is it can predict tumor growth fraction, which we here validate by endpoint histology. The new framework offers further significant advantages such as being able to account for spatial gradients in drug distribution. Together, these demonstrate that mechanistic models can be made fit-for-purpose to be incorporated into the current preclinical trial workflow and can enable the incorporation of data with spatial resolution as this is more routinely obtained.

Mathematical models as used in preclinical studies typically can be expressed as a rate equation for the tumor volume of the form

with the initial tumor volume , a fit parameter. The function describes cell loss as a function of the drug concentration , most often which represents exponential loss of tumor volume. The function represents the net growth rate, with a range of different forms used in preclinical studies (refs. 3, 6, 18; Table 1) with the form often optimized for data fitting. For example, where growth is observed to be linear , whereas exponential growth is captured by constant. Although all the models from Table 1 have successfully been used to fit preclinical data, it is often challenging to distinguish them in terms of quality of fitting (5, 6, 8–10). However, there are some key structural differences. The first five (linear, exponential, logistic, Gompertz, and exponential linear) are commonly used empirical growth laws. The surface growth model (19) and the proliferative rim model (20), are less widely adopted, and while mechanistic in foundation, each assume proliferation is restricted to a fixed unchanging region of the tumor. As such these models can only be expected to mechanistically approximate the tumor dynamics in specific circumstances and for limited times. A significant drawback of all the growth models in Table 1, except the Gompertz and logistic models, is that they predict unbounded growth as time increases, which is not experimentally observed. However, the listed models may still adequately describe animal-derived data, as the data cannot be obtained over long enough timescales for any growth retardation to be observed. However, in clinical studies where the observation timescales are longer, it is less clear that this would hold.

An alternative approach to modeling tumor growth outside of pharmaceutical drug development has been to model the spatial diffusion of nutrients within tumors directly. It is assumed that nutrient availability determines the status of the cell, with cells experiencing low nutrient levels undergoing necrosis, a form of cell death. This results in a map of the tumor with different sized cellular compartments for, for example, proliferative and nonproliferative necrotic (dying) compartments determined responsively from the environmental conditions (17, 21), see Fig. 1, as is observed experimentally (22). The size of the compartments thus alters dynamically and is determined organically from the tumor growth, generating predictions for the size of each compartment over time. This basic framework has been extensively built upon in the mathematical literature (23, 24). The increasingly complex models are typically expressed in terms of partial differential equations which are either solved numerically themselves or coupled into finer scaled simulations.

We adopt this diffusion-based framework and make simplifying assumptions appropriate to the preclinical context. We thus reduce the model to a growth law of the form Eq. (1). The method is based on assuming a spherical tumor with the nutrient penetration within the tumor determined from the diffusion equation. From the diffusion equation, we then obtain the growth fraction of the tumor , that is, the volume fraction that is proliferating, which varies in time and is determined from the tumor volume. When this represents a fully proliferating tumor. A key parameter is the volume at which the tumor first experiences necrosis. The resulting model is similar to Eq. (1), but with only the growing fraction contributing to the tumor growth and with the drug affecting only living cells, thus

with when . The growth equation, Eq. (2), is coupled to the equation for growing fraction , Eq. (3), and these two equations are solved together. In this model and the initial condition is . Specifically, is the cell replication rate, is the rate at which cells die in the proliferating region, and is the rate of breakdown and removal within the necrotic region. While the volume of the tumor , the entire tumor proliferates, and the model predicts exponential growth. As tumor size increases eventually and the model switches to diffusion-limited growth, in which the tumor has a necrotic core surrounded by a proliferating growing rim, with the growth fraction of the tumor determined from Eq. (3; see Supplementary Information for further details).

The formulation Eqs. (2) and (3), which we term the diffusion-limited model, is of the same mathematical form as currently implemented. However, it is now possible to directly predict the growth fraction of the tumor. We will demonstrate the functionality of the description Eqs. (2) and (3) both for control datasets and for drug trial data. We focus here on preclinical trials of drug treatments. We take , which is similar to the standard and broadly applicable exponential loss term (25); however, loss now only occurs in the proliferating compartment. We see that the mechanistic approach allows us to incorporate that cell cycle–specific pharmaceutical agents target proliferating cells. The pharmacokinetic description is in general compound specific. Here we consider CPT-11 (Irinotecan) which has been shown to be adequately described by an exponential model (14). We thus set for the fitting , with the known administered dose and the decay rate, leaving as the single parameter to be fitted describing drug effects. With the half-life of Irinotecan being 12 hours, .

Equations (2) and (3) form a growth law and two constraints which can be numerically solved as an algebraic differential equation. However, given current fitting protocols it is easier to work with a system of differential equations. This may be obtained from Equations (2) and (3) by differentiating the constraint Equation (3). All routines are implemented within MATLAB2019a (MathWorks) and are made available in GitHub at https://github.com/DrAdamNasim/Diffusion_Limited_Cancer_Growth_Model. Because of the significant intersubject heterogeneity, we use a nonlinear mixed effects fitting approach (NLME) to fit both the population and individual subject parameters using the stochastic approximation expectation-maximization algorithm (see Supplementary Data).

Control data are from AstraZeneca mouse-derived CDXs for two cell-lines SW620 (29 mice), epithelial colorectal adenocarcinoma, and Calu6 (178 mice), epithelial lung adenocarcinoma. Longitudinal time series of subcutaneous tumor volume were calculated by manual measurement of length and width using calipers (). With regards treatment data, we consider SW620 CDXs (95 mice) which were treated with weekly doses of 50 mg/kg of CPT-11 (Irinotecan) either three times on days 1, 8, 15 (protocol 1, 68 mice), or four times during the experiment on days 1, 8, 15, 22 (protocol 2, 18 mice) or 4, 11, 18, 25 (protocol 3, 9 mice).

In addition, we use the Novartis PDX dataset, which is the largest publicly available database of PDX control data (26). The Novartis data contains 226 mice with six different tumor types: breast carcinoma 39 mice, non–small cell lung carcinoma 28 mice, gastric cancer 44 mice, colorectal cancer 45 mice, cutaneous melanoma 33 mice, and pancreatic ductal carcinoma with 37 mice.

All animal studies in the United Kingdom were conducted in accordance with the UK Home Office legislation, the Animal Scientific Procedures Act 1986, and with AstraZeneca Global Bioethics Policy.

All codes are made freely available in GitHub at https://github.com/DrAdamNasim/Diffusion_Limited_Cancer_Growth_Model.

To demonstrate the fitting of the diffusion-limited model to longitudinal data, we consider the AstraZeneca CDX data in Fig. 2. In Fig. 2A–C, we show the individual fits of the diffusion-limited model to the experimental data separated by protocol. The diffusion-limited model is observed to fit well across the range of curves. Similarly, qualitatively good fits are observed for the other CDX and PDX control datasets (Supplementary Fig. S1). To further demonstrate the quality of fit, we perform a visual predictive check (VPC; ref. 27). This confirms that the diffusion-limited model captures the full range of dynamics (Fig. 2E and F). Similar results are also obtained from VPCs for the other control PDX and CDX datasets (Supplementary Fig. S2).

Looking quantitatively at the NLME results (Table 2), we see a small residual mean square error. We further compare the diffusion-limited fit to results obtained from fitting one of the most widely used models, the exponential-linear model (2, 14, 28), which has the same number of parameters as the diffusion-limited model. The parameter estimation for both the exponential-linear and diffusion-limited models is given in Table 2. Figure 2D shows five representative fits. We clearly see the similarity in the model behavior. This similarity in the quality of fit is confirmed by the closeness of the Akaike information criterion number (−217 and −221, for the diffusion-limited and exponential-linear, respectively) and root mean square error (0.146 and 0.142, respectively; Table 2). A similar result is achieved when considering the control datasets. These results are summarized in Supplementary Tables S1 to S8. Given the typically limited time course of preclinical data, it has been shown across numerous studies that no model may be optimal across all datasets (3, 6, 9). However, the results indicate that the diffusion-limited model provides at least a similar standard of fit for typical data.

As discussed, the rate equation predicts the growth fraction through Eq. (3). Thus, in addition to being able to fit the model to longitudinal measurements of tumor size, we can now predict how the growth fraction changes, including in response to treatment. We show, for example, in Fig. 3A that those tumors showing greatest response (dosing strengths 75 and 100 mg/kg) are predicted to have the largest growth fraction, with consequent implications for treatment. The growth fraction dynamically alters as the tumor volume changes, that is, as the tumor volume decreases so the diffusion equation predicts that more of the tumor will have access to nutrient responsively increasing the growing fraction accordingly.

The growth fraction of tumors is normally unavailable due to the difficulties of accessing this information in vivo. However, at the termination of xenograft experiments, it is possible to obtain growth fraction data through histology, although this is not routinely done. We consider independent datasets for the necrotic area for untreated SW620 and Calu6 cell lines obtained from histology of bisected tumors (Ki-67 staining). The total volumetric growth fraction , which varies in time, and this and the necrotic area of the cross-section are related by . We simulate tumor dynamics for a tumor using the population parameters obtained from the mixed effects fitting of CDX data (Supplementary Tables S1 and S2). These simulations predict a necrotic area at the end of the experiment for both the SW620 cell line and the Calu6 cell lines that fits the experimental data (lying within the range of the data in Fig. 3B and C). In this case, the necrotic area data are an independent dataset. We could alternatively have used this data to further calibrate the model parameters. This highlights one of the benefits of a mechanistic model, that we can utilize other sources of information to better calibrate the model by way of multimodal fitting.

In fitting the CDX treatment data (Fig. 2), we assumed that the drug was available throughout the tumor at the same concentration, although of course the killing effect was only experienced by the proliferating cells. The drug concentration thus decreased in time but without spatial variation. For small molecule drugs, such as Irinotecan, this is likely a good approximation with good fits being obtained to the data. Indeed, nonspatial tumor kill models have been shown to provide fit to clinical treatment data well indicating broad applicability (25). However, for significantly larger molecule drugs, or for engineered drug delivery systems, the ability to deliver the drug to the tumor core by diffusion is likely to be restricted. Indeed, research into antibody–drug conjugates (ADC) already report this shielding effect as a potentially significant issue in treatment effectiveness (29).

The spatial distribution of drugs in the tumor can be incorporated in the model framework using the same diffusive framework as for nutrients. This extended model, while more complex, is equally straightforward to numerically implement requiring a single additional integration. As a specific example, we consider the case where cell loss only occurs in proliferating cells. In this case, the cell loss is integrated over the proliferating compartment so that

Where is the tumor radius, with for a spherical tumor . When the drug concentration is the same throughout the tumor we recover from Equation (4). For a drug diffusing into the tumor, however, the distribution is not constant and in this case the concentration in the proliferating region (where the drug acts) is (see Supplementary Data)

For , where is the time varying applied drug concentration. The additional parameters and are known and defined in terms of and , see Supplementary Information. This reduces to when and . The parameter effectively quantifies the depth of penetration of the drug. A larger corresponds to more limited drug penetration, while as we regain a uniform drug distribution throughout the tumor. large corresponds to a surface kill. In terms of the other parameters , where is the drug decay rate, the drug efficiency, and the drug diffusivity.

The solution is then substituted into Eq. (4) to determine . Although more complex than the previous model, it is easily implemented numerically. Heatmaps showing drug concentration within the tumor clearly demonstrate the effect of (Fig. 4A) As increases the amount of drug in the tumor drops significantly, indeed by the time the concentration at the inner necrotic radius drops to 34% percent of its concentration at the outside edge (Fig. 4B). Exploring the effect of spatial drug distribution, we plot a representative solution of the spatial drug model (Fig. 4C). We see in the tumor growth curves (blue lines, Fig. 4C) that as increases and the drug penetration is reduced (going from bottom to top) that the effect of the drug on the tumor growth is reduced, with less response at each dose. With greater drug penetration ( small), a greater effect is observed with faster dynamic rebounds. Overall, it is clear tumor response will vary depending on the drug distribution.

We present a mathematical model for tumor growth and treatment that is based on a mechanistic description of nutrient limited growth. Despite being derived from spatial partial differential equations, we show that these equations can be reduced to a growth law which is not significantly more complex than those currently used for preclinical drug trials. The work extends the growth modeling options available for preclinical drug trials enabling spatial effects to be explicitly taken into account.

The development of a mechanistic approach has several advantages over more phenomenological growth laws. Perhaps most significantly, in its current form the model allows for the dynamic prediction of the growing fraction of the tumor, accounting for cell loss from, for example, hypoxia. As most cytotoxic agents target only actively proliferating cells, tracking the growth fraction has potentially significant implications for treatment dynamics. This is only enhanced by the increasingly clear role hypoxia has in inducing downstream biological processes which directly promote tumor resistance to treatment (30–34). We have focused on preclinical trials of drug treatments and similar therapeutics. To model other types of cancer therapies, such as immunotherapies, the model would have to be adapted. However, the underpinning growth model could be incorporated into other modeling frameworks with different models of treatment.

By considering a range of different datasets including PDX and CDX data, we have shown that the diffusion-limited model fits the longitudinal volume data comparably with current industry-standard models. For independent datasets containing growth fraction data at termination, we also predict a growth fraction commensurate with that observed. This without presenting significant increased complexity in the numerical solution. In addition, the diffusion-limited model parameters are all based on observable phenomena: net proliferation rate, maximum volume before necrosis begins, initial volume and cell loss through necrosis. The introduction of spatial modeling into a growth law framework has also enabled us to present an additional extension to the model which enables the spatial distribution of drug throughout the tumor to be explicitly modeled. The indicative results obtained demonstrate the importance of accounting for this in preclinical studies where there is evidence of reduced penetration into the tumor. Indeed, the increasing importance of mAb treatments, with their reduced diffusivity, and ADCs (35) will only increase the importance of spatial modeling in preclinical trials.

A. Nasim reports grants from AstraZeneca during the conduct of the study. No disclosures were reported by the other authors.

A. Nasim: Formal analysis, investigation, methodology, writing-original draft, writing-review and editing. J. Yates: Conceptualization, resources, formal analysis, methodology, writing-original draft, writing-review and editing. G. Derks: Conceptualization, formal analysis, methodology, writing-original draft, writing-review and editing. C. Dunlop: Conceptualization, formal analysis, methodology, writing-original draft, writing-review and editing.

A. Nasim thanks the University of Surrey's Doctoral College Grant and additionally AstraZeneca for funding. All authors thank AstraZeneca for generous access to data and the EPSRC funded QSP-UK network (EP/N005481/1) for facilitating initial discussions.