Morphological stability for in silico models of avascular tumors

The landscape of computational modeling in cancer systems biology is diverse, offering a spectrum of models and frameworks, each with its own trade-offs and advantages. Ideally, models are meant to be useful in refining hypotheses, to sharpen experimental procedures and, in the longer run, even for applications in personalized medicine. One of the greatest challenges is to balance model realism and detail with experimental data to eventually produce useful data-driven models. We contribute to this quest by developing a transparent, highly parsimonious, first principles silico model of a growing avascular tumor. We initially formulate the physiological considerations and the specific model within a stochastic cell-based framework. We next formulate a corresponding mean-field model using partial differential equations which is amenable to mathematical analysis. Despite a few notable differences between the two models, we are in this way able to successfully detail the impact of all parameters in the stability of the growth process and on the eventual tumor fate of the stochastic model. This facilitates the deduction of Bayesian priors for a given situation, but also provides important insights into the underlying mechanism of tumor growth and progression. Although the resulting model framework is relatively simple and transparent, it can still reproduce the full range of known emergent behavior. We identify a novel model instability arising from nutrient starvation and we also discuss additional insight concerning possible model additions and the effects of those. Thanks to the framework's flexibility, such additions can be readily included whenever the relevant data become available.


Introduction
Tumors are highly complicated biological systems, yet constitute a concrete example of cellular self-organization processes amenable to modeling in silico [7].In a cancerous tumor, the cells have undergone several significant mutations and obtained distinct hallmarks providing the population with remarkable growth capabilities [24,25].Furthermore, populations comprising large numbers of cells interact on multiple scales, yielding a range of emergent phenomena [13], which can be studied using computational models based on knowledge of single-cell behavior.
To the modeler's aid in this regard, biological data streams nowadays contain detailed features at the individual cell level such as cell size and -type, mutation-and growth rate, molecular constituents, and gene expression [2,9,32].
Complementing biological experiments, mathematical models can in addition provide explanations to observed data, concerning, e.g., drug-response in tumor growth, with potential applications in precision medicine [4,40].Progress in cell biology has led to a good understanding of intra-cellular processes which unlocks the possibility to model these systems from fundamental principles, to translate 'word models' formulated from biological experiments into mathematical and computational models, and to test the features of these models [31].Often quoted uses of computational models include the testing of hypotheses, the investigation of causality, and the integration of knowledge when comparing in vitro and in vivo data [6].Bayesian inference methods present a means to quantitatively investigate these matters, provided there exist appropriate data and meaningful priors associated with the model parameters.
Several cell population models exist in the literature, ranging from continuous to agentbased to hybrid models, and taking place at various scales [10,14,17,34,36].Such models may reach a predictive power, where agreement/disagreement with biological data can advance our understanding of mechanistic relations within the biological systems [5,15,19,26].Pertinent to the present work, previous research shows how analyzing the emergent morphology of cell population models can provide insight into the role of the model parameters [1,20,21], promoting future use of Bayesian methods.Such analysis has, for example, enabled modelers to analyze the behavior of the invasive fronts of tumor models and their response to parameter changes representing vascularization, nutrient availability [11,27], and cell-cell adhesion effects [8,35].
Motivated in part by improvements of in vitro techniques for obtaining detailed time-series tumor and single-cell data and the current trend in computational science towards data-driven modeling, we present and analyze a basic continuous mathematical model of avascular tumor growth, here derived from a previously developed cell-based model [16].Our aim is that the model should be highly parsimonious in order to cope with issues of model identifiability.For this purpose the in silico tumor's fate should be well understood when regarded as a map from parameters to simulation end-result.Initial results from an earlier version of the model [16] display boundary instabilities, akin to those discussed in [22], which we analyze thoroughly.The self-regulating properties of avascular tumors concerning size that have been observed in vitro [18] and in silico [23] motivate a careful investigation into the model capabilities in this regard.
We have structured the paper as follows.In §2 we summarize the stochastic cell-based tumor model as well as the associated mean-field space-continuous version.We analyze the latter in §3, assuming radially symmetric solutions first, and then via linear stability analysis.In §4 we investigate via numerical examples the key aspects of the analysis as well as its relevance for the stochastic model.A concluding discussion is found in §5.

Stochastic modeling of avascular tumors
An advantage held by stochastic models is that they implicitly define a consistent likelihood and thus formally have the potential to be employed in Bayesian modeling when confronted with data.We summarize our stochastic framework in §2.1 and the basic stochastic tumor model in §2.2.We next derive a corresponding mean-field version in §2.3, which has the distinct advantage of being open to mathematical analyses.

