 Research article
 Open Access
 Published:
Multiscalar cellular automaton simulates invivo tumourstroma patterns calibrated from invitro assay data
BMC Medical Informatics and Decision Making volume 17, Article number: 70 (2017)
Abstract
Background
The tumour stroma or tumour microenvironment is an important constituent of solid cancers and it is thought to be one of the main obstacles to quantitative translation of drug activity between the preclinical and clinical phases of drug development. The tumourstroma relationship has been described as being both pro and antitumour in multiple studies. However, the causality of this complex biological relationship between the tumour and stroma has not yet been explored in a quantitative manner in complex tumour morphologies.
Methods
To understand how these stromal and microenvironmental factors contribute to tumour physiology and how oxygen distributes within them, we have developed a latticebased multiscalar cellular automaton model. This model uses principles of cytokine and oxygen diffusion as well as cell motility and plasticity to describe tumourstroma landscapes. Furthermore, to calibrate the model, we propose an innovative modelling platform to extract model parameters from multiple invitro assays. This platform provides a novel way to extract metadata that can be used to complement invivo studies and can be further applied in other contexts.
Results
Here we show the necessity of the tumourstroma opposing relationship for the model simulations to successfully describe the invivo stromal patterns of the human lung cancer cell lines Calu3 and Calu6, as models of clinical and preclinical tumourstromal topologies. This is especially relevant to drugs that target the tumour microenvironment, such as antiangiogenics, compounds targeting the hedgehog pathway or immune checkpoint inhibitors, and is potentially a key platform to understand the mechanistic drivers for these drugs.
Conclusion
The tumourstroma automaton model presented here enables the interpretation of complex invitro data and uses it to parametrise a model for invivo tumourstromal relationships.
Background
Tumours function as complex, interactive, multicellular aberrant organs. Within the tumour microenvironment, malignant cells exist amongst extracellular matrix, immune cells, vasculature, lymphatic vessels, fibroblasts, and other stromal cells. This tumour microenvironment can exert pro and/or antitumorigenic actions, depending on context, while malignant cells create a permissive and supportive tumour matrix by secreting stromamodulating growth factors, including PDGF, EGFR, VEGF, and TGFβ [1]. These factors can activate the surrounding stroma, causing the secretion of additional soluble molecules that mediate extensive crosstalk between the tumour cells and stromal components [2].
Stromal composition, architecture, and quantity varies between patient tumour types, while tumour heterogeneity is a hallmark of patient tumours. It is now wellrecognised that complex tumourstromal interactions underlie tumour growth, progression, and invasion. Recently, it has been demonstrated that the tumour microenvironment can affect tumour pathophysiology [3], therapeutic sensitivity, and response [4, 5]. The recognition of the importance of the stromal compartment in treatment outcome has led to the development of high throughput screens that incorporate a stromal component in addition to a tumour cell component for finding novel therapies. As such, it is necessary to develop further insight into tumourstroma interactions to aid in the modelling of preclinical drugresponse data in patient tumours [6,7,8]. Furthermore, tumour stromal morphology at tissue level was identified as a potential driver of drug efficacy in patientderived xenografts [9], also identified as determinant of tumour grade in prostate cancers [10]. The origination of these patterns (described in [11]) is in great measure the motivation of this work.
Despite the important role the tumour microenvironment, few preclinical tumour models include an extensive desmoplastic stroma or three dimensional interaction [12]. Here is when simulation work becomes a key translational tool. Recent work in the field of modelling has explored the interactions between tumourimmune cells and nutrients in a multiscalar manner [13,14,15], also incorporating an indepth study of intracellular mechanisms [16] and oxygen [17], as well as the role of peripheral anatomical features in disease progression, such as ductal niches [18]. However, to our knowledge, there has been no attempt to predict invivo tumourstroma growth and progression using merely invitro assays (see Fig. 1), we consider this a step forward in translational science. For simplicity, we describe two models: 1) the parameter estimation model (PEM) is used to calibrate 2) the tumourstroma model (TSM), which is a cellular automaton.
Methods
Tumourstroma model (TSM)
The TSM is a predictive multiscale latticebased cellular automaton model described by eight parameters (Figs. 1 and 2b) whose values can be estimated with the parameter estimation model (PEM) defined below. The TSM was inspired by a previously developed oxygencentred model [19, 20] and the subsequent discovery that tumourstroma morphology affects therapeutic response [9]. The model considers two cell types: stroma and tumour, whereby the latter can exist in three states: viable, hypoxic, or necrotic. These states depend on the cells’ access to oxygen and a time component (Fig. 2a).
The model is built on a fixed 2D regular lattice in which each square element (or voxel) represents groups in the range of thousands of cells of one kind. These cells, based on their position in space, behave metabolically different. Further, the model considers delays inherent to the biological processes, such as cell cycle progression, which are orchestrated by an internal clock defined by the variable tau (τ) in the time domain.
Our virtual tumour is composed of four types of cells which behave in four ways:

tumour cells may proliferate at a rate proportional to their access to oxygen and nutrients modelled as a characteristic doubling time (\( {\upalpha}_{\mathrm{T}} \)) and never until reaching a minimum time to divide (\( \uptau >{\upbeta}_{\mathrm{T}} \)) [20] (see equation (11)),

tumour cells may undergo reversible hypoxia if their oxygen concentrations sink below a threshold (\( {\mathrm{h}}_{\mathrm{H}} \)) for long enough (\( \uptau >{\upbeta}_{\mathrm{H}} \)) (see equation (14))

