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

Mathematical models used in preclinical drug discovery tend to be empirical growth laws. Such models are well suited to fitting the data available, mostly longitudinal studies of tumor volume; however, they typically have little connection with the underlying physiologic processes. This lack of a mechanistic underpinning restricts their flexibility and potentially inhibits their translation across studies including from animal to human. Here we present a mathematical model describing tumor growth for the evaluation of single-agent cytotoxic compounds that is based on mechanistic principles. The model can predict spatial distributions of cell subpopulations and account for spatial drug distribution effects within tumors. Importantly, we demonstrate that the model can be reduced to a growth law similar in form to the ones currently implemented in pharmaceutical drug development for preclinical trials so that it can integrated into the current workflow. We validate this approach for both cell-derived xenograft and patient-derived xenograft (PDX) data. This shows that our theoretical model fits as well as the best performing and most widely used models. However, in addition, the model is also able to accurately predict the observed growing fraction of tumours. Our work opens up current preclinical modeling studies to also incorporating spatially resolved and multimodal data without significant added complexity and creates the opportunity to improve translation and tumor response predictions. Significance: This theoretical model has the same mathematical structure as that currently used for drug development. However, its mechanistic basis enables prediction of growing fraction and spatial variations in drug distribution.

The paradigm di↵usion-limited growth model of tumour growth [S1, S2, S3] is based on the solution of a di↵usion model describing nutrient distribution within the tumour coupled to the equations of conservation of mass under the assumption of radial symmetry. We summarize the key points here.
The conservation of mass equation (with constant cell density) on the tumour domain is r · u(t, x) = ⇤(x, t, w), (S.1) where u is the velocity field describing the motion of cells within the tumour and ⇤ is a growth function representing a mass source in proliferating regions and mass removal in necrotic regions. Cell proliferation increases the cellular mass which translates into tumour expansion and cell movement through the vector u. For the drug treatment models below ⇤ depends on the drug concentration w as well.
In conjunction with the conservation of mass equation we solve the di↵usion equation for the concentration of nutrient in the tumour which is denoted by c.
As nutrient di↵usion is a fast process compared with tissue growth, the quasisteady approximation is valid leading to the di↵usion equation in the form D c r 2 c = k(c), (S.2) S-1 to be solved on the tumour domain, where k(c) is the cellular uptake rate.
Assuming a critical concentration c = c crit or equivalently a critical volume V ⇤ above which necrosis is initiated due to lack of nutrient, the cellular uptake is constant when V < V ⇤ but above this size a necrotic region forms so that cell uptake only occurs for growing cells.
The parameter V ⇤ can be expressed in terms of the other model parameters.
To see this, consider the case where the tumour is small enough that all the cell proliferate and consume nutrient. Assuming radial symmetry, the di↵usion equation (S.2) can be integrated twice to give the concentration of nutrient in the tumour as where k cons is the assumed constant rate of consumption, R T is the tumour radius, and c 0 is the concentration of nutrient outside the tumour. If necrosis occurs when c = c crit , we see this will first occur at the centre of the tumour (r = 0), when the tumour radius is R ⇤ = p 6(c 0 c crit )D c /k cons . Hence V ⇤ = 4⇡(R ⇤ ) 3 /3. The parameters D c and k cons enter the solution through V ⇤ and only V ⇤ is needed as a fit parameter of the model. Consider the full model, after assuming radial symmetry, we introduce two dynamic radii: the overall tumour radius R T (t) and the radius at which the nutrient level drops low enough for necrosis to occur R N (t) (only relevant if R T corresponds to a volume larger than V ⇤ ), see Fig. 1. The radius R N can also be expressed in terms of the growth fraction: Under the assumption of spherical symmetry, the conservation of mass and S-2 di↵usion equations reduce to 1 r 2 @(r 2 u r ) @r = 8 > < > : We can directly integrate the di↵usion equation (S.5) on the domain r  R T , with continuity conditions on c and @c/@r across the inner boundary, boundedness at r = 0. Two final conditions come from that at r = R T the concentration is c 0 which is the known concentration of nutrient outside the tumour and that at R N , c = c crit which is the minimal concentration before necrosis occurs (which as shown above may alternatively be expressed in terms of V ⇤ the minimal volume). The solution for the concentration is then Importantly, evaluating the solution for concentration at which can be be rearranged using Equation (S.3) and the expression for V ⇤ to give the growth fraction, which varies in time, as a function of the tumour volume through the equation Finally, integrating the growth equation (S.4), in a similar manner gives a solution for the velocity u. Evaluating this solution on r = r T where the velocity is dR T dt (t) gives the growth law for the overall volume V , which after eliminating R N using equation S.3 can be expressed in terms of the growth fraction as