Stochastic framework
The Darcy Law Cell Mechanics (DLCM) framework [16] is a cell-based stochastic modeling framework where the cells are explicitly represented and the rates of their state updates, e.g., movement, proliferation, death, etc., are determined and govern the corresponding events in a continuous-time Markov chain.Movements of the cells generally follow Darcy's law for fluid flow in a porous environment, but since the framework takes place in continuous time, other types of cell transport are easily incorporated.
The spatial domain is discretized into i = 1, 2, ..., N vox voxels v i , and populated by a total of N cells cells.The DLCM framework can be used over any grid for which a consistent discrete Laplace operator can be derived.Each voxel may be empty or contain some number of cells and if this number exceeds the voxel's carrying capacity, the cells will exert a pressure onto the cells in the surrounding voxels, see Fig. 2.1.The pressure propagates through the considered domain and the local pressure gradient induces a cell flow.The simplest implementation allows each voxel to be populated by u i ∈ {0, 1, 2} cells at any time, thus with a carrying capacity of 1.Note that carrying capacity here does not refer to the maximum possible number of cells in a voxel, but rather to the capacity beyond which a voxel is no longer in a mechanically relaxed state.
Although the state u i takes on discrete values in each voxel and at any given time, the governing model is derived from a continuous assumptions where the corresponding cell density is then u = u(x, t).We let u be governed by the continuity equation where I is the flux of u.There are three main assumptions in the DLCM framework, with the first assumption pointing to the central role of the flux I: Assumption 2.1.Consider the discrete tissue formed from the population of cells distributed over the grid.We assume that: 1.The tissue is in mechanical equilibrium when all cells are placed in a voxel of their own.
2. The cellular pressure of the tissue relaxes rapidly to equilibrium in comparison with any other mechanical processes of the system.
3. The cells in a voxel occupied by n cells may only move into a neighboring voxel if it is occupied by less than n cells.From Assumption 2.1(1) it is clear that random (e.g., Brownian) motion around a voxel center is ignored, and, therefore, that only voxels with cells above the carrying capacity are considered pressure sources.As in [16], the flux is determined from the pressure gradient in the form of Darcy's law which can be derived as a limit for flow through porous media [39]: where p is the pressure and the Darcy constant D can be interpreted as the ratio of the medium permeability κ to its dynamic viscosity µ, D := κ/µ.The relation between pressure and cell population is completed by assuming a constitutive relation in the form of a heat equation for the pressure and using Assumption 2.1(2) to arrive at the stationary relation where s(u) is the pressure source which equals to one for voxels above the carrying capacity, u i > 1, and zero otherwise.Finally, we detail the flux parameter D in (2.2).Let R(e) denote the rate of an event e, and let i → j denote the event that a cell moves from voxel v i to v j .With unit carrying capacity only two distinct movements rates are possible according to Assumption 2.1(3): one for cells moving into an empty voxel, and one for cells moving into an already occupied voxel, where D 1 and D 2 are (possibly equal) conversion factors from units of pressure gradient to movement rate for the respective case.Here, I(i → j) is the pressure gradient integrated over the boundary shared between the voxels v i and v j .To enable a comparision of the DLCM model with a corresponding PDE model the effect of surface tension needs to be included.As a consequence of this, we also need to include migration between neighboring singly occupied voxels on boundaries where surface tension is added, and then with a rate coefficient equal to D 1 .This is an event which is formally not allowed in the original framework, cf.[16].
To sum up, a population of cells occupying a grid may be evolved in time by first solving for the pressure in (2.3) using the Laplacian on the grid, and then converting the pressure gradient into rates via (2.2) and (2.4).The rates are now interpreted as competing Poissonian events for the corresponding cellular movements to be simulated as a continuous-time Markov chain.Any other dynamics taking place in continuous time are thus readily incorporated in a consistent way.

Framework tumor model
As a candidate for a 'minimal' avascular tumor model we consider the one presented in [16] which consists of a single cancerous cell type in three different states: proliferating, quiescent (i.e., dormant), and necrotic cells.As a matter of convenient implementation the range of u i can be extended to include u i = −1 which represents a voxel containing a dead necrotic cell so that u i ∈ {−1, 0, 1, 2}.Also, let Ω denote the tumor domain with u i ̸ = 0 and let Ω ext denote the entire computational domain.
An avascular tumor has to rely on oxygen and nutrients to diffuse through the surrounding tissue to reach the tumor, a process assumed to be much faster than cell migration, growth, and death.As such, nutrients are readily modeled by a stationary heat equation with a boundary condition on the external boundary ∂Ω ext (far away from the tumor boundary ∂Ω) as where c is the concentration variable understood as a proxy variable for oxygen and any other nutrients required for the cellular metabolism.Further, λ is the ratio of the oxygen consumption rate to the oxygen diffusion rate, and a(u i ) is the number of living cells in the voxel i, i.e., a(u i ) = max(u i , 0).The rates describing the tumor growth are then defined as follows: cells are in the proliferating state if c i ≥ κ prol and then divide at rate µ prol , where κ prol is the minimum oxygen concentration required for cell proliferation.A cell dies and then becomes necrotic at rate µ death if c i < κ death , where κ death is the minimum oxygen concentration required for individual cell survival.Finally, necrotic cells degrade at rate µ deg to free up the voxel they are in.Cells in voxels at intermediate oxygen levels are in the quiescent state.Note that cells instantly switch between all living states provided the oxygen concentration allows for it.At the tumor boundary a pressure condition needs to be imposed in order to capture the net effect of cell-cell adhesion as well as the interactions between cancerous and healthy cells.We let the phenomenological constant σ represent this via a Young-Laplace pressure drop proportional to the boundary curvature C. Denoting by p (ext)  An alternative approach to (2.6) for including surface tension effects directly on the microscopic level can be found in [16], where the local adhesive forces work passively to resist cell migration directly by negative contributions to (2.4).Local implementations of adhesion that do not change the population pressure field as with (2.6) are unfortunately not fully consistent with a pressure-driven migration law, thus obscuring a mechanistic understanding.The implementation of (2.6) along with other modifications to the DLCM framework are further discussed in §B.A summary of the parameters of the proposed model is found in Tab.2.1.