hypoxic cells may irreversibly become necrotic cells if they are hypoxic for long enough (\( \uptau >{\upbeta}_{\mathrm{N}} \)) (see equation (15)),

stromal cells may migrate in a random walk manner (controlled by \( {\upmu}_{\mathrm{S}} \)) or may be recruited close to the tumour (regulated by \( {\mathrm{k}}_{\mathrm{S}} \)), see equations (12)(13).
Up until now the relationship between tumour and stromal tissue assumed in the model is cooperativeonly. However, the stroma confers stiffness to the tumour matrix, thereby constraining growth. Since the model is built on a fixed lattice, there is a natural impossibility to translocate or divide if there is no empty space around and so, this constraint is built into the model, where h is the sum of empty spaces in the Moore neighbourhood of the cell in question (see equation (11) and (13)).
The model equations were written assuming that,

1)
the tumour mass is composed of tumour tissue (in three states: viable, hypoxic, and necrotic) and stroma (containing fibroblasts, immune infiltrate, and vasculature);

2)
all processes contain pure time delays, i.e. the cells have to spend a minimum time in the cell cycle before dividing; we avoid hereby the nonrealistic instantaneous cell division;

3)
hypoxia and necrosis appear sequentially following increasing exposure to oxygen deprivation;

4)
hypoxia is reversible, providing oxygen levels are recovered;

5)
stromal tissue has a purely synergistic biochemical relationship with tumour tissue, supplying oxygen and nutrients; and

