Oxygen-Driven Tumour Growth Model: A Pathology-Relevant Mathematical Approach

Xenografts -as simplified animal models of cancer- differ substantially in vasculature and stromal architecture when compared to clinical tumours. This makes mathematical model-based predictions of clinical outcome challenging. Our objective is to further understand differences in tumour progression and physiology between animal models and the clinic. To achieve that, we propose a mathematical model based upon tumour pathophysiology, where oxygen -as a surrogate for endocrine delivery- is our main focus. The Oxygen-Driven Model (ODM), using oxygen diffusion equations, describes tumour growth, hypoxia and necrosis. The ODM describes two key physiological parameters. Apparent oxygen uptake rate (kR′) represents the amount of oxygen cells seem to need to proliferate. The more oxygen they appear to need, the more the oxygen transport. kR′ gathers variability from the vasculature, stroma and tumour morphology. Proliferating rate (k p) deals with cell line specific factors to promote growth. The K H,K N describe the switch of hypoxia and necrosis. Retrospectively, using archived data, we looked at longitudinal tumour volume datasets for 38 xenografted cell lines and 5 patient-derived xenograft-like models. Exploration of the parameter space allows us to distinguish 2 groups of parameters. Group 1 of cell lines shows a spread in values of kR′ and lower k p, indicating that tumours are poorly perfused and slow growing. Group 2 share the value of the oxygen uptake rate (kR′) and vary greatly in k p, which we interpret as having similar oxygen transport, but more tumour intrinsic variability in growth. However, the ODM has some limitations when tested in explant-like animal models, whose complex tumour-stromal morphology may not be captured in the current version of the model. Incorporation of stroma in the ODM will help explain these discrepancies. We have provided an example. The ODM is a very simple -and versatile- model suitable for the design of preclinical experiments, which can be modified and enhanced whilst maintaining confidence in its predictions.


S1. ODM history description S.1.1 Initial ODM
Our motivation comes from previous mathematical models based on hypoxia and concretely a model developed by Evans et al. [7], in which the compartments adopt a pseudophysical position in space. In general, our approach is a compilation of preexisting modified models. Under extreme exposure to starvation and low oxygen levels, cells undergo nec rosis which is typically localised in the core of the tumour. Necrotic death is hereby described in our model by the death rate ( ). These necrotic cells may be engulfed by the surrounding tissue, although to our knowledge, there has been no evidence reported in the literature. Further, although necrotic tissue may trigger an inflammatory response with subsequent macrophage recruitment, such a response has rarely been observed in immunocompromised animal models (cell components should migrate across the endothelium, see [s1, 2]). In this model, necrotic cells are not removed by the immune system and contribute to the tumour volume.
On the other hand, highly proliferating cells undergo numerous intrinsic processes involving mutations and DNA damage, which may activate the mitochondrial apoptotic pathway [s3]. However, the actual fate of the apoptotic cells is still unclear. They could remain, be reabsorbed or transported to other parts of the organism [s4]. In our model, the parameter ( ) represents the apoptotic rate of cells, which were not depleted by phagocytosis from the host immune system.
In our model, the communication between cell cycle and quiescent-hypoxic state is represented by a reversible reaction, where "deactivation", ( ) and further "reactivation", ( 1 ) are reciprocally possible. This allows future implementation of synergic therapeutic combinations, namely radio-chemotherapy [s5].
Bringing all of these points together means the time rate of change of the proliferating volume is given by state equation (S1) whilst that of the quiescent region is state equation (S2) and necrosis is state equation (S3). Fig A shows a graphical representation of the model. Note that the "reactivated" quiescent cells are evenly distributed across all the proliferating layers, thus NP appears on equation S1.
The parameters of the model for tumour growth description are therefore = Disclaimer: Two different indexes have been defined here: -i: refers to the natural counter of the shells -j: refers to the index of the re-discretised layers.
In this model, the layers lose their integrity in every iteration, so P, works as a single compartment, which will be subsequently reorganised in different layers. So the interchange succeeds between 2 compartments P and Q effectively.
Knowing the algebraic equation describing oxygen partial pressure (eq. 6 of the main paper), the model described in eq. (S1) can transform into, where in this case the proliferating constant has been regrouped as " = ′ • / with units h -1 . Inversely the oxygen levels in the surrounding tissue have become apparent A new term has also been added, substituting the whole tumour volume by the sum of its internal terms as follows: The apparent oxygen uptake rate now has adopted the form: The ODM model presented here is more flexible and provides with better results than the ODM presented in the main paper, but the relatively large number of parameters, makes the system not practically identifiable. However, in certain circumstances, this model would be more useful than the one presented in the main text of the paper.