Mean-field PDE tumor model
We next derive the corresponding 'minimal' partial differential equations (PDE) model, constructed to very closely mimic the mean-field of the stochastic model.Certain aspects of the DLCM model, e.g., the exclusion principle in (2.4), which excludes migration events between neighboring singly occupied voxels, but also (2.5), where a nonlinear interaction between nutrients and the cell population takes place, invalidate the assumption of independence between the states of different voxels, cf.[12].This complicates formulating a mean-field PDE exactly.However, the underlying continuous physics of the DLCM model provides an appropriate starting point for the construction of a corresponding PDE model.The derivation is essentially based on mass balance with a cellular growth-and death rates, and a velocity field proportional to the pressure gradient.
Seeing as cells comprise mostly water we assume that they are incompressible.This implies that the material derivative is zero, and hence the conservation law governing the tumor cell density ρ in a velocity field v becomes where Γ is cell growth and loss due to proliferation or death, respectively, and remains to be defined.Assuming that the cells move as a viscous fluid with low Reynolds number through r p r q r n Proliferative Quiescent Necrotic Figure 2.2: Schematic of circular tumor with three distinct regions.Arrows show the distinct radii, r p , r q , and r n which correspond to the distance from the origin to respective regional interface.
a porous medium, the velocity field for the cell density is governed by Darcy's law [39] where the porous media that the cells reside in is the extra-cellular matrix (ECM), whose permeability is one of two factors determining D. Combining (2.7) and (2.8) and assuming a homogeneous permeability within the tissue domain we arrive at (2.9) For the source term Γ we mimic the stochastic model and let the cellular growth and death rates be constant and determined by oxygen thresholds.We thus define Γ as which defines the proliferative, quiescent, and necrotic region, respectively (see Fig. 2.2).These relations close the pressure relation (2.9).
Similar to the cell-based formulation, we assume a Young-Laplace pressure drop at the outer tumor boundary, which obeys (2.6).This surface tension effect arises from various cell-cell adhesion effects -the loss of which, due to loss of E-cadherin function in cells, is associated with tumor metastasis and tissue invasion [24].Let Ω here denote the region with ρ > 0 akin to the tumor domain of the stochastic model.By Darcy's law, we then have that ∂Ω − moves with the velocity where n is the interface normal vector and ∂Ω − (t) denotes the boundary as approached from inside Ω.The condition (2.11) connects the velocity field with the movement of the tumor boundary.
The pressure boundary conditions concern the external medium that the tumor grows within.Recall that Ω ext is the domain containing both Ω and the external medium.By assuming that the external medium obeys laws similar to the tumor tissue, we can summarize the medium's impact on the tumor growth through the tumor boundary conditions.We start from the assumption that the pressure propagates freely throughout the external medium outside the tumor and, for consistency of the complete two-tissue system, the external medium is assumed to abide by the same assumptions as the tumor tissue (i.e., incompressibility, Darcy's law, and conservation of mass).We further assume that growth and death of the external tissue is negligible and, thus, arrive at the governing equation for the outer pressure, p (ext) in Ω ext \ Ω: where the solution is undetermined up to a specified boundary condition.The velocities of the tumor and the external tissue are both determined by Darcy's law, but are allowed different Darcy coefficients, D and D ext , respectively.The velocities are assumed compatible at the interface and the complete set of boundary conditions on ∂Ω thus reads as ) Thus, we assume that the physical extent of the region between the tissues (where cell mixing might occur) is negligible in comparison with the spatial scale of the model.The nondimensional coefficient is expressed as Dext = D ext /D, but for brevity we omit the hat in the following analysis.The compatibility condition (2.14) allows for an investigation into the effects of varying the stiffness between the tumor and the external tissue due to, e.g., an increase in collagen density [37].For simplicity, we assume that the ECM permeability is homogeneous and time-independent across the domain of tissues, hence implicitly assuming that breakdown and remodeling of ECM that could affect tumor progression [30] are negligible.Finally, oxygen diffuses in towards the tumor through the surrounding tissue from a source (e.g., a vessel) far away with regards to the spatial scale of the system.Akin to the cell-based model, we consider a stationary heat equation but with a different source term as with c = c out on ∂Ω ext , and where λ is ratio of the consumption rate per cell density to oxygen diffusion rate in the PDE setting.We assume that ∂Ω ext is radially symmetric and lies at a distance R from the domain origin and we also make the simplifying assumption that the external tissue consumes negligible oxygen.
A suitable choice for the characteristic length is l c = R as it remains constant during growth.Nondimensionalization of the model assuming a radially symmetric tumor then yields the characteristic units with the dimensionless parameters σ = D/(µ prol R 3 ) × σ and μdeath = µ death /µ prol .We nondimensionalize the oxygen parameters independently of the pressure and arrive at the characteristic units and additional dimensionless parameter, respectively, as (2.17) Subsequently, we use the characteristic units and nondimensional parameters, but we drop the hats.The units are set to one such that R = 1, µ prol = 1, c out = 1, and D = 1.
While the PDE model is intended to closely match the mean-field of the DLCM model, there exist certain differences between the two and we view the PDE model as an effective model of the DLCM model.For the simulation of the PDE model, we thus use effective parameters that we derive from the outcome of the DLCM simulations (details in §4).

Analysis
We analyze the morphological properties of the tumor model in two spatial dimensions.This simplifies the analysis while still allowing for a qualitative comparison with in vitro data from tumor spheroids embedded in matrigel, where nutrients enter through the spherical surface (corresponding to the circular boundary in 2D).The case of a radially symmetric growth is discussed in §3.1 and a spatial linear stability analysis in §3.2.In essence, the outcome of the latter analysis include conditions for when the former radially symmetric case is a valid ansatz.Finally, in §3.3 we uncover how the morphological instabilities of the model develop during its different growth phases and a possible role of the external medium in exacerbating or reducing such effects.