6)
the high stiffness of the stroma may limit tumour progression.
These simple assumptions are enough to demonstrate the morphological and biomechanical roles of the stroma in different tumour phenotypes.
Definition of oxygen distribution
Define now the simplest problem in which within a tumour space V, there is one single element of stroma carrying functional vasculature that serves as a source of oxygen located at the point x _{so}{x_{so}, y_{so}} ⊂ V (see Fig. 2a). The equations of oxygen diffusion are based on the Fick’s law of molecular diffusion,
where, \( {\mathrm{P}}_{\mathrm{O}2} \) is the oxygen partial pressure and \( \mathrm{D} \) and \( {\mathrm{k}}_{\mathrm{R}} \) are the diffusion and oxygen uptake rate by the cells (Rate R) and \( t \) is a time variable.
Assuming first that the diffusion process occurs much more rapidly than the accumulation of HIF1α and the structural changes at tissue level, we can approximate Eq. (1) to a quasisteadystate \( \left(\frac{\mathrm{d}{\mathrm{P}}_{\mathrm{O}2}}{\mathrm{d}\mathrm{t}}\approx 0\right) \). With this and assuming that the diffusion coefficient is independent of the location in the tissue \( \left(\mathrm{D}\ne \mathrm{f}\left(\mathrm{space}\right)\right) \), we obtain,
Define P_{O2}(x _{ s }) ∈ ℝ^{2} as the map of oxygen computed for the first single point (x _{s} ∈ V). Now let us define the boundary conditions as (1) Diritchlet boundary at the stromal cell, where the concentration of oxygen should match the oxygen concentration in blood giving the constant P_{O2}(x _{s}) = P_{O2} ^{blood} and (2) Neumann bounds at the stroma cell, where there is a maximum (∇P_{O2}(x _{s}) = 0), and oxygen decreases gradually until infinity (\( \nabla {\mathrm{P}}_{\mathrm{O}2}\left(\infty \right)=0 \)).
There is no analytic closedform solution of the Eq. (2) for the given boundary conditions when the number of stromal cells grow. However, a reasonable equation contained in the solution space should contain the terms,
defined in the whole of V and where  x − x _{s} represents the Euclidean norm to the stromal cell in question. Now simply for the boundary conditions to be true, the expression \( {\mathrm{C}}_3=0 \) has to hold. This simplifies things significantly giving the expression
for one single stromal cell. Applying the boundary conditions described before and assuming that the constant is a function of the uptake rate we obtain \( {\mathrm{C}}_2=\sqrt{\frac{{\mathrm{k}}_{\mathrm{R}}}{9\mathrm{D}}} \) and \( {c}_1={P}_{O2}^{blood} \) for a single stromal cell.
However, in our reality we consider a growing number of stromal cells, whereby the analytic solution of this is not generalisable. To extend this approach to the multisource distribution, unlike for linear systems, the superposition principle cannot be applied. In this case, we will use a convenient approach in which the total equivalent distance will be calculated as the reciprocal sum of the distance to each stromal cell:
Despite this being just a merely convenient approach, the approximation is valid to the degree of accuracy required as demonstrated in Additional file 1.
Lastly, we ought to find the function \( {\mathrm{C}}_{2,\mathrm{eq}}=\mathrm{f}\left(\frac{{\mathrm{k}}_{\mathrm{R}}}{\mathrm{D}}\right) \) for the multiple source case, which we obtained empirically as
with \( {\mathrm{f}}_{\mathrm{corr}}\in \mathbb{R} \) and fixed value of \( {\mathrm{f}}_{\mathrm{corr}}=400 \) obtained by minimising the disparity between the analytic solution produced by Eq. (5) and the numeric solution computed with a GaussSeidel algorithm as further explained in Additional file 1. The Eqs. (4), (5), and (6) summarise in
which was taken as the final expression used for the simulations of this model, with k_{R}'=√k_{R}/D as the apparent oxygen uptake rate.
Definition of tumour growth
For the tumour cells, the oxygen uptake is directly proportional to the division rate. The proliferation constant is represented by a time constant \( {\upalpha}_{\mathrm{T}} \) and a pure time delay \( {\upbeta}_{\mathrm{T}} \), therefore the proliferation rate can be expressed as
Now substituting the oxygen levels on (7) with the Eq. (8) and integrating between \( {\displaystyle \int \frac{\mathrm{dT}\left(\uptau \right)}{\mathrm{T}\left(\uptau {\upbeta}_{\mathrm{T}}\right)}={\displaystyle {\int}_{\uptau ={\upbeta}_{\mathrm{T}}}^{\uptau}\frac{{\mathrm{P}}_{\mathrm{O}2}(x)}{\upalpha_{\mathrm{T}}}\mathrm{d}\uptau}} \) we obtain
where \( {\upbeta}_{\mathrm{T}} \) is the time at which the probability for transition is nonzero. This will be represented with a step function \( \mathcal{H}\left(\uptau {\upbeta}_{\mathrm{T}}\right) \) with \( \mathcal{H} \) being a one dimensional Heaviside function. Now since the logical expression: \( \mathrm{T}\left(\uptau \right)> T\left({\upbeta}_{\mathrm{T}}\right)\forall \uptau >{\upbeta}_{\mathrm{T}} \) holds, we invert the Eq. (9) to determine the fraction of cells that divided at any given time compared to the initial ones,
where \( {\mathrm{X}}_{\mathrm{T}\to 2\mathrm{T}}\left(\uptau \right) \) symbolises the fraction of cells that divide.
Subsequently, if we take the continuous deterministic fraction of dividing cells just defined in Eq. (10) in discrete terms, it follows that for each single cell with oxygen level defined by its position, the probability for it to commit to the cell cycle is statistically proportional to the fraction of proliferating cells (\( \Pr \left(\mathrm{T}\to 2\mathrm{T}\right)\propto {\mathrm{X}}_{\mathrm{T}\to 2\mathrm{T}} \)). However, in this case this probability is truncated by the pure time delay (defined as a Heaviside step function) and the availability of space to divide, where the expression (10) substituting for Eq. (7) gives the final expression (11).
In turn, stromal cells may be recruited with probability (12) and move with probability (13). Lastly, tumour cells may become hypoxic and subsequently necrotic with probabilities governed by the age of cells and oxygen reach as described in (13) and (14).
The equations and a graphical display are summarised below and in Fig. 2, whereas the codes can be found in Additional file 2.
Summary of equations
Tumour  \( \Pr \left(\mathrm{T}\to 2\mathrm{T}\right)=\left(1,,{\mathrm{e}}^{{\mathrm{e}}^{\frac{20}{3}\cdot {\mathrm{k}}_{\mathrm{R}}^{{\textstyle \hbox{'}}}\cdot {\left\left x{x}_{\mathrm{s}}\right\right}_{\mathrm{e}\mathrm{q}}}\cdot \left(\tau {\beta}_{\mathrm{T}}\right)/\left({\upalpha}_{\mathrm{T}}\right)}\right)\cdot \mathcal{H}\left(\tau {\beta}_{\mathrm{T}}\right),\kern0.5em \mathrm{f}\mathrm{o}\mathrm{r}\mathrm{h}\ne 0, \)  (11) 
Stroma  Recruitment \( \mathrm{\triangle}\mathrm{S}\left(\mathrm{t}\right)={\mathrm{k}}_{\mathrm{S}}\cdot \mathrm{T}\left(\mathrm{t}\right)\cdot \mathrm{\triangle}\mathrm{t} \),  (12) 
Motility \( \Pr \left({\mathrm{S}}_{0\to 1}\right)={\mu}_{\mathrm{S}}\cdot \Delta \mathrm{t}\cdot \mathrm{S}\left(\uptau \right)/\mathrm{T}\left(\uptau \right) \), for \( \mathrm{h}\ne 0 \),  (13)  
Hypoxia  \( \Pr \left(\mathrm{T}\to \mathrm{H}\right) = \mathcal{H}\left(\uptau {\upbeta}_{\mathrm{H}}\right)\cdot \mathcal{H}\left({\mathrm{P}}_{\mathrm{O}2}{\mathrm{h}}_{\mathrm{H}}\right) \)  (14) 
Necrosis  \( \Pr \left(\mathrm{H}\to \mathrm{N}\right) = \mathcal{H}\left(\uptau {\upbeta}_{\mathrm{N}}\right) \)  (15) 
The TSM was run stochastically at a time step of \( \Delta \mathrm{t}=3\mathrm{h} \) on a square 2D lattice of 50×50 voxels. The initial conditions are described in Additional file 1. Once the probabilities have been calculated, the algorithm should check the events one by one, comparing a random number to their probability. The steps can be summarised as:

1.
selection of time step \( \Delta \mathrm{t} \),

2.
update time,

3.
randomly select each possible event and compare the probability of each event \( {\mathrm{f}}_{\mathrm{e}}\left(\uptau, \mathrm{r}\right) \) to a random number,

4.
update the state,

5.
repeat 3–4 until all possible events have been checked,

6.
repeat from 2 until simulation time is reached.
Parameter estimation model
The PEM transforms both 2Dimensional cell culture and 3Dimensional tissue slice culture biomarker expressions into the invivo parameter values of the TSM using a number of mathematical strategies graphically summarised in Fig. 2a. Briefly, the PEM is subdivided into five modules: growth curve module, immunohistochemistry (IHC) image processing module, heterogeneity module, Western blot module, and the oxygendiffusion module. A mathematically formal description of these modules is provided below.
Growth curve module
Tumour parameters \( {\upalpha}_{\mathrm{T}} \) and \( {\upbeta}_{\mathrm{T}} \) will be calculated from invitro cell confluence assays on the cell line MCF7. We use the mean of the populations of cells as the expected value for a quasiexponential cell growth with constant \( 1/{\upalpha}_{\mathrm{T}} \) (see Fig. 3a, fit) and the standard deviation as a measure of the diversity in the delay for the cells to commit to the cell cycle (\( {\upbeta}_{\mathrm{T}} \)). Let \( {\mathrm{C}}_{\mathrm{n}}\left(\mathrm{t}\right) \) be the % of confluence, \( \mathrm{t} \) the time measured, and \( {\upsigma}_{\mathrm{Cn}}\left(\mathrm{t}\right) \) the standard deviation. We can then define our parameters as,
where the linear matrix product in Eq. (16) was solved by Gaussian elimination. All units are in hours.
The results for both parameters for the cell line MCF7 are \( {\upalpha}_{\mathrm{T}}=9.23\mathrm{h} \) and \( {\upbeta}_{\mathrm{T}}=18.17\mathrm{h} \), which together sum up to 27.4 h. Surprisingly, this value is very close to the doubling times reported in the literature (20–40 h [21]) and somewhat higher to the doubling time calculated from the same data set (20.64 h).
Immunohistochemistry (IHC) module
We now concentrate on the spatial distribution of hypoxia biomarkers to estimate oxygen diffusion parameters. IHC was carried out on paraffin embedded tissue slice sections as has been previously described in Davies et al. [22]. To calculate delays in the onsets of hypoxia and necrosis we used tissue culture experiments on thirteen different slices, which after isolation, sectioning, and staining we quantified digitally and discretised evenly with the socalled Tissue Culture Profiling (TCP) algorithm (as described in the Additional file 1), giving the results shown in Fig. 3c. Judging by the mismatch between the two curves at 24 h and 48 h (see Fig. 3c) there is a delay in the onset of HIF1α, also previously described as inherent to the cell biology [23]. We calculated these delays in the onset of hypoxia with the socalled thresholdbased method which utilises a binary onset of hypoxia at a certain absolute oxygen level threshold value, although other methods have been explored as summarised in [11].
This method utilises the absolute curve shift between the curves at 24 and 48 h for 20% expression of HIF1α and 10% of necrosis as the standards to calculate the delays. The parameters can be calculated with
where \( {\upbeta}_{\mathrm{H}} \) and \( {\upbeta}_{\mathrm{N}} \) are the delays of hypoxia and necrosis for the thresholdbased method. Their values are \( {\upbeta}_{\mathrm{H}}=139\mathrm{h} \) and \( {\upbeta}_{\mathrm{N}}=184\mathrm{h} \), which correspond to 43 μm and 32 μm of curve shift respectively.
Heterogeneity module
The heterogeneity module uses imaging data to output a measure of stromal heterogeneity, which we define based on the fractal dimension of the 2D image. We estimate fractal dimensions via a countingbox algorithm, which counts boxes of decreasing sizes needed to cover the region of interest of our image. The total count of boxes (\( \epsilon \)) is then fitted to a line in a doublelogarithmic scale and extrapolating it to size zero. Thus,
From the fractal dimensions, a single value of the heterogeneity has been taken into account (assuming for a 2D image), where the heterogeneity is a deviation from the “integer” dimension, thus
where \( {\Psi}_F \) represents the heterogeneity calculated from the fractal dimension.
Finally, the stromal parameters are described directly by the total quantification and by the heterogeneity factor, as
where \( {\Psi}_{\mathrm{F}} \) is again the fractal heterogeneity and S/T is the stroma to tumour ratio, usually expressed as a percentage. Note that T refers only to the tumour fraction and not the total, hence \( \mathrm{T}=\left(\mathrm{Total}\mathrm{S}\right) \). The value of 10 is a convenient number added to counterbalance the error introduced by the analysis of a crosssection in lieu of the whole tissue culture invivo.
Now accounting for the available data we obtain,
for both parameters respectively which results in values of 0.064 h^{−1} (Str/T)^{−1} and 0.0031 h^{−1} (Str/T)^{−1}. This calculation has been summarised in Fig. 3b.
Western blot module
To find the oxygen concentration at which HIF1α is expressed, the MCF7 cell line was cultivated invitro at different oxygen concentrations (1, 5, and 21%) and then Western blots [24] were run for HIF1α. Standard Western blot protocols were used, HIF1α (#610959, BD Biosciences, 1:250 dilution). The intensity of the chemoluminescense is proportional to the amount of HIF1α in the sample (see Fig. 3d). The experiments were done on two separate occasions: first 1 and 21% O_{2} atmospheres were used, then 5 and 21% O_{2} with a total of 12 measurable points. Since the chemiluminescent exposure may have been slightly different, the sample results have been linearly normalised for the reference average of 21% O_{2}.
The results show a fair exponential distribution (as expected) with formula
for which an arbitrary expression of 20% of HIF1α would need around 0.064 mmol/L of oxygen in tissue or around 5% of atmospheric oxygen in culture (see Fig. 3d). The estimation has been done with a naïvepooled algorithm with the 12 points depicted in Fig. 3d.
Oxygen diffusion module
In this module, we attempt to calculate oxygen diffusion parameters as described in Eq. (2) by means of a tissue slice experiment. The geometry of the slice leads to the one dimensional resolution of Eq. (2) solved in Cartesian coordinates, with Dirichlet boundary conditions for both sides of \( {\mathrm{P}}_{\mathrm{O}2}\left({\mathrm{x}}_{\mathrm{filter}}\right)=0.28\;\mathrm{mmol}/\mathrm{L} \) and \( {\mathrm{P}}_{\mathrm{O}2}\left({\mathrm{x}}_{\mathrm{air}}\right)= \) \( 1.19\;\mathrm{mmol}/\mathrm{L} \), being these calculated from themodynamic equilibrium expressions (see section Additional file 1). These quantities are comparable to the cava return veins and alveolar concentrations [25]. The approximate analytic expression results in
whereby \( {\mathrm{k}}_{\mathrm{R}}^{\prime }=\sqrt{{\mathrm{k}}_{\mathrm{R}}/\mathrm{D}} \) and \( \mathrm{x},{\mathrm{x}}_{\mathrm{air}} \) are the 1D positions in the transversal direction of the slice. Note that in this case the origin of coordinates was set at the filter (\( {\mathrm{x}}_{\mathrm{filter}}=0 \)). A two dimensional version of the algorithm was also discussed in [11].
Summary of results of the PEM
The complete parameter set has been estimated for the cell line MCF7 only, which results are presented in Table 1. To expand the use of the model and test the formation of complex stromal patterns, we decided to explore the cell lines Calu3 and Calu6, which demonstrate very different stromal morphological features between them. However, data availability is limited for the last two cell lines. To solve this problem, the values of parameters \( {\beta}_H,{\beta}_N,{k}_R^{\prime },\mathrm{and}\;{h}_H \) were assumed to be the same as for MCF7, whereas the values of \( {\upalpha}_{\mathrm{T}},{\upbeta}_{\mathrm{T}},{\mathrm{k}}_{\mathrm{S}},\mathrm{and}\;{\upmu}_{\mathrm{S}} \) were estimated using the ratio \( {\upalpha}_{\mathrm{T}}/{\upbeta}_{\mathrm{T}} \) unchanged but corrected for the doubling times of the cell lines (35 h for Calu3 and 26 h for Calu6 [26]) and the whole tumour histological images to correct for the stromal recruitment. For the codes and examples on the PEM see Additional file 3.
Results & discussion
Once the algorithms have been presented and the parameters estimated, it follows to analyse the potential of their tumourstromal model in the different phenotypes.
Simulations on complete parameter set from parameter estimation model
Using the parameters obtained from the invitro platforms for MCF7 summarised in Table 1, we have made simulations to build Fig. 4a. In this figure, we observe that the necrotic and hypoxic regions appear concentrically as it does in reality (see IHC image on Fig. 4a and Ref. [27]); the oxygen levels are very low; and the viable rim is of almost constant thickness, which is consistent with numerous observations [28]. All these properties are congruent with the literature and our data; and are futher equivalent to those of an avascular tumour.
However, the values reported for MCF7 obtained from the fitting of other invivo models [19] resulted to be \( {k}_R^{\prime }=7.54 c{m}^{1} \) which is 24.5 times smaller than that reported in the here estimated invitro system. This –as reported in by the authors in Ref. [19]  is due to the increased oxygenation proportionated by the angiogenenic vessels. If we now simulate the same system changing the value of \( {k}_R^{\prime } \), we observe that there is no apparent difference in the distribution of stroma. However, the viability of the tissue has changed, as the whole tissue is well perfused and viable (Fig. 4a). The actual linear oxygen reach now would have increased to 6.5 mm. This does not reflect the real tissue distribution as shown by the IHC (Fig. 4f) but will be useful to emulate a number of clinically possible scenarios as we will see later on.
The model predictions for \( {k}_R^{\prime } \) value of \( 185 c{m}^{1} \) fail to represent the growth rate accurately, although the estimations are within a close range to the data (see Fig. 4d), whereas for the value of \( 7.54 c{m}^{1} \) the predicted growth mimics almost perfectly the growth rate (see Fig. 4d). The overall growth rates differ by 82.4% at the end of the study, which is a remarkable difference (Fig. 4cd). The moderate stromal growth (~25% at 0 days to 17.5% and 11.7% at 30 days) is slightly higher than that of the real observations, where the proportion of the stroma at the end of the study lies around \( 4.1\%\pm 3.4\% \) at day 30 (n=8, Fig. 4f). It is important to point out that those values were calculated by manually outlining the visual stromal regions on images stained for HIF1α. This means that the values may (and probably will) vary greatly when stained for activated stroma (αSMA). However, due to the lack of time and staff those experiments were not conducted. The predicted proportions of hypoxia (24.3%) and necrosis (48.6%) are likewise within the ranges of the observed values with means of \( 36.8\%\pm 7.1\% \) (n=10) and \( 36.7\%\pm 14.1\% \) (n=8) respectively (Fig. 4f).
It is interesting that the modification of the oxygen uptake rate alone is sufficient to simulate vastly different possible clinical scenarios, however none of those is completely congruent with both the growth rate and stromal distribution phenotypes (see Fig. 4). On the one hand, a low \( {k}_R^{\prime } \) will provide with high rates of proliferation to the cells with space to proliferate, which is consistent with the overall growth rates observed. On the other hand, we see a large necrotic core forming in the tumour, which is indicative that \( {k}_R^{\prime } \) is actually a function of the tumour position with radial geometry.
Modelling the different phenotypes: Calu3 and Calu6
The parameters estimated for Calu3 and Calu6 (Table 1) were used in the equations (11)(15) to simulate the results shown in Fig. 5. Firstly, Calu6 (Fig. 5b) shows a very similar profile as the one shown on the MCF7 cell line (owing to the similar parameter values, compare Table 1), which again is congruent with reality, as Calu6 and MCF7 present with similar histological profiles.
However, when we look at Calu3 the situation changes. Note that the only parameters that differ are the proliferation rates \( {\upalpha}_{\mathrm{T}} \) and \( {\upbeta}_{\mathrm{T}} \), and the stromal recruitment \( {\mathrm{k}}_{\mathrm{S}} \) (Table 1). Now, if we look at the oxygen maps, these have changed significantly. The growing representation of the stromal tissue raises oxygen levels keeping the tissue viable. We observe this in the increased viable rim on Fig. 5ab. The heterogeneity of Calu3 seems to be also lower, breaking with the concentric symmetry observed in the other two cell lines Calu6 and MCF7.
After validation of the patterns observed for Calu3 and Calu6 it is worth examining the volumetric growth rates predicted by the model (see Fig. 5c). The key observations here are that the proportion of stroma in Calu3 is significantly higher; at 15 days the hypoxic fraction has disappeared in Calu3 while it grows steadily for Calu6; and consequently the necrosis fraction reaches a plateau for Calu3 whereas for Calu6 it continues expanding (Fig. 5ab).
Tumour volumes calculated with the TSM are higher than the data (Fig. 5c), although within similar ranges. Most importantly, the global growth speed of either invivo model is incongruent with the data, because the simulated Calu3 grows faster than simulated Calu6, which is in conflict with empirical observations.
In terms of stroma, when Calu3 and Calu6 tumour xenografts are analysed by IHC we observe that not only the amount of stroma is significantly higher in Calu3, but also the Mean Vascular Density (MVD), which negatively correlates to the amount of necrosis in these images, see Additional file 1). The quality of the vessels observed in average was similar in both models with no significant differences in vessel thickness, vascular area, lumen area or vessel perimeter (see Additional file 1).
The simulations produced by the original parameter set (\( {\mathrm{k}}_{\mathrm{R}}^{\prime }=185\mathrm{c}{\mathrm{m}}^{1} \)) fail to reproduce the heterogeneity at the end of the study for MCF7, whereas for the adapted parameter set (\( {\mathrm{k}}_{\mathrm{R}}^{\prime }=7.54\mathrm{c}{\mathrm{m}}^{1} \)) the heterogeneity values are very close to the data (Fig. 4e). Moreover, the values of heterogeneity observed in the simulations for Calu3 and Calu6 mirrors almost perfectly the data (Fig. 5d). All trends are downwards with the progression of time, which is consistent with the visual inspection of the images and the general rise in entropy of the distribution of stromal cells (Fig. 5ab and Fig. 4ab). Interestingly, the evolution of heterogeneity for the adapted parameter set for MCF7 corresponds to the one observed for Calu3.
Simulations with parameter sets for both cell lines show remarkable similarities between the tumour, stroma, necrotic, and hypoxic cell distributions and the tumour histopathology observed for Calu6 (Fig. 5a) and Calu3 (Fig. 5b) when grown invivo. It is notable that the viable rim in each simulated phenotype varies significantly, as demonstrated by the oxygen distribution maps (Fig. 5ab, top row). However, results on the total tumour volume show that the growth dynamics of the cell line Calu6 are well captured (Fig. 5c), whereas the cell line Calu3 appears to grow faster in the model simulations than in reality (Fig. 5c). This effect reveals that the tumourstroma relationship might not only be synergistic by providing growth factors and nutrients to the tumour, but it might actually constraint its growth mechanically (fibrous stroma confers rigid extracellular matrices with lesser room for cell proliferation [29, 30]), and chemically (cytotoxic components of the immune infiltrate play an active role in the retardation of tumour growth [31]).
Lastly, a comprehensive analysis of various simulation settings in terms of time step, lattice size, or lattice dimensionality (3D) will be an interesting exercise which will provide with other insights of the TSM. Especially interesting will be the usage of a cell resolution lattice and a sufficiently short time step to properly capture diffusion of macromolecules and changes in cell polarity. These computations require high specification hardware and computation time.
Conclusions
The PEM was a useful platform to generate invivorelevant parameter values from invitro data; nevertheless its usage relies upon the availability of the invitro assays mentioned above. To circumvent this, we proposed a feasible parameter extrapolation method which resulted valid and may be useful to other scientists. The results generated by the PEMcalibrated TSM for the two cell lines Calu3 and Calu6 demonstrate that the tumourstroma mutuallypositive relationship assumed here does not represent both the stromal patterns and growth rates simultaneously. This means that while the stroma constrains growth by limiting space to grow, these effects do not explain the slower growth rate of Calu3 (Fig. 5c). Experimentally we observe an increase in thickness of the viable tumour rim of the Calu3 model, which is explained by the exacerbated blood vesselcarrying stromal recruitment described with the TSM. Overall, these models suggest that at a tissue level, the stroma limits tumour growth, and at a molecular level, the stroma fosters tumour viability. Empirically, it has been shown invitro that stromal components provide growthlimiting molecular signals in the form of cellsequestering factors (not built into our model) [1], which may explain the disparity in predicted and observed tumour growth.
Our invitrocalibrated mathematical models describe aspects of invivo tumour biology extrinsic to the tumour cell, which may provide a valuable addition to existing predictive models, such as cancer pharmacokineticpharmacodynamicefficacy models developed using cell linederived xenograft models. This adds another dimension to simplistic tumour models which may allow us to better predict the effect of drugs in morphologicallycomplex clinicallyrelevant pathophysiologies [3, 9].
Change history
23 June 2017
An erratum to this article has been published.
Abbreviations
 H&E:

Hematoxylin and eosin
 HIF1α:

Hypoxia inducible factor 1α
 IHC:

Immunohistochemistry
 MVD:

Mean vascular density
 PEM:

Parameter estimation model
 SM:

Supplementary materials
 TCP:

Tissue culture profiling algorithm
 TSM:

Tumourstroma model
References
 1.
Mueller MM, Fusenig NE. Friends or foes—bipolar effects of the tumour stroma in cancer. Nat Rev Cancer. 2004;4(11):839–49.
 2.
Whipple CA. Tumor talk: understanding the conversation between the tumor and its microenvironment. Cancer Cell Microenviron. 2015;2(2):e773.
 3.
Smith NR, Baker D, Farren M, Pommier A, Swann R, Wang X, Mistry S, McDaid K, Kendrew J, Womack C. Tumor stromal architecture can define the intrinsic tumor response to VEGFtargeted therapy. Clin Cancer Res. 2013;19(24):6943–56.
 4.
McMillin DW, Negri JM, Mitsiades CS. The role of tumour–stromal interactions in modifying drug response: challenges and opportunities. Nat Rev Drug Discov. 2013;12(3):217–28.
 5.
Rampias T, Favicchio R, Stebbing J, Giamas G. Targeting tumor–stroma crosstalk: the example of the NT157 inhibitor. Oncogene. 2015;35(20):25624.
 6.
Knowlton S, Joshi A, Yenilmez B, Ozbolat IT, Chua CK, Khademhosseini A, Tasoglu S. Advancing cancer research using bioprinting for tumoronachip platforms. Int J Bioprinting. 2016;2(2):3–8.
 7.
Das V, Bruzzese F, Konečný P, Iannelli F, Budillon A, Hajdúch M. Pathophysiologically relevant in vitro tumor models for drug screening. Drug Discov Today. 2015;20(7):848–55.
 8.
Choi SYC, Lin D, Gout PW, Collins CC, Xu Y, Wang Y. Lessons from patientderived xenografts for better in vitro modeling of human cancer. Adv Drug Deliv Rev. 2014;79:222–37.
 9.
Delgado San Martin J, Hare J, Yates J, Barry S. Tumour stromal morphology impacts nanomedicine cytotoxicity in patientderived xenografts. Nanomedicine. 2015;11(5):1247–52.
 10.
Pérez JL, Puente ET, Kulich EI, MAr JAB, Nistal M, Peramato PG, García MR, Villar JMN, De Miguel M. Relationship between tumor grade and geometrical complexity in prostate cancer. bioRxiv. 2015:015016.
 11.
Delgado San Martin JA. Mathematical models for preclinical heterogeneous cancers. http://ethos.bl.uk/OrderDetails.do?uin=uk.bl.ethos.690595: University of Aberdeen; 2016. Accessed May 2017.
 12.
Unger C, Kramer N, Walzl A, Scherzer M, Hengstschläger M, Dolznig H. Modeling human carcinomas: physiologically relevant 3D models to improve anticancer drug development. Adv Drug Deliv Rev. 2014;79:50–67.
 13.
RobertsonTessi M, Gillies RJ, Gatenby RA, Anderson AR. Nonlinear tumorimmune interactions arising from spatial metabolic heterogeneity. bioRxiv. 2016:038273.
 14.
Liang S, Hoskins M, Khanna P, Kunz RF, Dong C. Effects of the tumorleukocyte microenvironment on melanoma–neutrophil adhesion to the endothelium in a shear flow. Cell Mol Bioeng. 2008;1(2–3):189–200.
 15.
Starruß J, de Back W, Brusch L, Deutsch A. Morpheus: a userfriendly modeling environment for multiscale and multicellular systems biology. Bioinformatics. 2014; 30(9):13312
 16.
Wang Z, Butner JD, Kerketta R, Cristini V, Deisboeck TS. Simulating cancer growth with multiscale agentbased modeling. Semin Cancer Biol. 2015;30:70–78.
 17.
Chakrabarti A, Verbridge S, Stroock AD, Fischbach C, Varner JD. Multiscale models of breast cancer progression. Ann Biomed Eng. 2012;40(11):2488–500.
 18.
Picco N, Gatenby R, Anderson A. Nichedriven stem cell plasticity and its role in cancer progression. bioRxiv. 2016:056762.
 19.
DelgadoSanMartin J, Hare J, Moura A, Yates J. Oxygendriven model: a pathologyrelevant model of tumour growth in xenografts. PLoS Comput Biol. 2015;11(10):e1004550.
 20.
Weber TS, Jaehnert I, Schichor C, OrGuil M, Carneiro J. Quantifying the length and variance of the eukaryotic cell cycle phases by a stochastic model and dual nucleoside pulse labelling. 2014.
 21.
Sutherland RL, Hall RE, Taylor IW. Cell proliferation kinetics of MCF7 human mammary carcinoma cells in culture and effects of tamoxifen on exponentially growing and plateauphase cells. Cancer Res. 1983;43(9):3998–4006.
 22.
Davies EJ, Dong M, Gutekunst M, Närhi K, van Zoggel HJ, Blom S, Nagaraj A, Metsalu T, Oswald E, ErkensSchulze S. Capturing complex tumour biology in vitro: histological and molecular characterisation of precision cut slices. Sci Rep. 2015;5:17187.
 23.
Russell J, Carlin S, Burke SA, Wen B, Yang KM, Ling CC. Immunohistochemical detection of changes in tumor hypoxia. Int J Radiat Oncol Biol Phys. 2009;73(4):1177–86.
 24.
Liu ZQ, Mahmood T, Yang PC. Western blot: technique, theory and trouble shooting. N Am J Med Sci. 2014;6(3):160.
 25.
Pittman RN. Regulation of tissue oxygenation. Colloquium series on integrated systems physiology: from molecule to function. 2011;3(3):1–100.
 26.
Brower M, Carney DN, Oie HK, Gazdar AF, Minna JD. Growth of cell lines and clinical specimens of human nonsmall cell lung cancer in a serumfree defined medium. Cancer Res. 1986;46(2):798–806.
 27.
Ribba B, Watkin E, Tod M, Girard P, Grenier E, You B, Giraudo E, Freyer G. A model of vascular tumour growth in mice combining longitudinal tumour size data with histological biomarkers. Eur J Cancer. 2011;47(3):479–90.
 28.
Evans ND, Dimelow RJ, Yates JW. Modelling of tumour growth and cytotoxic effect of docetaxel in xenografts. Comput Methods Prog Biomed. 2014;114(3):e3–13.
 29.
Netti PA, Berk DA, Swartz MA, Grodzinsky AJ, Jain RK. Role of extracellular matrix assembly in interstitial transport in solid tumors. Cancer Res. 2000;60(9):2497–503.
 30.
Van Liedekerke P, Palm M, Jagiella N, Drasdo D. Simulating tissue mechanics with agentbased models: concepts, perspectives and some novel results. Comput Part Mech. 2015;2(4):401–44.
 31.
Fridman WH, Pages F, SautesFridman C, Galon J. The immune contexture in human tumours: impact on clinical outcome. Nat Rev Cancer. 2012;12(4):298–306.
Acknowledgements
We thank the NICHE MarieCurie consortium and AstraZeneca for the funding and data presented. We acknowledge the work carried out within the PREDECT consortium for development of the tissue slice culture invitro system. Also, thank you to Dr Alessandro de Moura, Dr Mark Penney, and Dr Ekkehard Ullner for their advice in the mathematics and Mandy Lawson for providing invitro MCF7 data.
Funding
The work has been partly funded by the European Marie Curie ITN Action number FP7PEOPLE2011ITN289384 and AstraZeneca, UK.
Availability of data and materials
All codes can be found within the supplementary material of this publication. All data utilised and the derived calculations is presented in graphic from in the publication. The data that supports the findings of this study is available from AstraZeneca but restrictions apply to the availability of this data which were used under license for the current study and is not publicly available. Data are available from the authors upon reasonable request and with permission of AstraZeneca.
Authors’ contributions
JD & JY conceived and designed the mathematical models. JD carried out all the drylab work, whereas JH & ED completed the wetlab work. All authors contributed to writing the paper, read, and approved the final version of the paper.
Competing interests
The authors, funding bodies, and AstraZeneca disclose no competing interest or intellectual property.
Consent for publication
Not applicable.
Ethics approval
All animal studies were conducted in accordance with UK Home Office legislation, the Animal Scientific Procedures Act 1986 (ASPA) and with AstraZeneca Global Bioethics policy. The analysis in this paper is retrospective, utilising control/untreated animals of different oncology projects within AstraZeneca between 2004 and 2015. All tumour volumes, animal weights and welfare were maintained within the margins fixed by UK and European regulations. No data was generated specifically for this paper.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Author information
Affiliations
Corresponding author
Additional information
This article was updated and an erratum published to correct typesetting discrepancies in eqs. (13) and (14).
An erratum to this article is available at https://doi.org/10.1186/s1291101704927.
Additional files
Additional file 1:
Supplementary material: This file contains additional information on the experiments conducted and the raw results as well as algorithmic detail of the scripts that quantify the distribution of IHC stain within a tissue culture, and the calculation of the oxygen concentrations at either interface of the tissue slice. (PDF 927 kb)
Additional file 2:
mainTSM_Code: Main script for the TumourStroma Model (TSM). This script is a standalone script that simulates the model given certain parameter values and initial conditions. (CPP 9.34 kb)
Additional file 3:
PEM_Code: Compressed folder containing the codes for the Parameter Estimation Model algorithms and example figures. These scrips are standalone and contain all the necessary information to replicate the data in the manuscript. (7Z 8793 kb)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
DelgadoSanMartin, J.A., Hare, J.I., Davies, E.J. et al. Multiscalar cellular automaton simulates invivo tumourstroma patterns calibrated from invitro assay data. BMC Med Inform Decis Mak 17, 70 (2017). https://doi.org/10.1186/s1291101704611
Received:
Accepted:
Published:
Keywords
 Cancer
 Stroma
 Cellular automaton
 Oxygen
 Hypoxia
 Immunohistochemistry
 Cell line