S1.2 Advanced ODM
However, this rather complicated formulation can be simplified if we assume that hypoxia occurs at "a state of low cellular oxygen". Since we calculate oxygen tension spatially, we can set up a "threshold" at which hypoxia occurs (in practical terms this is the oxygen levels at which hypoxia markers start to be expressed). We defined this "switch" as a smooth sigmoidal function. If we extend this definition to necrosis we can simulate the 1D maps of oxygen and hypoxia/necrosis ( Following, the version of the ODM presented in section S1.2 is just an extended version of the discretisation throughout the whole tumour combined with the sigmoid switches for hypoxia and necrosis, resulting in the model presented in the main paper.

S2. Assumptions of the ODM
We revisit here the main assumptions of the model presented in the main paper. This time we want to explore the mathematical spaces of the assumptions and ground our decisions.

S2.1 Oxygen Diffusion
We will first describe part of the rationale of oxygen diffusion in the model. There are two possible kinds of hypoxia characterisation: acute (perfusion related) or chronic (diffusion related). It is important to note that, different markers of hypoxia might be looking at different kind of hypoxia [s6], which suggests the need to analyse the problems on a case specific basis.
Let us define concentration as oxygen partial pressure 2 = and consider the quasisteady state approximation 2 = 0. If diffusion is independent of position, D ≠ f(r) we derive the eq. (S7).
Above, D, PO2 and RO2 are respectively the diffusion coefficient, oxygen partial pressure and oxygen uptake rate by the cells. Rearranging the equation and assuming spherical symmetry, we obtain the equation 10 of the main paper.

S2.2 Limited Oxygen Distribution
Having introduced the oxygen diffusion problem we proceed with the ODM specif ic  vivo data. We see after the first 4 days, the growth in vivo and in vitro are equivalent as shown in Fig D. However, somewhere between day 4 and 7 there is an inflection point, which in the literature is often referred to as "the angiogenic switch". In our opinion, this is the switch from a cell suspension to a real organ, where oxygen delivery occurs from diffusion limited blood vessels rather than directly from the medium . This involves, not only development of a tumour vasculature, but also recruitment of stromal cells creating a favourable microenvironment.
In a xenograft, the tumour is only in contact with the well vascularised epithelial tissue of the lower ectodermic layers from one side and the body wall from the other side. As opposed to the vascularised somatic tissue, the angiogenic vasculature of the tumour may not carry blood cells, and often transports low oxygen/nutrient plasma. Therefore, we can consider the effective access to blood carrying vessels is corrected through a shrinking lateral area-volume ratio (in other words, access to oxygen decreases with tumour growth). In the clinic, the situation may be reversed, thus the vasculature may be mature and tissue invasiveness may lead to the opposite trend (access to oxygen increases with tumour development).
Above, VTo and VT are the initial and final volumes of tumour respectively. The "shrinking supply" phenomenon has been vastly reported in the literature.

S2.3 Spatial Discretisation
In this model, as will be shown later in this section, the discretisation plays two main roles. One is the reduction of the total error associated with it and the other is transforming the model into a structurally identifiable model.
Hereby, there are 3 possibilities to discretise the model dividing evenly:  Volume: thickness would be uneven and proportional to its mean oxygen concentration.
 Radius: Error still increases with the gradient of PO2  Oxygen Partial Pressure Change: constant error function.
In this sense, we applied the last and more consistent method as shown in Fig D / panel C. The plot shows the division of 5 intervals according to the even oxygen d ecay.
To explore the minimum number of shells needed to obtain an acceptable accuracy, we explored the surface of oxygen diffusion. Results show that n=20 provides a sufficient accuracy (see Fig. D / panel D), at a tolerance of 0.3% in volume.

S2.4 Oxygen uptake rate
Oxygenation in tumours is often lower than in other organs (up to 6 times less [s8]).
Furthermore, hypoxic cells have little access to oxygen, for which metabolism slows down and to a pseudo-dormant state. In this first assumption, we consider 3 possible forms of the oxygen uptake rate, described by eqs. (S10)-(S12).
Zero Order 2 = ′ • 2 0 (S10) First Order 2 = ′ • 2 1 (S11) Second Order 2 = ′ • 2 2 (S12) Solving the problem described in eq. 6 of the main text with a given value of kR'= 2 cm -1 combined with the different expressions of the uptake rate described in eqs. (S10(S12), we obtain oxygen curves as described in Fig D /  give you the first order expression. In this case, we have no previous information to calculate the error function. Although the second order expression yields a similar oxygen distribution, it cannot be solved analytically. Also, the oxygen uptake has been reported to follow first order dynamics numerous times. These arguments drove us to the choice of the first order expression by eq. (S11).

S2.5 Proliferation constant as function of the Oxygen Diffusion
The proliferation constant should be oxygen dependent in our models (k p (PO2)). We

S3. Indentifiability problem
Parameter identifiability has been defined as: "Given a model of the system and specific input-output experiments, we ask, if the data were error free, could the parameters of the model be uniquely determined?" [36]. The concept of identifiability is directly linked to the model inference. In other words, we want to know whether our estimations really represent the model reality.
Two different kinds of numerical identifiability are commonly handled in the literature: structural and practical. We analysed our model from the perspective of both identifiability tests, prioritising the statistical confidence of our model inference.

S3.1 Structural Identifiability
As the name indicates, the structural identifiability depends on the intrinsic structure of the model and was first introduced in compartmental models by Bergman and Astrom Here we applied the Taylor's series problem to our observables (eq. 11 in main paper).

(S16)
Let us approximate VT using a Taylor series expanded around the initial tumour volume. (S19) For simplicity and due to the fact that the constants do not play any role in the structural indentifiability, we left lumped constants together. According to Taylor's theorem each coefficient is linearly independent from the rest and therefore uniquely identifiable. First and second coefficients (eqs. (S18) and (S19) show 2 linearily independent equations with two parameters. If we rearrange and combine the two equations in terms of k R', we obtain: To solve equation (S20), we face a similar challenge to solving equations of type: = • .
If the domain is real and positive volumes: ∀ , ; > 0 ∈ ℝ and also the parameter values > 0 ∈ ℝ, there is a unique solution confined in our physically feasible ranges. This makes the problem locally identifiable, which is a sufficient condition. However, special care has to be taken in the practical identifiability and significance of the parameters.

S3.2 Practical Identifiability
We refer to practical identifiability of nonlinear systems to the data-sensitive unique parameter estimations. There are a number of methods described in the literature, although . Values for κ and γ higher than a fixed threshold 1000 and 10 respectively was considered unidentifiable. We then can determine the number of elements of Σ which fulfil the above condition, i.e. the number of "identifiable parameters" nip. The remaining parameters are the "unidentifiable parameters" nup, which in that case would be fixed.

S3.3 Standard errors
The standard errors (SE) of the ODM has been calculated by means of the "Fisher information matrix", FIM: Above − ̅ represents the residual of the model estimation.
This method would provide some idea about the estimates, although it would be more accurately represented by bootstrapping or other resampling methods for better accuracy.

S4. Oxygen perfusion -observations
In this section we will report some key features related to oxygen diffusion observed in preclinical tumours and reported in the literature, supporting some of the assumptions made in the model.
Tumours are highly dynamic irregular complex structures of cells, but most importantly they are 3D structures. One interesting phenomenon to support this idea, as commented in the main text is that some necrotic areas originate in areas close to the vessels (CD31 positive, Fig E). This reminds us that immunohistochemical (IHC) data is data of a cross section of a 3D tumour, where vessels may be visibly occluded or apparently wide open but occluded somewhere else.   (Table A). Own observations 218 To elucidate the effect of vasculature in kR' we drew a cartoon to show that hypoxia appears at approximately 250-400 μm from non-vascularised tissue, triggering secretion of VEGF, PDGF and other angiogenesis promoters. This is well known as the "angiogenic switch".
Thereafter, as the tumour grows it develops its own vasculature which is more distributed and angiogenic the lower this parameter is (Fig F / panel A). Reduction of the oxygen uptake rate translates into an apparent longer oxygen penetration in the tissue (Fig F / panel B) and thus greater proliferation rates.