Radial symmetry and the stationary state
A characteristic property of the avascular tumor is that it reaches a stationary growth phase due to limited oxygen/nutrition availability [18].We therefore first derive analytical relations that provide insight into which regions of the model parameter space map to such a stationary state under the preliminary assumption that the tumor is radially symmetric.
The constant cellular growth and death rates define distinct characteristic regions of the model tumor according to (2.10): the proliferative, quiescent, and necrotic region.For a radially symmetric tumor at time t ≥ 0, motivated by the form of the oxygen field governed by (2.15), we let r p (t) denote the tumor radius, r q (t) the radius of the interface between the proliferative and quiescent region, and r n (t) the radius of the interface between the quiescent and necrotic region, cf.Fig. 2.2.Given λ > 0, the assumed radial symmetry implies that 0 ≤ r n ≤ r q ≤ r p , and under the chosen units, r p < R = 1.We note that if r p is sufficiently close to the oxygen source at R, the model assumption of avascularity breaks down.
We simplify the problem by assuming that ρ = 1 across the entire tumor domain since this allows the oxygen field to be explicitly solved.By incompressibility and slow migration of cells, this is a reasonable approximation and, besides, a constant cell density in the PDE model is a close match to the discrete stochastic model that we wish to investigate.
We thus solve (2.15) under radial symmetry while imposing C 1 -continuity for the oxygen c across the interfaces between the characteristic regions.Under radial symmetry, the divergence theorem applied to (2.15) implies an inhomogeneous Neumann boundary condition across each radial boundary whose value is proportional to the volume of oxygen sinks contained within.The full problem for the oxygen under radial symmetry thus reads where the last three relations follow from application of the divergence theorem under radial symmetry and at each interface separately.We solve (3.1) and find where the expression does not differ in the proliferative and quiescent regions since the oxygen consumption rates are identical there.By definition, the oxygen level at r q is κ prol and at r n it is κ death .Using this, (3.2) implies the following algebraic relations between the characteristic regions in terms of the reduced parameter set {K prol , K death }.
Given radial symmetry and ρ = 1, the tumor volume change is derived from mass balance as V = µ prol V p − µ death V n , where V p and V n are the volumes of the proliferative and necrotic region, respectively.Thus, in two spatial dimensions we get that under the nondimensionalization where µ prol = 1.Hence an initial state r p (t = 0) together with the reduced parameter set {K prol , K death , µ death } fully determine the dynamics of a radially symmetric tumor under (3.3)- (3.4).Assume now that (r eq n , r eq q , r eq p ) is a stationary solution of (3.3)-(3.4).Writing r 2 q = r 2 q (r 2 p ) and r 2 n = r 2 n (r 2 p ), and by implicitly differentiating (3.3) we can linearize (3.4) around this solution and retrieve the single eigenvalue Proposition 3.1 (Stability of radially symmetric equilibrium).Given µ death > 0, assume that 0 ≤ r eq n ≤ r eq q ≤ r eq p is a stationary solution of (3.3)- (3.4).Then Λ r < 0 in (3.5) whenever r eq p ≤ exp(−1) ≈ 0.368.
Proof.Put (r eq n ) 2 = η(r eq q ) 2 for some η ∈ (0, 1) and note that (r eq p ) 2 = (1 + µ death η)(r eq q ) 2 by stationarity (the cases η ∈ {0, 1} are treated as limits).The eigenvalue becomes By inspection we find that as a function of µ death , the expression on the right is monotonically decreasing and hence it is bounded by its behavior as µ death → 0+: In turn, as a function of r eq p , this expression is monotonically increasing such that, in particular, for r eq p ≤ exp(−1) we have that From the elementary inequality −y ≤ (exp(−y) − 1) • (1 + y/2) for y ≥ 0 we conclude, taking y = − log η, that Λ r < 0. The same inequality applies also to the limits η → {0 + , 1 − }.
Proposition 3.1 provides a sufficient condition for stability based only on the size of the tumor and can be thought of as a modeling cut-off: either a radially symmetric tumor is small enough to be stable, or it has grown too large relative to the oxygen source for stability to be guaranteed for all η.While there might exist values of η for which a larger tumor is stable, the size alone for r eq p > exp(−1) is not sufficient to deduce that the tumor is stable, but the converse is true for smaller tumors.In the numerical experiments in §4, we use Proposition 3.1 to ensure that the tumor is small enough that a radially symmetric solution is expected to reach a stable stationary state.For suggested stationary radii (r eq n , r eq q , r eq p ), with r eq p small enough, µ death is defined by setting (3.4) to 0, and similarly K prol and K death are found from (3.3), which are used to determine the stability of the stationary state and the full dynamics of the tumor's growth rate.

Morphological stability
We analyze the stability of the PDE model by studying the system's response to perturbations of a radially symmetric solution.The main result depends on the three Lemmas in §A and reads as follows: Theorem 3.2 (Linear stability).Let the outer tumor boundary r p be perturbed by for some |ϵ| ≪ 1. Write the induced inner perturbations on the same form, each defined as the interface between the regions of different cellular growth rates according to (2.10) with the oxygen field governed by (2.15) (see also Fig. 2.2).We let the pressure field and the cell density advection be defined as in §2.3.Then to first order in ϵ, the kth perturbation mode in (3.6) grows as α (p) k (t) ∝ e Λ(k)t according to the dispersion relation , in which the interface perturbation coefficients α are given by (A.11) and where the radial growth follows from (3.4), by Darcy's law.
Proof.We first solve for the pressure field (2.9) under the assumption of radial symmetry where, again, the divergence theorem implies Neumann interface conditions.The result is where the negative pressure gradient at r p (3.11) recovers the velocity of the tumor's outer rim (3.4) as expected.
We continue by perturbing the tumor boundary according to (3.6) and, assuming independence between modes for the induced perturbations, we have (3.7) and let the pressure perturbation in each region i ∈ {ext, p, q, n} be Using similar arguments to those in the proof of Lemma A.3, we can show that the pressure perturbation coefficients are of the form (A.13).We next use Lemmas A.1 and A.2 to find the continuity relations for the coefficients.However, the pressure discontinuity at the tumor boundary, (2.13) and (2.14), must be treated separately.For the former, the first order approximation in ϵ of σC(r * p ) is evaluated, and for the latter we use that D and D ext are constant within their respective domains.Thus, the continuity relations become (3.13) and for the derivatives, As in [21], we find the form of this dispersion relation from the velocity at the tumor boundary by applying (2.11) to the perturbed solution.Considering only the first order terms in ϵ, we find that ∂α Finally, evaluating (3.15) using Lemma A.3 for the coefficients α k yields the dispersion relation (3.8).
Note that (3.8) is independent of the coefficients α The first part in (3.8) is the Saffman-Taylor instability [33] term.This type of instability has previously been discussed in the context of growing cell populations [28].Further, the induced perturbations on the oxygen field act to dampen any morphological instability as seen by the inner region perturbation term, which is negative and can only decrease the value of Λ(k).However, for large k this dampening vanishes in general as is seen from the following reasoning.Let r 2 n = θ n r 2 p , r 2 q = θ q r 2 p , with 0 ≤ θ n < θ q < 1, i.e., the regions do not overlap.Then the inner region perturbation term becomes which for r p < 1 tends to zero as k grows.Finally, surface tension also reduces the amplitude and range of unstable perturbation modes.Similar effects are observed due to cell adhesion in glioblastoma models in silico and in vitro [29].