Spatial drug distribution model
In a similar manner to the nutrient equation the spatial distribution of drug in the tumour is found by solving the di↵usion equation for the drug concentration w on the tumour domain. We again assume a quasi-steady state due to drug di↵usion being fast compared to tissue growth, a widely accepted assumption for cytotoxic chemotherapeutic agents [S4].
Assuming that drug uptake only occurs in the proliferative cells but with drug degradation continuing throughout the tumour, we get for the drug di↵u- Here D w is the di↵usion constant of the particular compound in tumour tissue, l is the decay rate of the drug within the tumour. The function g(w, x) is the form for the drug action which we assume here to be linear g(w, x) = K kill h w(t, x), within the proliferation part of the tumour, where K kill is the drug kill, and h > 0 is the e ciency of the drug, this is a measure of how much drug is used in the killing process: h ! 0 indicates that large amounts of drug are required in the killing process whereas large h is the opposite. As for the nutrient equation, we obtain (S.10) We require that the drug concentration as well as its derivative are continuous throughout the tumour domain as well as that the concentration is finite at the tumour centre, and denote the external drug concentration at the tumour by , which is assumed known. This boundary value problem can be explicitly solved in terms of spherical harmonic functions [S5] which gives the distribution of drug concentration within the tumour when V > V ⇤ as hDw . The parametersĀ, A and B are defined by Note that as the drug kill e↵ect needs only be considered in the proliferative compartment we only require the solution for w(r, t) only in R N < r < R T . When V < V ⇤ and GF = 1 the Substituting for w from equation (S.11), proliferative compartment only, into the conservation of mass equation (S.4) and integrating as before gives the modified system Note that if F = 0 (full drug distribution, hence w(t, r) = w ext (t)), the integral in (S.14) can be explicitly integrated and gives the relation for K V as in (4).

Numerical implementation
The di↵usion limited model Equations (2) to (3) (main paper) may be implemented as an algebraic di↵erential equation (ADE), however we choose to work with it as a system of ODEs. To this end we take the time derivative of the nutrient constraint equation (3), this gives We solve (S.16)-(S.18) using MATLAB(2019a) using the ODE45 solver for the exponential phase and the ODE15S solver for the non-exponential phase using S-5 an event solver to manage the switch between phases. A sti↵ solver is required due to the high sti↵ness that may occur during integration from high dose strengths.
The mixed e↵ects data fitting was implemented within MATLAB(2019a) using the NLMEFITSA optimisation routine. (With the necrotic radius R ⇤ used as a fit parameter rather than V ⇤ .) This performs Monte Carlo simulation designed to converge to the maximum likelihood estimates of the parameters.
Multiple initial parameter estimates were tested with no significant di↵erences in either the parameter estimation, curve fitting or run times being observed.
The initial parameter values used are listed in the GitHub code [S6]. A linear approach to the approximation of the log likelihood was chosen. This was due to it being much faster to converge to a solution compared to other approximate techniques (Gaussian quadrature and Importance sampling were also tested, but they did not give significantly di↵erent results and took far longer to run).
The default Matlab constant error model is used. We used a probit transform for the parameter transform which is a natural way to bound the parameters (a logit transform could also have been used). In particular, this is useful to ensure large values for the kill parameter are not tested, which occurs when using the common log transform. Choosing large, unphysical values of the kill term result in increased sti↵ness of the system and longer run times. The probit transform is defined as where erf is the error function. This probit transform is appropriate for all parameters except the rate of damage parameter k 1 for the exponential-linear model. This parameter is required to be greater than one and so it was scaled as indicated in Table 2. Where standard deviations are calculated these are             S-12