Notable special cases
The dispersion relation (3.8) provides rich insight into the morphological dynamics of our model of a growing avascular tumor.We outline notable regimes of these dynamics below.
The Saffman-Taylor instability When the tumor grows in a medium that flows on a significantly smaller timescale than the migration rate of the tumor cells, corresponding to D ext ≫ D (≡ 1 by nondimensionalization), the dispersion relation becomes The same result is obtained by assuming a homogeneous outer pressure, p (ext) (r) = p 0 , for some constant p 0 .Following the same arguments we get again (3.13) but with γ (ext) k and p (ext)′ equal to zero, and that (3.14) holds except for the first relation, since C 1 -continuity is no longer valid on ∂Ω making Lemma A.2 inapplicable there.The dispersion relation is still given by (3.15) and Lemma A.3 also holds (the oxygen field does not explicitly depend on the outer pressure), which combined yields (3.17).The Saffman-Taylor term (the first term) is here at its most stabilizing since k(1 − D ext )/(1 + D ext ) ≥ −k for D ext ≥ 0. On the other side of the spectrum, we have the situation when the external tissue is significantly more viscous and practically immovable within the temporal scale of the growing tumor.This corresponds to the condition D ext ≪ 1, and (3.8) becomes and every mode k ≥ 2 is unstable during growth with no stabilizing effect from the surface tension.The case D ext < 1 is the common form of the Saffman-Taylor instability.
Growth phases The tumor's morphological stability depends on which growth phase the tumor is in.We identify the following quantity as a discriminant of the tumor's growth phase: i.e., the relative tumor boundary velocity.When the tumor is initially growing in a nutrientrich environment and r n , r q are small, then from (3.9) we see that r ′ p ≈ r p /2 which implies that ∆ θ ≈ 1 2 and hence (3.8) becomes Hence as the tumor grows exponentially, we expect all modes k > 1 to be stable for D ext > 1 and unstable for D ext < 1 in the absence of surface tension effects.
As the tumor grows larger and the nutrient availability can no longer sustain the entire tumor, the tumor front slows down and ∆ θ → 0 + , and the effect from the Saffman-Taylor part diminishes.From (3.16) we see that the inner region perturbation term in (3.8) is negative, and hence close to the tumor's stationary state we have that Since the inner region perturbation term tends to zero for increasing k, the upper bound becomes a good approximation of Λ(k) for large k.Clearly, (3.21) shows that a positive value of σ is necessary for the stationary stability.

Surface tension and stationarity
The stability relation provides an estimate of the surface tension required to maintain a radially symmetric growth as t → ∞.Considering only the case when D ext ≫ 1, we see from (3.17) that the least stable case with respect to the discriminant is obtained when ∆ θ = 0. We find the necessary surface tension by requiring Λ(k) = 0 for all k.Again, using that the inner region perturbation term is negative, we obtain from (3.17) the bound where σ stable,k is the lower bound of σ required for stabilizing small perturbations of mode k.Thus, to stabilize all modes k ≥ 2 it is sufficient to have that σ = r 3 p /6, depending only on the total tumor volume.Again, since the inner perturbation term tends to zero as k grows, the upper bound is a good approximation of σ stable,k for large k.
Creeping instability We finally remark on the interesting mode k = 1, the only mode unaffected by surface tension.Geometrically, this mode corresponds to movement of the tumor's center of mass: the tumor begins to creep towards the oxygen source given a small perturbation.From (3.8), Thus, the tendency for creeping always exists when r q and/or r n are positive, and the model would require additional features regarding, e.g., the external tissue's response to invasion, in order to inhibit this effect.Note that this tendency is reduced for tumors growing within an environment more viscous and/or less permeable than its own.As showed previously, however, such conditions make all the other modes k ≥ 2 less stable.

Numerical examples
In the following section, we present numerical simulations of both the stochastic model from §2.2 and the PDE model from §2.3.We focus on the case D ext ≫ 1 which we showed had an inherent stabilizing effect during growth in §3.3.In §4.1 we assess the validity of the assumption that growth is radially symmetric and how the stability responds to surface tension.In §4.2, we explore the relation between surface tension and the emergent morphology and compare the outcomes between the stochastic and the mean-field PDE model.Due to certain differences between the models, we use effective parameter values for the PDE simulations for the parameters µ prol , µ death and λ, and denote those with a bar, e.g., μprol .The effective parameters are derived from the stochastic model simulations via basic scaling considerations or preliminary simulations as detailed in §B.
The set of parameters used in both the DLCM and the PDE simulations are found in Tab.4.1.

Radial symmetry and surface tension
We first solve the PDE under the assumption of radial symmetry.We solve the reduced 1D problem comprising (3.3) and (3.4) as derived in §3.1 and evaluate the perturbation growth rates (3.17) during tumor growth and study their response to surface tension.For quantitative measurements of the regional characteristics during tumor growth, we consider the volumetric quantities V p = πr 2 p , V q = πr 2 q , and V n = πr 2 n .Fig. 4.1 shows the evolution of the regional characteristics for a simulation using the standard parameters for the PDE found in Tab.4.1, accompanied by the perturbation growth rates (3.17) at the stationary state.In Fig. 4.1a we observe the emergence of the characteristic sigmoidal growth of the total volume, with an initially exponential growth followed by a growth rate that plateaus.Fig. 4.1b shows the perturbation growth rates versus mode close to the stationary state for a range of σ.It is clear that the assumption of radial symmetry for low values of σ does not hold when oxygen is not sufficient to sustain the growth of the entire population.From (3.17) we find that σ stable,2 ≈ 3.1 • 10 −3 , i.e., this is the value needed to   stabilize modes k ≥ 2. From Fig. 4.1b we see that σ lower than this prompt nontrivial spatial behavior with instabilities that may occur over different timescales (investigated further in §4.2).As suggested by (3.23), creeping is expected for long enough times.We test this by simulating the model using σ = 3.2 • 10 −3 to ensure that modes k ≥ 2 are stabilized (see details about the numerical methods and the implementation of surface tension in §B).The results are shown in Fig. 4.2 where we see that the tumor reaches close to its stationary state at around t = 10, in accordance with the solution to the 1D equations in Fig. 4.1.We observe a notable collective migration from the domain origin starting at t = 80.

The emergent morphology and its response to stochasticity
We begin by a brief investigation into how well the dispersion relation (3.17) describes the morphological stability of the DLCM model tumor.To this end, we conduct simulations of the complete DLCM model, initialized close to the estimated equilibrium volumes (see Fig. 4.1a) using the standard parameters and setting σ = 10 −4 .We impose an initial perturbation to the tumor geometry of the form (3.6) with ϵ = 0.05 and for modes k = 1, ..., 8.We regard the mean of (3.17) during a selected time interval as an analytical prediction and measure the growth of each mode by fitting the amplitude to an exponential growth law.Fig. 4.3 indicates that the rates agree fairly well, thus supporting the use of the PDE-based stability analysis in predicting the behavior also of the DLCM model.Motivated by the quantitative evidence for a correspondence between the PDE analysis and the DLCM model, we compare the morphology and the growth patterns of the stochastic model and the PDE model for similar parameters.To avoid perturbations that are biased by the discretization method we add small amounts of white noise to the cell density updates.These and further details of the PDE solver, including the implementation of surface tension, are discussed in §B.We evaluate the tumor boundary roundness defined over a 2D region as which ranges from 0 to 1, where 1 means the shape is perfectly circular and smaller values measure its deviation from circularity.We first investigate numerically the effects of surface tension on the tumor morphology.Specifically, we study the onset of second mode instability according to (3.17) for a tumor growing in two spatial dimensions.For this purpose, we compare the growth using σ = 2•10 −3 and σ = 5 • 10 −4 until t = 30, during which the former value is stable for k = 2 although it is not stable over larger times t.For both experiments, we compare morphology and growth using a single simulation per model (a discussion on the impact of stochasticity is offered in §5).Fig. 4.4 shows the evolution of the tumor's characteristic volumes for both models together with the tumor roundness (4.1) using the larger value of σ.Fig. 4.5 shows snapshots of the solution corresponding to Fig. 4.4.Both solutions remain close to being radially symmetric during the full simulations since the small second mode instability does not show during these relatively short time intervals.Notably, the creeping effect becomes apparent earlier for the DLCM simulations.
Similarly, Fig. 4.6 and Fig. 4.7 shows regional evolution and spatial solutions, respectively, for the lower value of σ.Both models display a significant decrease in roundness (accompanied by a total volume increase) as the tumor begins to split in two some time after the growth has plateaued.A notable difference between the growth curves during this process is that the DLCM tumor does not reach a fully stationary state before the splitting, possibly due to the higher exposure to noise in the stochastic model.Finally, Fig. 4.8 shows the resulting morphology of both models using four different and smaller values of σ.For these experiments we use use a lower κ death = 0.92 for a thinner proliferating rim which results in effective parameters μdeath = 1.0 and λ = 1.1.Again, the radially symmetric tumor displays significant creeping only for the DLCM model at the selected final time (cf.Fig. 4.8a and Fig. 4.8e).The morphologies of the tumors are similar in terms of emergent unstable modes and sizes of the characteristic regions, e.g., the tumors beginning to separate into two is seen in both Fig. 4.8b and Fig. 4.8f.For the most unstable case using σ = 0 in Fig. 4.8d, the DLCM tumor grows somewhat larger and with different morphological and regional characteristics.Finally, Fig. 4.8a and Fig. 4.8c display small cell clusters detaching from the tumor to grow on their own, a phenomena that we never observe in our simulations of the PDE model (cf.Fig. 4.8e and Fig. 4.8g).

Discussion
We have analyzed the morphological stability of a PDE model of avascular tumors.The PDE was derived to closely represent the mean-field of a stochastic model expressed in the DLCM framework.Assuming radial symmetry in the PDE model, we first characterized the growth dynamics as well as the stationary state.A linear stability analysis in two dimensions was subsequently carried out and we found a dispersion relation that describes the stability of the morphology of the tumor.Finally, we compared the analytical predictions with numerical simulations of both the stochastic DLCM and the PDE model.The observed morphology of the stochastic model was found to be in line with the predictions from the PDE analysis, including also the proposed relations required for a stable stationary state.
The Saffman-Taylor instability acts on the tumor boundary, where the determining factors are the porous medium permeability and the tissue viscosity as summarized in the coefficients was shown analytically in §3.3 and experimentally in §4.1.This also explains the asymmetry and unlimited growth of the original DLCM tumors presented in [16], which did not implement an explicit surface tension effect.Morphological instability arising in conditions when nutrients are scarce is in line with analyses and simulations of other models of tumor growth and cell colonies.The model in [10] is a fluid-based PDE model including Darcy flow, with cell growth proportional to nutrient level, and a constant apoptosis rate.The model boundary velocity is explicitly dependent on the nutrient gradient, in contrast to the model analyzed herein that nevertheless displays a similar nutrient dependent growth instability.Further examples show similar instabilities in nutrient-deprived conditions both in agent-based models of cell colonies [20] and in hybrid models of tumor growth including cell cycles and ECM [1].The former model tumor expands solely due to cell proliferation and not to pressure-driven migration, thus suggesting the persistence of this type of instability across a range of models.
Interestingly, our model predicts a creeping effect in which even an otherwise stable tumor as a whole migrates towards the oxygen source.Experiments in §4.2 suggest that the stochastic model has a somewhat stronger tendency for creeping than the PDE model.The larger noise levels of the DLCM model is a good candidate explanation for this difference.Moreover, the creeping effect cannot be inhibited by surface tension, must be controlled through additional mechanisms of the model such as an elastic response from the external tissue (see [38] for a review on minimal morphoelastic tumor models).Thus, the creeping phenomena prompts the following questions: Does creeping occur in vitro or in vivo and if so, over what timescales?If not, what mechanisms keep it from occurring; alternatively, what assumptions could make our model more realistic in this regard?
An additional observation from the DLCM model was the detachment of cell clusters even in the case of a surface tension large enough to support a radially symmetric solution.This is most likely due to the discrete and random nature of individual cells in the model.Detachment is therefore a distinct feature of on-lattice stochastic modeling in this context, which appears to be a more realistic representation of tumors growing under noisy conditions than a purely fluid mechanical continuous model can offer.That said, the stochasticity of the DLCM model does not have a significant impact on the region volumes and final size of a radially symmetric tumor.Rather, assuming stability, these outcomes are fairly accurately governed by the deterministic relations (3.3)- (3.4).The process noise does, however, have an impact on the morphology of the tumor under less surface tension, which implicitly affects the tumor size first when the boundary has been significantly distorted.A deeper investigation of this is outside the scope of the experiments reported here.
One fundamental difference between the two models is the spatial exclusion principle which is implemented in the DLCM framework via the carrying capacity.A consequence of this can be seen when comparing Fig. 4.7c and Fig. 4.7f, where the PDE tumor is close to separating into two pieces, while the DLCM tumor in contrast retains an oval shape.The latter is due to necrotic cells which degrade while still occupying voxels, thereby slowing down the mass flow.Similarly, Fig. 4.8d and Fig. 4.8h display significant differences in morphology and size, where the former model supports a larger necrotic region.These examples highlight emergent differences between the two ways of modeling cell extent and migration and call for an input of biological observations to approach a higher level of realism.
We end by briefly mentioning some potential modifications that may improve on the expressive power of the model.Considering the PDE model first, we see from (3.11) that the ambient pressure becomes very large for a stiff external medium.Thus, when modeling such scenarios the addition of pressure-dependent effects such as pressure-driven oxygen flow or a pressure-based proliferation rate become relevant.Such additions carry over to the stochastic framework in a fairly straightforward manner.To further improve on the realism of the nutrient modeling, a limit on the diffusion flux of the oxygen into the tumor across its boundary could also be considered (cf.[18]).While these model modifications are readily implemented, the resulting emergent behavior is not obvious and a precise mathematical analysis is more involved due to the increase of nonlinear feedback mechanisms.
In conclusion, our basic stochastic avascular tumor model turned out to be a fruitful tool in leading to some new insights as well as investigations of the effects of various couplingand feedback relations.Building and analyzing the mean-field PDE in tandem with the discrete stochastic model forced us to think more deeply in terms of trade-offs for continuum models of discrete cellular agents; this approach limits model refinements to a certain extent since the discrete and the continuous versions need to be consistent.We anticipate that the combination of a mean-field PDE and a stochastic model built from first principles will enable the development of filtering tools aimed specifically at integrative Bayesian approaches to data-driven applications.For example, the stochastic model can be used to characterize the process noise of a Kalman filter that leverages the PDE model as its state transition.

Availability and reproducibility
The computational results can be reproduced with release 1.4 of the URDME open-source simulation framework [3], available for download at www.urdme.org(see the avascular tumor examples and the associated README in the DLCM workflow).

PDE solution
The pressure and oxygen fields (2.9) and (2.15), respectively, are readily solved by a finite element method (FEM).We solve for each quantity using linear basis functions and, for simplicity, using a lumped mass matrix.Solving for the oxygen, however, requires special considerations as the oxygen consumption depends on the oxygen level at each time step.To address this, we solve for the oxygen iteratively within a smaller time scale, denoted τ , while holding t n and all other variables that depend on t n constant.We solve the resulting time-dependent heat equation implicitly using the Euler backward method, while treating the source-term explicitly.We iterate in pseudo-time τ until |c τ +1 − c τ | is less than some tolerance after which we set c(t n+1 ) = c τ +1 .The inner iterations are given by where Ω h (t n ) is the discretized tumor domain at time t n .The advection of cell density according to (2.7) is solved by an explicit finite volume (FV) scheme using the pressure distribution at a given time step.We consider a square Cartesian grid and an initial uniform cell density of arbitrary shape placed within an otherwise empty domain where ρ = 0. We keep track of the moving boundary by finding volumes where ρ < ρ thresh .We apply the pressure boundary condition (2.6) on such volumes that are simultaneously adjacent to volumes with ρ ≥ ρ thresh (details on the implementation of surface tension are found in the next paragraph).The threshold was determined to ensure a front speed consistent with the compatibility condition (2.8) during tumor growth and the value ρ thresh = 0.9 was chosen to reduce smearing effects close to the boundary and maintain ρ ≈ 1 across the tumor domain.Furthermore, cell densities outside the tumor boundary should not consume oxygen, and we therefore set λ = 0 for volumes with ρ < ρ thresh .Using a first order upwind scheme to solve (2.7) including a noise term amounts to the following FV scheme , where V i is the volume of the volume element, ∂V i is the boundary of the volume, e ij is the unique edge adjacent to volumes V i and V j .The quantities ρ n i and Γ n i are the volume average quantities of ρ and Γ, respectively, at time t n taken over volume V i .Finally, to excite all perturbation modes of the boundary without discretization bias, we have added scaled Gaussian noise to the cell density change as the last term in (B.2) dictates, with N ∼ N (0, 1) and ω = 0.025.The form of the noise is suggested by a Wiener process approximation of a Poissonian interpretation of the drift term in (B.2) at a system size ω −1/2 .For our simulations we used a Cartesian discretization of the mesh with edge length ∆x = 0.02 and ∆t = min(∆x, 0.1/ max i (|F n i |)), which yielded stable solutions across our experiments.

Surface tension
The curvature C in the pressure boundary condition (2.6) is evaluated numerically on the tumor boundary by estimating the derivatives on the corresponding con-tour lines.At each time step, we find the contour lines of the tumor boundary where ρ = ρ thresh for the PDE model, and where u i = 1 for the DLCM model.These contour lines are then parameterized by (x c (s), y c (s)) for s ∈ [0, S] assuming periodicity, i.e., (x c (s), y c (s)) = (x c (s + S), y c (s + S)), since they are closed curves.We find the cubic spline of each coordinate such that after differentiation we obtain the signed curvature We then calculate the Young-Laplace pressure difference at each boundary volume by using the curvature given by (B.3) at the contour point closest to the volume.We include a simple threshold in our implementation to avoid that a pressure boundary condition is applied to undesired contours.Specifically, we neglect contour lines that consist of a number of points less than a fraction f min = 0.95 of the largest contour.Thus, holes within the tumor or clusters of cells outside the tumor do not produce any surface tension.This prevents spurious surface tension and pressure values, which would otherwise be present in particular in the stochastic model.

DLCM modifications
The physics of the PDE model motivates two features of the stochastic mode not present in the original presentation [16], namely surface tension and pressure sinks due to mass loss.Surface tension is implemented in the same manner in DLCM as for the PDE solver.However, the exclusion of certain events in DLCM due to the voxel carrying capacity introduces a bias associated with the strength of surface tension.Therefore, we allow cells in singly occupied voxels on the boundary to migrate into the population with rate coefficient equal to D 1 .Without inwards migration, there is nothing to balance out the sporadic outwards migration induced by surface tension on the inherently noisy boundary.In particular, allowing inwards migration ensures that higher levels of surface tension do not introduce a bias in the average front speed of the tumor due to this noise.We note that this bias is unavoidably present in the initial growth phase of the tumor when no attractive, necrotic region exists, as can be observed in the higher initial rate of growth of Fig. 4.4a compared with Fig. 4.6a.
The PDE model expresses that mass loss due to cell death produces pressure sinks that drive cell migration to replace the lost mass, and that the pressure sinks are proportional to the rate of mass loss, see (2.10).In the stochastic setting, cell mass loss occurs for necrotic cells that are degrading at the rate µ deg , and for simplicity we use that the pressure sinks and pressure sources are equal to µ death and µ prol , respectively, akin to the PDE.Further, to ensure that the degrading cells are pushed by the population pressure to the center of the necrotic region before they have fully degraded, we derive an approximate value of the degradation rate µ deg depending on µ death for this to happen.Consider a fixed, radially symmetric necrotic region and the motion of a particle from the region's boundary into its center governed by (2.8).The pressure gradient at any point in this region is found from (3.12) and letting the particle start at r(t = 0) = r n , its position is given by r(t) = r n e −µ death t/2 . (B.4) The expected time τ it takes for the particle to reach r(τ ) = δr n , 0 ≤ δ < 1 of the necrotic region's center is thus readily obtained from (B.4).The condition on µ death that ensures that the expected time for the particle to degrade is equal to τ is then given by 1/τ as µ deg = µ death 2 log(1/δ) .
For our implementation, we assume that δ = 0.01 is enough to ensure that migration of necrotic cells to the tumor center is possible and occurs frequently.

PDE effective parameters
The effective parameters of the PDE model corresponding to the DLCM model parameters are derived by comparing the models radially symmetric solutions.There are three parameters that differ in value due to model differences: µ prol , µ death , and λ.Firstly, there is a difference in volumetric growth rates between the models since the DLCM tumor must undergo both a proliferation event and a migration event separately to expand in size.From a DLCM model simulation using the parameters in Tab.4.1 and σ = 0, we consider the initial exponential growth phase (r n = r q = 0) such that the solution to (3.9) reduces to r p (t) = r p (0) exp(µ prol /2 t).We fit observations of r p during this phase by minimizing the mean-square error to get an estimate of the corresponding effective growth rate of the PDE model, and we find that μprol ≈ µ prol /2.7.Since the dimensionless parameter μdeath is proportional to μ−1 prol , the former is scaled accordingly.Secondly, we note that since doubly occupied voxels consume twice the amount of oxygen in DLCM, the effective λ can be expected to be greater for the corresponding PDE assuming identical characteristic regions.The effective parameter λ is estimated by taking the mean of (3.3) using the fixed DLCM model parameter κ prol and observed values of (r n , r q , r p ) during a stationary phase of the DLCM simulation.This gives the estimate λ ≈ 1.15.The same procedure yields μprol = 1.0 and λ = 1.1 for the other set of experiments that uses different parameters shown in Fig. 4.8.

Figure 2 . 1 :
Figure 2.1: Cell population representation in the DLCM framework using a Voronoi tessellation.Green voxels represent singly occupied voxels and red voxels represent doubly occupied voxels as shown explicitly in Fig. 2.1a, where the ellipses indicate the corresponding underlying population of cells.The grey area highlights a region of empty voxels.The doubly occupied voxels exceed the carrying capacity and exert pressure on the surrounding cells.The movement rates are visualized in Fig. 2.1b, where the arrows represent movement rate and direction for a subset of the possible movements.The grey voxels represent empty voxels that cells may migrate into.Adapted from [16, Fig. 2.1].

k
are proportional to α (p) k , as seen in(A.11).It follows that Λ(k) is unambiguously determined by the set of parameters and values {r p , r q , r n , µ death , D ext , σ}.

Figure 4 . 1 :
Figure 4.1: Regional characteristics and perturbation growth rates under the assumption of radial symmetry.Fig. 4.1a shows the evolution of the three characteristic regions.Fig. 4.1bshows the perturbation growth rates (3.17) for the first few modes at the tumor's stationary state and at varying surface tension σ.

Figure 4 . 3 :
Figure 4.3: Comparison between analytical (3.17) and experimental perturbation growth rates for the first few modes when σ = 10 −4 .The error bars show the standard deviation of the analytical values during the time interval of measurement, and the shaded region corresponds to the standard deviation of the growth rate estimate during the same interval (20 independent runs).

Figure 4 . 4 :
Figure 4.4: Evolution of characteristic volumes for DLCM (a) and PDE (b), respectively, using the standard parameters in Tab.4.1 and σ = 2 • 10 −3 .The 2D solutions (dashed) are compared with the 1D solution (3.3)-(3.4)(solid) for the same parameters.Blue shows the roundness (4.1) of the tumor boundary over time with moving window standard deviation in shaded.The vertical lines indicate the times t ∈ [0, 25, 30] of the solutions shown in Fig. 4.5.