A Mathematical Model for the Onset of Avascular Tumor Growth in Response to the Loss of P53 Function

Abstract

We present a mathematical model for the formation of an avascular tumor based on the loss by gene mutation of the tumor suppressor function of p53. The wild type p53 protein regulates apoptosis, cell expression of growth factor and matrix metalloproteinase, which are regulatory functions that many mutant p53 proteins do not possess. The focus is on a description of cell movement as the transport of cell population density rather than as the movement of individual cells. In contrast to earlier works on solid tumor growth, a model is proposed for the initiation of tumor growth. The central idea, taken from the mathematical theory of dynamical systems, is to view the loss of p53 function in a few cells as a small instability in a rest state for an appropriate system of differential equations describing cell movement. This instability is shown (numerically) to lead to a second, spatially inhomogeneous, solution that can be thought of as a solid tumor whose growth is nutrient diffusion limited. In this formulation, one is led to a system of nine partial differential equations. We show computationally that there can be tumor states that coexist with benign states and that are highly unstable in the sense that a slight increase in tumor size results in the tumor occupying the sample region while a slight decrease in tumor size results in its ultimate disappearance.

Introduction

There is solid evidence for the existence of a tumor suppressor gene known as p53. The p53 protein functions at several levels to control cell growth. It can inhibit the transcriptional activity of other proteins such as Sp1 and HIF-α, both of which activate growth factor (GF) gene expression [51, 50, 57, 60, 65]. It can also inhibit cell mitosis and regulate cell apoptosis [50]. In particular, one can view p53 as a tumor growth inhibitor and Sp1 as a tumor growth stimulator. (There are other cell-generated inhibitors and stimulaters but we shall confine our attention to these two.)

The purpose of this paper is to present a biochemically based model for the onset of solid tumor formation. The model describes the onset of solid tumor formation with diffusion limited growth on a scale consistent with diffusion limited avascular tumor sizes. A schematic summary of the underlying biochemical pathway is given in Figure 1. This paper does not go beyond this point to treat the onset of vascularization due to the diffusion of the carcinogenic waste products of tumor necrosis (tumor angiogenesis). This idea is not new in the sense that in [14] the authors used the kinetics of cell metabolism, to model the growth of tumor cells in the micro environment.

Figure 1

Biochemical schematics: This figure may be regarded as the biochemical underpinning of our model. (Such figures are sometimes informally called “wiring diagrams” or “network diagrams”.) For example, benign cells express the wild type p53, which in its turn inhibits the transcription factor SP1 from initiating TGFα expression. Hence the blunt arrow leading from p53 to Sp1. No such inhibition is possible from the mutant p53 and hence is not shown in the figure or included in the model. Similarly, the double arrows express the fact that both cell types will secrete TGFα and in turn are influenced by this growth factor. The pathway from TGFα to the ECM proteins is really a schematic for the cell production of the matrix metallo-protein, which in its turn degrades the ECM via “standard” enzyme kinetics. Such diagrams, popular in the biochemical and molecular biology community, contain the seeds of systems of differential equations that begin with a consideration of the Law of Mass Action. In this sense, our model is in the spirit of [14] where the authors modeled the citric acid cycle in order to predict the onset of solid tumor growth in the micro environment.

It is not our intention to assert that this is the only way avascular tumor growth can arise from gene mutation.

Among many examples are cancers such as lymphoma that are driven by an overexpressed and mutated MYC gene [7, 32]. The wild type MYC protein both inactivates the cell cycle by suppressing p21 and suppresses apoptosis by inducing BIM. The mutant forms of the MYC protein found in these cancer fail to induce BIM [see ref below]. This pathway runs in parallel with the p53 pathway, which is also activated by MYC ([32]). The mathematical modeling of this biochemical pathway is very similar to the the pathway we consider here in many respects.

The goal here is to combine the loss of p53 regulatory function in benign cells with the consequent increase in the rate of benign cell loss (via mutation into tumor cells) and the effect this mutation has on the rate of formation tumor cells in order to model solid tumor growth in the micro environment. The model attempts to duplicate nature in the sense that, in the model, it is assumed that a few benign cells undergo a mutation by which they lose the ability to express the wild type p53 gene. More precisely, it is assumed that a few benign cells undergo a mutation by which there is a loss of wild type p53 expression over a small spatial region for a fixed time period. (Following the usual biochemical convention, when we refer to p53, we will only mean the wild type because the mutant type does not regulate Sp1 (and hence growth factor) for one of several reasons including (a) an incorrect protein transcribed from the mutant RNA, (b) no mutant p53 at all or (c) a mutant p53 with no suppresser capability.) This loss of tumor suppressor function will, in the model, result in excess Sp1 production. Sp1 functions as a transcription factor for, among other genes, the gene for MMP-1, a protease that is known to degrade tissue structural proteins such as collagen. This idea of genetic instability is not new (see [69]) for example. However, we think the modeling approach, as it combines fundamental biochemistry with chemotaxis and modification of classical competition dynamics, is novel.

We begin with an initial state in which we have a region that consists only of benign cells and matrix. This state constitutes a stationary solution for a dynamical system based upon the biochemistry described above in the absence of loss of function. Our results show that the length of the period of mutation (the time interval over which the loss of function occurs), the spatial size (the volume over which loss of function occurs) and the intensity of the mutation (percentage of cells in which loss of function occurs) determine whether the tumor cells take over the region of interest. In biological terms, the period of the mutation refers to the period during which the tissue is exposed to a mutagen, the spatial size refers to the volume of tissue that is exposed to the mutagen, and the intensity of the mutation refers to the concentration and potency of the mutagen. For example, the frequency of p53 loss of function in the skin is influenced by the period of exposure to uv light, the area of skin exposed to the light and the intensity and spectrum of the light.

Mathematically, the loss of function is modeled by the simple device of lowering the effective p53 rate constant for a specified portion of time and over a specified region in the benign cell region for a specified percentage (intensity). (See equation (2.1) and notice the modification in the coefficient k_p in equation (2.2)). The mutation in cell type is then modeled by a corresponding loss of benign cell growth rate coupled with an increase in tumor cell cell growth rate whenever the concentration of p53 falls below a threshold level. (See equation (3.5) and notice the last term on the right hand side of each rate law.)

In particular, the model shows that:

When the mutation is sufficiently intense, widespread or of sufficiently long duration, the malignant cell mass that develops during the loss of function period contains sufficient malignant cell mass to take over the region under consideration.

On the other hand, if the mutation is insufficiently intense, very local or of insufficient duration, the initial transfer of some cells from the benign state into the tumor state induced by the mutation is not sustainable and the small tumor perturbation dies out, with the system returning to its initial state.

The model postulates that if the intensity is small enough, and the region over which the mutation occurs is fixed, the duration time can never be long enough to initiate a growing tumor.

The model postulates the existence of tumor states that can coexist with benign states for very long times (i.e. form coexisting quasi-steady states.)

Some comments on the mathematical literature: There is a large literature on the mathematical modeling of various aspects of solid tumor growth. See [1, 2, 3, 4, 13, 12, 15, 35, 37, 62] for some representative examples. Some of this literature reflects the controversial nature of modeling solid tumor growth. See for example [30] and the rebuttal [58]. There is also an extensive literature concerned with the spontaneous remission of solid tumor growth, i.e. tumors that disappeared without treatment of any kind. A key word search on the phrase “spontaneous remission or regression of cancer” in PubMed led to over 2200 screens of articles. Many of the articles listed were case studies of various identified cancers. A few dealt with various mathematical aspects of the problem [36, 70]. We suggest that our results demonstrate one possible mechanism for this remission.

There are a number of papers dealing with competition models and Lotka-Volterra systems, as systems of ordinary differential equations and as systems of reaction diffusion equations, see for example [11, 18, 21, 22, 19, 20, 34, 38]. However we could find no papers in which these systems were combined with Keller-Segal type equations of chemotaxis. (We did a keyword search on MathSciNet and found no matches for “competition species” AND “chemotaxis” or “competition species” AND “Keller-Segal”.) Our model does combine chemotaxis for cell movement together with the “competition species” notion for competing species. Given the dearth of literature on these combined notions, perhaps we can claim some novelty for our model. It is not only mathematically naive but also false, to draw the inference that conclusions about the dynamical behavior of solutions of a system ordinary differential equations such as the Lotka-Volterra system for competition systems, will lead to valid conclusions about analogous systems of partial differential equations when these equations are combined with diffusion terms or diffusion with chemotactic terms. Turing instability is probably the classical example of the failure of the stability properties of nonlinear ordinary differential equations to be inherited by systems of partial differential equations with the nonlinearities arising from the ode system. See also the two very nice papers [45, 63] for alternate approaches to this issue.

Here is an outline of the remainder of the paper. Sections 2–8 constitute the main body of the paper.

•

Section 2: We present the biochemical and cell biological underpinnings of the model. The details of the biochemistry are set forth in Appendix A.

•

Section 3: We discuss the cell movement equations. Here we use the notion that cells move up or down chemical gradients (chemotaxis) together with the idea that mitosis of each cell type interferes with the mitosis of the other.

•

Section 4: We set forth the mechanistic equations that define the notion of diffusion limited tumor growth. That is, the entire system is governed by the consumption of nutrients delivered to the tumor region via diffusion. Waste products are generated as a consequence of protein degradation and cell apotosis. Because we do not use the waste products further as sources for tumor growth, we ignore their evolution here.

•

Section 5: We give the appropriate boundary conditions and initial conditions for the problem dynamics.

•

Section 6: We discuss the dynamics predicted by a subsystem of ordinary differential equations.

•

Section 7: We describe how the normalizing constants are selected in order to nondimensionalize the system and how the initial steady state is found. First the system (without mutational considerations) is modified in such a way that constant solutions of the modified system can be found explicitly. These constant values are used to non-dimensionalize the full problem. By using the constant values as initial conditions for the system (without mutational considerations) the putative non-constant, benign solution is obtained (numerically) as the time independent limit of the system.

•

Section 8: We present the results of our simulations and our conclusions. We consider our problem with spherical symmetry to keep the computational complexity to a minimum.

•

Appendix A: We give a brief overview of the molecular biology and the resultant chemical kinetics upon which we base our model.

•

Appendix B: We record the numerical values of the constants we used in the simulations.

The reader is warned that quite a bit of mathematics lies ahead. We do not apologize for this, since in order to use mathematics to describe a complex biochemical process, one might expect that the mathematical description itself might likewise be involved. However, to make the reader's task somewhat more palatable, we offer the following observations. Nothing more complex than elementary enzyme kinetics, mass action and Fick's law is involved in the description of the dynamics of the time evolution of the protein and nutrient concentrations as well as the local cell population densities. The nutrient equation controls the flow of nutrients reaching every point of the extracellular matrix.

Notation and Kinetics

We consider the interplay among the proteins p53, GF the transcription factor, Sp1, the extra-cellular structural proteins F, and the matrix metalloproteinase, MMP-1 as described in Figure 1. In this model, it is assumed that transport of subcellular species such as growth factor, p53 etc. is passive, i.e. is controlled by cell movement. We imagine a region D in two or three dimensional space. Throughout, X = (x, y) or X = (x, y z) according to the dimension in which D lies. We use the notation in Table 1. We have eleven species, nine of which are chemical, two of which are cellular. (In the sequel the constants N_B, N_T denote the carrying capacity densities for benign and tumor cells respectively.) Throughout this paper we adopt the convention that the Greek μ is reserved for decay (turnover) constants while V is reserved for equilibrium constants.

Table 1. Notation.

The standard chemist's notation for the concentration of a species X is [X]. Using such a notation for functions which depend on spatial as well as temporal variables, especially in combination with sub and superscripts, is clumsy. To simplify this, we use the notation in the table below. Furthermore, if a molecular species is associated with a benign cell or with a tumor cell we use the subscripts B, T with its label. For example, S_B,S_T refer to the concentrations of Sp1 in the benign and in the tumor cells respectively.

Species	Notation	Function form	Equation number
p53 protein	P	p(x,t)	(2.1)
Sp1, benign cell derived GF transcription factor	T_r	S_B(x,t)	(2.6)
Sp1, tumor cell derived GF transcription factor	T_r	S_T(x,t)	(2.7)
Sp1 in activated state	$T_{r}^{a}$	$s_{B}^{a} (x, t)$	(2.8)
Sp1 inhibited by p53	$T_{r}^{i}$	$s_{B}^{i} (x, t)$	(2.8)
GF, growth factor	V	v(x,t)	(2.11)
MMP-1, matrix metalloproteinase-1	M	m(x,t)	(2.12)
Tissue collagen/fibronectin proteins	F	f(x,t)	(2.13)
nutrients (amino acids, lipids, oxygen, sugars, etc.)	Y	Y(x,t)	(4.1)
benign cell density	B	η _B (x,t)	(3.5)
malignant cell density	T_m	η _T (x,t)	(3.5)

We begin with p53. A relationship for the rate of formation of this transcription factor based on the number of benign cells present and level of resources is given by

\frac{\partial p}{\partial t} = \frac{K_{p} Y (x, t)}{K_{p} + Y (x, t)} \frac{η_{B} (x, t)}{N_{B}} - μ_{p} p (x, t) .

(2.1)

The assumption here is that only benign cells express the wild type gene, wherever they appear in the tissue. The cell density has been re-scaled so that k_p now includes the maximum possible benign cell density as a factor. (In Appendix A, we give a rationale for the term $\frac{K_{p} Y (x,t)}{K_{p} + Y (x, t)} \frac{η_{B (x, t)}}{N_{B}}$ based on chemical kinetics.) We refer to (2.1) as the equation of p53 evolution without loss of function.

In order to express the idea of evolution of p53 with loss of function, consider a subregion D′ ⊂ D, a time interval (more generally, several time intervals) I = (T₁ T₂) and a function Ψ(x, t) with support D′ × I and having values in (0,1]. A mutation is said to occur in the region D′ ⊂ D over a time interval I if there has been a loss of p53 expression in the set D′ × I. That is, we replace (2.1) by

\begin{matrix} \frac{\partial p}{\partial t} = \frac{K_{p} [1 - Ψ (x, t)] Y (x, t)}{K_{p} + Y (x, t)} \frac{η_{B} (x, t)}{N_{B}} \\ - μ_{p} p (x, t) . \end{matrix}

(2.2)

We refer to this equation as the evolution equation for p53 with loss of function or as the mutated p53 equation. The function Ψ is called the loss of function coefficient. The intensity of the loss of function is defined as

ℐ = \frac{1}{V_{D^{'}} L (I)} \int_{D^{'}} \int_{I} Ψ (x, t) d t d x

(2.3)

where V_D′L(I) is the product of the volume of D′ and the length of the time interval over which the loss of function occurs. Notice that 0 ℐ ≤ 1. The strength of the mutation is defined as

S = V_{D^{'}} L (I) ℐ .

(2.4)

Corresponding to such a loss of function coefficient will be a drop (pointwise) in p53 concentration. When this concentration falls below a certain level, p_c say, some benign cells will become tumor cells. The model is not reversible. That is, tumor cells can only disappear via apoptosis. For the time being, it is enough to remark that each equation (2.1), (2.2) will play a role as part of a larger system. The system consisting of (2.1) and the remaining dynamics (discussed below) will be used to compute a non constant steady state in which there are no tumor cells present. This steady state will then form the initial condition for the system consisting of (2.2) and the dynamics discussed below.

Remark 1

There is experimental evidence that p53 in normal tissue is much more stable than in malignant tissue. For example, in [6] it was reported that although this protein has a half life of about three hours in normal mammary epithelial cells, it was approximately 15 minutes in E5 immortalized cells, ie. cells that were immortalized by the E6 gene of HPV-16, the human papillomavirus commonly associated with cervical cancers. We have not included this observation as it is clear intuitively that its effect should drive up the expression of active Sp1 and hence the level of growth factor over what we already have computed in the regions where tumor cells are present.

From the literature [51] we know that p53 inhibits the transcription factor Sp1. We assume that this occurs via the equilibrium:

T_{r}^{a} + P \overset{V_{e}^{s}}{⇄} T_{r}^{i}

(2.5)

where $T_{r}^{a}, T_{r}^{i}$ represent the active and p53 inhibited form of Sp1. We write:

\begin{matrix} \frac{\partial_{s_{B}}}{\partial t} = \frac{k_{s} Y (x, t)}{K_{s} + Y (x, t)} \frac{η_{B} (x, t)}{N_{B}} \\ - μ_{s_{B}} s_{B} (x, t) \end{matrix}

(2.6)

for Sp1 growth and decay in benign cells. In tumor cells we write

\begin{matrix} \frac{\partial_{s_{T}}}{\partial t} = \frac{k_{s} Y (x, t)}{K_{s} + Y (x, t)} (\frac{η_{T} (x, t)}{N_{T}}) \\ - μ_{s_{T}} s_{T} (x, t) \end{matrix}

(2.7)

where $μ_{s_{T}} \leq μ_{s_{B}}$ since the decay rate (the turnover rate) for Sp1 in benign cells is larger than that for tumor cells in general. We have, in terms of the notation in Table 1, for benign cells only

S_{B}^{a} + S_{B}^{i} = S_{B}, S_{B}^{i} = V_{e}^{s} s_{B}^{a} p,

(2.8)

and hence

S_{B}^{a} = \frac{S_{B}}{1 + V_{e}^{s} p} .

Before proceeding further, we need some terminology. By a “switch” we mean a Heaviside function, i.e. a function H(x) which is zero if x < 0 and one if x > 0.¹ In the single cell, a switch is either on or off. In a population model such as this one, the transition from “off to “on” is actually smooth and “smooth” versions of the Heaviside function were employed in the simulations. For example, one could use $H (x, K) = \max {x^{m} / (K + | x |^{m}), 0)}$ where K > 0 and m is an odd positive integer so that $H (K^{1 / m}, K) = 1 / 2.$ Alternatively H(x,m) = max{0, 1 – exp(–mx)} and H(ln 2/m,m) = 1/2.

The Heaviside function is only defined “almost everywhere”. Typically, one splits the difference and sets H(0) = 1/2.

We turn to the action of the growth factor. We consider not only the action of the transcription factor Sp1 on GF synthesis but also the degradation of growth factor during the course of the cell synthesis of MMP-1. The synthesis and degradation of growth factor is assumed to be regulated by Sp1 according to the discussion in the appendix. Thus we write:

\frac{\partial v}{\partial t} = \frac{k_{v} Y (x, t)}{K_{v} + Y (x, t)} (H (s_{B}^{a}) \frac{η_{B} (x, t)}{N_{B}} + H (s_{T}) \frac{η_{T} (x, t)}{N_{T}}) - μ_{v} v (x, t) .

(2.9)

Clearly, as the concentration of p53 falls, more inhibited Sp1 in benign cells is converted to the active form. This causes an increase in growth factor production.

A simplified scenario for the GF role in the production of MMP-1 might be the following: Once the cell has expressed a molecule of GF and released it to the ECM, the molecule binds to a GF receptor (R) and initiates signaling via a MAP-kinase pathway to induce transcription of the MMP-1 gene with the eventual translation to create the protease. A simplified mechanism for this is given in [44]. That mechanism, in its turn, is a bit more complicated than we need for our purposes here. Instead we use:

\begin{matrix} V + R ⇄_{k_{o f f}^{r}}^{K_{o n}^{r}} {V R} \\ Y + {V R} \overset{k_{m}^{r}}{\to} n M + V^{'} + R \end{matrix}

(2.10)

where V′ denotes the products of GF degradation. Michealis-Menten kinetics for this system yields, assuming V is in excess so that the concentration of the intermediate {VR} is constant:

\frac{d [M]}{d t} = \frac{n k_{o n}^{r} k_{o f f}^{r} [V] [Y] [R]}{k_{o f f}^{r} + k_{m}^{r} [Y]}

The number n is the cell amplification factor for MMP-1. It is thought to be fairly large and is a consequence of the amplification properties of the MAP-kinase signaling pathway and transcriptional response to growth factor. It is not constant, but depends upon the local concentration of growth factor, first increasing and then decreasing with growth factor concentration [61].

Likewise, there is a corresponding loss of growth factor that follows from the kinetics (2.10), namely

\frac{d [V]}{d t} = \frac{k_{o n}^{r} k_{o f f}^{r} [V] [Y] [R]}{k_{o f f}^{r} + k_{m}^{r} [Y]}

In the mechanism (2.10), we need to distinguish carefully between the receptors on the tumor cells and the receptors on the benign cells. Relating receptor density to cell density, and taking into account the transcription of GF as well as its decay and conversion to MMP-1, we may write:

\frac{\partial v}{\partial t} = \frac{k_{v} Y (x, t)}{K_{v} + Y (x, t)} (H (s_{B}^{a} (x, t)) \frac{η_{B} (x, t)}{N_{B}} + H (s_{T} (x, t)) \frac{η_{T} (x, t)}{N_{T}}) - μ_{v} v (x, t) - \frac{K_{c a t}^{r} v (x, t) Y (x, t)}{K_{m}^{r} + Y (x, t)} (\frac{η_{B} (x, t)}{N_{B}} + \frac{η_{T} (x, t)}{N_{T}}) .

(2.11)

In turn, growth factor induces the cellular expression of MMP-1:

\frac{\partial m}{\partial t} = Y (x, t) \frac{n K_{c a t}^{r} v (x, t)}{K_{m}^{r} + Y (x, t)} (\frac{η_{B} (x, t)}{N_{B}} + \frac{η_{T} (x, t)}{N_{T}}) - μ_{m} m (x, t) .

(2.12)

We have neglected the diffusion of growth factor in (2.11) because its half life is short in tissues.

Next we consider the rate equation for collagen/fibronectin degradation. It has been established experimentally that some tumor cell types have a limited ability to express fibronectin relative to benign cells [16, 54, 24, 64]. We assume that this rate is a small percent, ∈ _f of the rate of production of benign cells. Moreover, the molar rate of production of collagen/fibronectin must not only depend upon the local cell density but also upon the local concentration of resources, up to some saturable limit. To reflect these ideas, we write ²:

The chemical mechanism for collagen degradation is M + F ⇆ {MF} → M+F′ where F′ denotes the products of matrix degradation. The Michealis-Menten hypothesis yields the second term in (2.13). The logistic type term f(1 – f/f_M) must be modified to reflect cellular expression of matrix protein and hence the inclusion of the normalized cell concentrations. This can be made rigorous using kinetics but limitations of length preclude us doing so here.

\begin{matrix} \frac{\partial f}{\partial t} = \frac{(4 / T_{f} (x, t))}{K_{f} + Y (x, t)} f (x, t) (1 - \frac{f (x, t)}{f_{M}}) \\ (\frac{η_{B} (x, t)}{N_{B}} + ε f \frac{η_{T} (x, t)}{N_{T}}) \\ - \frac{K_{c a t}^{f} m (x, t) f (x, t)}{K_{m}^{f} + f (x, t)} . \end{matrix}

(2.13)

We have written 4/T_f as the time constant at saturation for convenience. In this form, at saturation and low concentrations of f, the doubling time for f is T_f ln 2/4. The factor of four is included since then 4x(1 – x) has a maximum value of unity.

Cell Movement

In order to track the movement of cells consider again the kinetics that follow from (2.10). The law of mass action when applied to the rate equation for the receptor density would yield ∂[R]/∂t = 0 in consequence of the Michealis-Menten hypothesis that [RV] is nearly constant. However, the issue is somewhat more complicated in this case because: (1) the receptor distribution is tied to the cell membrane movement; (2) the cell movement is dependent on the local concentrations of growth factor, enzyme and ECM protein in a chemotactic (and chemokinetic) as well as haptotactic manner rather than being completely random; and (3), the cells themselves undergo mitosis and apoptosis (programmed cell death). However, we can relate cell density to receptor concentration via the relationship [R](x, t) = cN(x, t) where the left hand side has the units of micro molarity (micro moles per liter) say while the cell density N is expressed in cells per liter so that the constant has units of micro moles per cell. Then a rate for [R] is proportional to a rate for N and we can consider cell movement rate laws in their own right coupled to the protein movement laws.

Suppose, for the moment, that N is the density of one of the two cell types under consideration here. From the continuity equation in the absence of sources or sinks (mitosis and cell death)

\partial_{t} N = - \vec{\nabla} . \vec{J}

(3.1)

where $\vec{J}$ is the local flux of cell density at (x, t).

Fick's law is then modified so that the flux of cell density will be influenced by the gradient in protease density m(x, t) (Chemotaxis), and gradients in collagen density f(x, t) (haptotaxis). These notions mean that each biochemical species influences the flux of the cells through a term that is proportional to the cell density and that depends on the gradient of the species, vis:

\vec{J} = D (N) {- \vec{\nabla} N + [M (m, f) \vec{\nabla} m + F (m, f) \vec{\nabla} f] N}

(3.2)

where M, F are some phenomenological functions of (m, f). These functions, sometimes called the chemotactic sensitivity functions, determine the influence of the specific species on the flux of cell densities. For example, where M > 0 the gradient of protease opposes the cell density gradient while where M < 0, it assists that gradient. If one makes the assumption that the vector $(M, F) = \nabla [T (m, f)]$ where ▽ is now the gradient in the variables (m, f), for some potential function T³, writes $T (m, f) \equiv$ = exp(T(m, f)) and uses this in the resulting flux vector, then the cell movement equation (without sources and sinks) takes the form

If $\vec{c u r l} [(M, F)] = \vec{0},$ such a function exists. In particular if M depends only on m and F only upon/this is always true.

\frac{\partial N}{\partial t} = \nabla . {D_{N} [N \nabla 1 n (\frac{N}{T (m, f)})]} .

(3.3)

This is perhaps a strange way to write the equation for chemotactic cell movement. However, the meaning is clear. It says that near steady state, “N should follow $T {(m, f)}^{"},$ an observation of Hans Weinberger (private communication).

The function $T (m, f)$ is called the probability transition function (PTF). We can understand its meaning at the intuitive level if we consider the stationary (time independent) version of the equation. Suppose N(x), m(x), f(x) are smooth functions connected by the relationship $N (x) = λ T (m (x), f (x))$ for some constant λ. Then it is clear that N(·) must be a stationary solution of (3.3) and we say that “the solution N follows $T$ ”.

Another way of looking at (3.3) is from the point of view of “cellular free energy”. That is, one views $T$ as a correction factor for chemotaxis/haptotaxis in the associated “cellular potential” here corresponding to the cell density. The “cellular potential” is defined as ρ(N) = – D ln(N/N₀) where N₀ is the “concentration” in some reference state. Without the correction factor, using Fick's first law, N_t = ▽J where J = – DN▽ρ(N) we obtain the ordinary diffusion equation for cell density. We can think of the cellular potential in much the same way as the chemist thinks of chemical potential for the isothermal change in free energy. Then we view $N / T (m, f)$ in much the same way as a chemist thinks of fugacity (mole fraction times pressure) at zero pressure in non-ideal gas dynamics or activity (which replaces ion concentration) in ionic solutions at finite dilution. See books on chemical thermodynamics such as [68] for more information. (The authors thank James Keener for bringing to our attention this interpretation of the appearance of the logarithm in (3.3).)

The choice of $T (m, f)$ is phenomenological. Other possible dependencies for the probability transition rate function $T$ are certainly possible. For example, $T (m, v, f)$ is a possibility. Any other external biochemical variable for which it is known that cells respond in a chemotactic/haptotactic manner may also be added to the argument list. For example it has been observed that endothelial cells will move up a protease gradient, but if the concentration of protease is too high, the enzyme will kill the cell, causing the cessation of movement. Likewise, EC cannot move through a matrigel bed if the density of matrigel is too high nor can they move along a surface unless there is some ECM protein such as collagen or fibronectin onto which it may attach its pseudopodia. Experimental evidence is given in [10]. These two qualitative observations suggest that this function should be biphasic in each variable (i. e. first increasing and then decreasing in each variable). The precise form we use in our simulations is described in Appendix B.

When we have two cell types present whose cell densities are denoted by η _B , η _T , say, (3.3) must be replaced by a system. One question that arises naturally is how to model the preference for cells of different types to fill regions vacated by dead cells. Generally speaking, as cells move, they consume the resources that are responsible for supplying the energy for cell movement. As this energy is exhausted, cell movement slows.

Additionally, the cells will behave differently under “crowding” conditions. That is, we expect the clumping of benign cells to slow their movement and push them into G₀, whereas tumor cells, which do not enter G₀, are more aggressive and not as influenced by the effects of crowding. In both cases, we expect movement to decrease as cell concentration rises, but we expect the rate of decrease of movement to be more rapid for benign cells than for tumor cells. Although we did not use this in our simulations, the model allows for this possibility. These ideas are much broader than we employ here [52]. Here we simply include the cell densities themselves as a part of the “chemotactic” sensitivity function and write

\begin{array}{l} \frac{\partial η_{B}}{\partial t} = \nabla . [D_{ℬ} (η_{B,} η_{T}) η_{B} \nabla 1 n (\frac{η_{B}}{T_{B} (η_{B}, η_{T}, m, f)})], \\ \frac{\partial η_{T}}{\partial t} = \nabla . [D_{T} (η_{B}, η_{T}) η_{T} \nabla 1 n (\frac{η_{T}}{T_{T} (η_{B}, η_{T}, m, f)})] . \end{array}

(3.4)

where we have omitted mitotic and apoptotic effects for the moment.

We next consider the issues of mitosis and apoptosis (sources and sinks) that must be included in both cell movement equations. These are viewed as forcing terms and are somewhat easier to explain. If we assume that the quantities of growth factor produced are close to the background rates, then cell proliferation rates may be taken to be nearly independent of growth factor. However this is not always true, especially in the case of angiogenesis. See [61] and [43] for an application of this fact.

It is known that p53 is a regulator of cell mitosis and proliferation. In particular, it enhances apoptosis and inhibits mitosis [50]. We make the assumption here that the distinguishing features between tumor cells and benign cells are the following:

A drop in p53 concentration below a critical value in benign cells results in the change of type (a mutation) of such cells to the malignant state. It is assumed that this occurs at a rate proportional to the benign cell density and is controlled by a switching mechanism. That is, the rate is proportional to H(p_c —p)η_B. The proportionality constant is denoted by λ _tr . Here p_c is a background critical value required to maintain a steady state population consisting only of benign cells in the region in question. The transfer is one way only; it allows for the growth of tumor cells at the expense of benign cells by permitting the conversion of the latter to the former through a modification in their respective growth rates.

The apoptosis to mitosis ratio of benign cells is higher than the corresponding ratio for tumor cells since the latter do not express p53 and the former do. (See [59] for some simple mathematical models of tumor growth based upon differing mitotic and apoptotic rates for tumor cells versus benign cells. See also the discussion in [7, 32].)

Both D_B, D_T are constant.

For ^* = B, T, the sensitivity functions $T_{*} (η_{B}, η_{T}, m, f)$ can be written as products $T_{*} (m, f) ψ_{*} (Z)$ where Z = η _B /N_B + η _T /N_T is the volume fraction of a given region occupied by cells. From physical considerations, we have 0 ≤ Z ≤ 1. (A simple formal argument using the logistic equation for total cell growth, namely dZ/dt = cZ(1 – Z) convinces us that if 0 < Z(0) < 1, then Z(t) remains in this interval. In the general case one can argue from the maximum principle.) In reality, Z ≤ Z_M < 1 where 1 – Z_M is the volume fraction of the extracellular region between the cells. However, we shall ignore this volume fraction in the model.

Based on the above considerations, the cell movement equations take the form:

\begin{array}{l} \frac{\partial η_{B}}{\partial t} = D_{ℬ} \nabla . [η_{B} \nabla 1 n (\frac{η_{B}}{T_{B} (m, f) Ψ_{B} (Z)})] \\ \begin{matrix} + [H (v - v_{B}) \frac{λ_{B} Y}{κ_{B} + Y} \\ (1 - \frac{η_{B}}{N_{B}} - L_{T} \frac{η_{T}}{N_{T}}) - μ_{B}] η_{B} \end{matrix} \\ - λ_{t r} H (p_{c} - p) η_{B'} \\ \frac{\partial η_{T}}{\partial t} = D_{T} \nabla . [η_{T} \nabla I n (\frac{η_{T}}{T_{T} (m, f) Φ_{T} (Z)})] \\ \begin{matrix} + [H (v) \frac{λ_{T} Y}{κ_{T} + Y} \\ (1 - \frac{η_{T}}{N_{T}} - L_{B} \frac{η_{B}}{N_{B}}) - μ_{T}] η_{T} \end{matrix} \\ + λ_{t r} H (p_{c} - p) η_{B} . \end{array}

(3.5)

The factor H(v - v_B) in the first of (3.5) says that when v > v_B benign cells proliferate, while when the growth factor is below this level, they cannot. The factorH(v) in the second of (3.5) has a similar meaning for tumor cells (with v_T = 0). These two factors reflect the fact benign cells are generally quiescent (in G₀) whereas tumor cells are usually constantly passing through the cell cycle. The threshold v_B is intended to model the lower proliferative response of benign cells to growth factor than that of tumor cells. Unfortunately, the number v_B is unknown. Therefore, by taking v_B = 0, in the model, we are giving the benign cells a greater survival advantage than they would ordinarily be expected to possess. The qualitative form of the results we give below will be unchanged by taking v_B > 0. The factors λ_BY/(K_B + Y) and λ_T Y/(K_T + Y) are included to reflect the idea that when nutrient levels are low, cell population production is low, but, as the nutrient levels rise, the population growth rates approach a saturation level given by the constants λ _B , λ_T.

The factors L_T, L_B are included as in [27, 28], We refer to them as mitosis inhibition or apoptosis enhancement factors as they can be viewed as either inhibiting the former or enhancing the latter. (We have written the equations as though they inhibit the mitosis.) In general, we expect that L_T ≥ L_B, that is tumor cells influence the benign cell net proliferation rate much more than vice-versa. This clearly gives the proliferation advantage to tumor cells. However, it does not imply that under all circumstances, the tumor cells will take over the region of interest from the benign cells.

Mass Balance Equations

In this model, nutrients function directly to support cell production of p53, Sp1, GF and MMP-1. These proteins also decay and are transported out of the region D by molecular diffusion (in the absence of a vasculature). Additionally, tumor and benign cells are undergoing mitosis and apoptosis in this region while collagen/fibronectin is being generated by the cells and degraded by MMP-1.

For example, in the rate equation for p53, resources are being consumed at a rate dictated by the first term on the right of (2.1) while the products of protein decay are being generated by the second term on the right of that equation. At steady state, as for example, in the healthy tissue, these two terms agree, otherwise they do not.

We have the mass rate loss for the resources:

where $S$ is the (positive) coefficient of Y in the preceding line. The constants l_B, l_T are conversion factors that permit us to convert the units of cells per unit volume to molecular weight per unit volume. That is, the constants $C_{B} \equiv N_{B} ℓ_{B} λ_{B}$ and $C_{T} \equiv N_{T} ℓ_{T} λ_{T}$ may be regarded as the maximum velocities of resource conversion into mass. These products are expected to be quite large, in any case much larger than the other velocities at saturation such as k_p, k_v, k_s, k_m etc. We have written the diffusion term, ▽·[D_Y ▽_Y] in this generality to indicate that the diffusion factor D_Y need not be constant.

If we prescribe the level of nutrients on the boundary of D and a given initial concentration of nutrients in D, with time, we expect the level of nutrients in the tissue to be low in regions “far” from the boundary. If the level of nutrient supply is sufficiently low, the tumor cells cannot proliferate except near the edge of the tumor. This is what is meant by “diffusion limited” tumor growth. This is the equation that controls the nutrient level in the tissue region, D.

Coupled to this should be an equation that reflects the diffusion of the products of cell metabolism in the tissue. Clearly the products of collagen degradation are among these products. Likewise, we should include the degradation products of p53, Sp1, MMP-1 and GF. However, since no further biochemical use of these waste materials is made in this article, we omit a detailed discussion of them.

Boundary and Initial conditions

We begin this section with a discussion of the boundary conditions. Since equations (2.1), (2.2), (2.6), (2.7), (2.11), (2.12), (2.13) are all ordinary differential equations in time, no boundary conditions are needed for them. On the other hand, (4.1), and (3.5) require boundary conditions. We take

\begin{array}{r} D_{ℬ} η_{ℬ} n . \nabla 1 n (\frac{η_{B}}{T_{B} (m, f) Ψ_{B} (z)}) = 0, \\ D_{T} η_{T} n . \nabla 1 n (\frac{η_{T}}{T_{T} (m, f) Ψ_{T} (Z)}) = 0, \\ D_{Y} n . \nabla Y + r_{y} [Y (x, t) - Y_{b}] = 0 \end{array}

(5.1)

where Z is again the local volume fraction occupied by the cell types and where it is understood that x ∈ ∂D, (the boundary of D) and where n denotes the outward directed unit normal to the boundary of D. The first two conditions say that the total flux of either cell type from the given region is zero. The flux constant r_y needs to be determined from phenomenological considerations. The constant Y_b is a prescribed level of nutrients on the outer wall.

We turn to a discussion of the initial conditions. Recall that a mutation is said to occur in a region D′ ⊂ D over at time interval I if there has been a loss of p53 expression in the set D'xI. That is, equation (2.1) is replaced by (2.2). The function Ψ is called the loss of function coefficient and the strength of the loss of function is defined in (2.3). We start such a perturbed system with a stationary solution of the system of equations for which there are no tumor cells present and no loss of function. Notice that we are not making any change in the initial values for either the proteins or the cell type densities. The only changes to be made are in the dynamics wherein the p53 equation takes the mutated form (2.2). The cell type equations are modified in that cell type change is driven only by the transfer rate λ_tr.

Dynamics of the Modified Lotka-Volterra System

In order to understand what might be expected from the full system, we consider the set of stationary solutions of (3.5) in the spatially homogeneous case for constant values of v, Y, f m. When the transfer coefficient Δ _tr = 0, we are led to consider a Lotka-Volterra competition model of the form:

\begin{array}{l} \frac{d x}{d t} = λ_{B}^{'} (x_{B} - x - L_{T} y) x, \\ \frac{d y}{d t} = λ_{T}^{'} (y_{T} - y - L_{B} x) y, \end{array}

(6.1)

when we freeze the coefficients that depend on the switches and resources.

A complete discussion of the dynamical behavior of the solutions of this system to be found in [46]. One is tempted to believe that the solutions of the full system might behave like those for the system of ordinary differential equations. However there are two very good reasons why this may not be the case.

First, the presence of the free energy terms (i.e. the terms involving cell movement, both random and chemotactic or haptotactic), as well as the fact that the various parameters may depend upon the other variables, does not make such an expectation a rigorous statement. (As a well known example of how partial differential operators can radically alter the behavior of solutions of systems of ordinary differential equations, we refer the reader to a discussion of Turing instability [46].) See, for other examples, [11, 18, 21, 22, 19, 20, 34, 38] where studies have been made of related systems in which the Lotka-Volterra terms (the right hand sides in (6.1)) function as reaction terms in a system of reaction diffusion equations of competition type.

Secondly, the system of ordinary differential equations that corresponds to the spatially homogeneous cell equations can be reduced to the following form when we freeze the coefficients that depend on the switches:

\begin{array}{l} \frac{d x}{d t} = λ_{B}^{'} (x_{B} - x - L_{T} y) x - Λ x, \\ \frac{d y}{d t} = λ_{T}^{'} (y_{T} - y - L_{B} x) y + Λ x . \end{array}

(6.2)

This reduces to the Lotka-Volterra system when Λ = 0.

The “constants” $λ_{B}^{'}, λ_{T}^{'}, x_{B}, y_{T}$ will depend upon $v, Y, μ_{B,} μ_{T,} λ_{B,} λ_{T,}$ , while Λ will depend upon λ _tr and p. The three relevant ordinary equations from our model take the form:

\begin{array}{l} \frac{d p}{d t} = \frac{k_{p} Y}{K_{p} + Y} \frac{η_{B}}{K_{B}} - μ_{p} p, \\ \frac{d η_{B}}{d t} = (λ_{B}^{'} H (v) \frac{y}{Λ_{B} + y} (1 - \frac{η_{B}}{N_{B}} - L_{T} \frac{η_{T}}{N_{T}}) - μ_{B}) \\ η_{B} - λ_{t r} H (p_{c} - p) η_{B}, \\ \frac{d η_{T}}{d t} = (λ_{T}^{'} H (v) \frac{y}{Λ_{T} + y} (1 - \frac{η_{T}}{N_{T}} - L_{B} \frac{η_{B}}{N_{B}}) - μ_{T}) \\ η_{T} + λ_{t r} H (p_{c} - p) η_{B} . \end{array}

In order to understand how p affects the dynamics, we use smooth forms of the Heaviside function and consider the case of no loss of function, i.e. Ψ ≡ 0. To put this system into standard form, we freeze the values of Y, v, let P_n be some normalizing constant for p53 and set

\begin{array}{l} x = \frac{η_{B}}{N_{B}}, y = \frac{η_{T}}{N_{T}}, z = p / P_{n}, \\ λ_{1} = λ_{B} H (v), K_{b v} \frac{Y}{Λ_{B} + Y}, μ_{1} = μ_{B}, \\ k_{1} = L_{T}, x_{0} = 1 - μ_{1} / λ_{1}, \\ λ_{2} = λ_{T} H (v), K_{t v} \frac{Y}{Λ_{T} + Y}, μ_{2} = μ_{T}, \\ k_{2} = L_{B}, y_{0} = 1 - μ_{2} / λ_{2}, \\ λ_{3} = \frac{(k_{p} / P_{n}) Y}{K_{p} + Y}, μ_{3} = μ_{p}, g (z) = λ_{t r} h (z) \end{array}

where h(z) is a smooth approximation of H(p_c/P_n–z). In terms of these quantities the system becomes

\begin{array}{l} x' = F (x, y, z) = λ_{1} x (x_{0} - x - k_{1} y) - x g (z) \\ y' = G (x, y, z) = λ_{2} y (y_{0} - y - k_{2} x) + x g (z) \\ z' = H (x, y, z) = λ_{3} x - μ_{3} z \end{array}

(6.3)

If the transfer rate λ _tr = 0, so that g = 0, then the first two equations decouple from the third equation, yielding a planar system for x and y which is a Lotka-Volterra competition model with four rest points: i) (0, 0), ii) (x₀, 0), iii) (0, y₀) and iv) (x₁ y₁), the point of intersection between the lines x₀ – x – k₁y = 0 and y₀ – y – k₂x = 0. The point (x₁, y₁) is not always physically relevant. Depending on k₁, k₂ it may not be in the first octant, or it may not satisfy x₁ + y₁ < 1. The dynamics of the system also varies with k₁, k₂. A thorough discussion can be found in [46]. In the case where the parameters are such that x₀ = y₀, k₁ = k₂ and k₁k₂ > 1, The points (x₀, 0) and (0, y₀) are asymptotically stable. These two, as well as a third, (x₁, y₁) are physically relevant (i.e. x₁ > 0, y₁ > 0, x₁ + y₁ ≤ 1). However, (x₁, y₁) is unstable.

The rest points of the full 3-dimensional system (6.3) satisfy F = G = H = 0. Since H = 0 if and only if z = αx, where α = λ₃/μ₃, the 3-d system can be reduced to a planar system for x and y; namely F (x, y αx) = G(x, y αx) = 0. This system can be analyzed graphically by plotting the curves F (x, y αx) = 0, G(x, y αx) = 0 and looking for points of intersection. For the choice of parameters and smooth approximation to the Heaviside function used in the simulations, this approach leads to Figure 2.

Figure 2

The null dines of the reduced planar system F(x,y,αx) = 0 (solid curves) and G(x,y,αx) = 0 (dashed curves). The rest points are the intersection points of a solid curve with a dashed curve. Notice that there is one with a negative y coordinate.

Figure 2 shows there are now 5 rest points, one of which is not phyically relevant. The solutions can be found numerically and a linear stability analysis can also be performed numerically. Table 2 gives the results of this process for the physically relevant rest points. These results show that for the reduced system (6.3) with smooth approximations to the Heaviside function, there can be an equilibrium state with co-existing benign and tumor cell densities that is unstable!

Table 2.

Numerically determined rest points of physical relevance (rounded to 3 decimal places) for the system (6.3) and their stability properties.

x	y	z	Dynamics
0.000	0.000	0.000	unstable, 1-Dunstable manifold
0.000	0.986	0.000	asymptotically stable
0.567	0.247	0.573	unstable, 1-Dunstable manifold
0.986	0.000	0.996	asymptotically stable

Steady State Solutions

Before we can perform simulations to explore the effects of modeling mutational loss of function in p53, steady states of the system must be determined. Since we cannot expect spatially homogeneous steady states, equilibrium solutions must be determined computationally, rather than algebraically. To do this, we use a shooting strategy which relies on asymptotic stability. That is, we choose initial conditions for the system and then let the system evolve. If the initial conditions are nearby an asymptotically stable steady state, then the evolving time-dependent solution will converge to the nearby equilibrium solution. This raises the question of how to choose appropriate initial conditions, which we discussed in the first subsection. Essentially we use constant solutions of a perturbed system. These values also serve as normalizing values for some of the variables. A second question arises along the way, that of an appropriate boundary value Y_b. This is discussed in the second subsection. Simulation results are also presented.

Initial Conditions and Normalizing Constants

In this subsection we treat Y_e as a parameter and seek spatially constant solutions of (2.1), (2.6), (2.7), (2.11), (2.12), (2.13), (3.5) and the following perturbation of (4.1):

\begin{array}{l} \frac{\partial Y}{\partial t} = \nabla . [D_{Y} \nabla Y] + Y_{r} - Y (x, t) \\ S (Y, f, s_{B}^{a}, s_{T}, v, η_{B}, η_{T}) \end{array}

(7.1)

where $S$ is defined implicitly in (4.1) and where Y_r is to be chosen so that there is a constant stationary solution for the system (2.1), (2.6), (2.7), (2.11), (2.12), (2.13), (7.1), and (3.5).

To determine this solution, we set the right hand sides of (2.1), (2.6), (2.7), (2.11), (2.12), (2.13), (3.5) to zero, assume that the variables in the resulting equations are constants, set all the Heaviside switches to unity except that we set H(p_c – p) = 0 in order to indicate that there is no transfer of cells from the benign state to the malignant state at p = p_e > p_c where p_e is an equilibrium value to be determined as explained below.

We use the underlying assumption that $η_{B}^{e} > 0$ and $η_{T}^{e} > 0$ Therefore $s_{T}^{e} = 0.$ = 0. For the benign variable we obtain a normalized solution, $η_{B}^{e} / N_{B},$ , as

N_{B}^{e} = \frac{η_{B}^{e}}{N_{B}} = 1 - \frac{μ_{B}}{Y_{B}^{e}}, (Y_{B}^{e} = \frac{λ_{B} Y_{e}}{κ_{B} + Y_{e}}) .

(7.2)

In order to have $N_{B}^{e} > 0$ we need $μ_{B} < Y_{B}^{e} .$ This is an inequality in Y_e that is equivalent to $Y_{e} > μ_{B} λ_{B} / (λ_{B} - μ_{B}) \equiv Y_{B}) .$ . It is also clear that $N_{B}^{e} \to 0^{+},$ as $Y_{e} \to Y_{B}^{+} .$

From the remaining equations we obtain:

\begin{matrix} p_{e} = \frac{Y_{p}^{e} N_{B}^{e}}{μ_{p}}, (Y_{p}^{e} = \frac{k_{P} Y_{e}}{K_{p} + Y}), \\ s_{e} \equiv s_{B}^{e} = \frac{Y_{s}^{e} N_{B}^{e}}{μ_{s B}}, (Y_{s}^{e} = \frac{k_{s} Y_{e}}{K_{s} + Y_{e}}), \\ s_{e}^{a} = \frac{s_{e}}{1 + v_{e}^{s} p_{e}} . \end{matrix}

(7.3)

Let

Y_{v}^{e} = \frac{k_{v} Y_{e}}{K_{v} + Y_{e}} and Y_{r}^{e} = \frac{K_{c a t}^{r} Y_{e}}{K_{m}^{r} + Y_{e}} .

Then for v_e and m_e, it follows that

v_{e} = \frac{Y_{v}^{e} N_{B}^{e}}{μ_{v} + Y_{r}^{e} N_{B}^{e}}, and m_{e} = \frac{n_{m} v_{e} Y_{r}^{e} N_{B}^{e}}{μ_{m}} .

(7.4)

To determine a constant solution f_e > 0 of (2.13), we set $F_{e} = f_{e} / f_{M}, {\tilde{K}}_{c a t}^{f} = K_{c a t}^{f} m_{e} / f_{m}, and {\tilde{K}}_{m}^{f} = K_{m}^{f} / f_{m} .$ Depending on parameter values there may or may not be a positive solution of the stationary equation for (2.13). If we set $x = F_{e}, α = Y_{f}^{e} N_{B}^{e}, β = {\tilde{K}}_{c a t}^{f},$ and $γ= {\tilde{K}}_{m}^{f} .$ It can be shown that if γ ≥ 1, there will be a positive solution if and only γ > β/α. It turns out that the stationary equation will have distinct real zeros if and only if β/α<(γ + 1)^2/4. In this case, the larger one is given by $x = \frac{1}{2} (1 - γ + \sqrt{{(γ+1)}^{2} - 4 (β / α)}) .$ If γ < 1 this root is positive. Since γ ≤ (γ + 1)²/4, it follows that the weaker requirement β/α < (γ + 1)²/4 is necessary and sufficient for x > 0 when γ < 1. If the constant Y_r is now determined by setting $Y_{r} = Y_{e} S (Y_{e}, f_{e}, s_{e}, 0, v_{e}, η_{B}^{e}, 0),$ then the above values give the desired constant solutions of (2.1), (2.6), (2.7), (2.11), (2.12), (2.13), (7.1), (3.5).

Remark 2

Notice that $p_{e}, …, η_{B}^{e}$ depend only upon the kinetic constants, the proliferation and apoptosis rates and Y_e, the level of resources.

The boundary value Y_b

The constant values found above are not steady states of the original system in which (4.1) appears rather than (7.1). In fact constant solutions are not possible. Clearly the distribution of resources cannot be spatially homogeneous, since $Y_{s} = Y S (Y, f, v, s_{B}^{a}, s_{T}, η_{B,} η_{T}) > 0,$ unless all variables are zero. However, there are spatially inhomogeneous equilibrium solutions of the system that are close to these constant values, provided that the flux of nutrients at the boundary diffuses through the region resulting in a distribution of resources approximating the distribution determined by Y_r. The key issue is to determine the proper boundary value for Y_b. This can be done analytically as follows.

Suppose that Y₀(x) is the resource component of an equilibrium solution of the model system (with Y_r = 0) in which there are no tumor cells present. Then Y₀ (x) satisfies

\begin{array}{l} 0 = \nabla . (D_{Y} \nabla Y) - Y_{s} (Y, f, v, s_{B}^{a}, 0, η_{B}, 0), x \in D, \\ 0 = D_{y} \frac{\partial Y}{\partial v} + r_{y} (Y - Y_{b}), x \in \partial D, \end{array}

(7.5)

where D is the region of interest, $\frac{\partial Y}{\partial v}$ denotes the outward normal derivative and $s_{B}^{a}, v, f, Y, η_{B}$ are the corresponding steady states of the full system. If Y₀(x) ≈ Y_e throughout D then we would expect their average values to be approximately the same. We use this observation to determine Y_b, and hence the flux at the boundary, by requiring

\frac{1}{| D |} \int_{D} Y_{0} (x) d x = Y_{e}, (| D | = \int_{D} 1 d x .)

Suppose that Ŷ (x) is the solution of (7.1) satisfying the homogeneous boundary condition $D_{Y} \frac{\partial Y}{\partial v} + r_{y} Y = 0, x \in \partial D .$ Since the equations (7.1) are linear and Y_b is constant, the principle of superposition yields Y₀(x) = Ŷ (x) + Y_b. Therefore $Y_{b} = \frac{1}{| D |} \int_{D} (Y_{0} (x) - \hat{Y} (x)) d x = Y_{e} - \frac{1}{| D |} \int_{D} \hat{Y} (x) d x .$ Thus, to determine Y_b, it is sufficient to determine the mean value of Ŷ(x). If the values of $s_{B}^{a}, v, f, Y, η_{B}$ are approximately the same as the normalizing values $s_{e}, v_{e}, f_{e}, Y_{e}, η_{B}^{e}$ then we should have $\hat{Y} (x) \approx {\hat{Y}}_{e} (x),$ where $\hat{Y_{e}} (x)$ satisfies $0 = \nabla . (D_{Y} \nabla Y) - Y_{s}^{e}, x \in D, D_{Y} \frac{\partial Y}{\partial v} + r_{y} Y = 0, x \in \partial D,$ where $Y_{s}^{e} = Y_{e} S (Y_{e}, f_{e}, s_{e}, 0, v_{e}, η_{B}^{e}, 0)$ is a constant. Clearly $Y_{s}^{e} = Y_{r}$ by the choice of Y_r that was made to obtain the normalizing constants.

In the case that $D = {x \in R^{3} : | | x | | < R}$ is the ball of radius R centered at the origin, then $\hat{Y_{e}} (x)$ can be explicitly computed to be $\hat{Y_{e}} (r) = \frac{Y_{s}^{e}}{3} (\frac{R}{r_{y}} + \frac{R^{2} - r^{2}}{2 D_{y}}) .$ Integration yields $\int_{D} {\hat{Y}}_{e} (x) d x = - \frac{Y_{s}^{e}}{3} (\frac{R}{r_{y}} \frac{4 π R^{3}}{3} + \frac{4 π R^{5}}{15 D_{Y}}),$ and hence $Y_{b} = Y_{e} + \frac{Y_{s}^{e}}{3} (\frac{R}{r_{y}} + \frac{R^{2}}{5 D_{Y}})$ as an appropriate resource level at the boundary.

Evolution to Equilibrium

All our computations were done with the system in non-dimensionalized form. The non-dimensionalization of the variables is given in Table 3. To computationally determine steady states we used a shooting strategy. We were able to determine an asymptotically stable steady state of the system, as predicted by the reduced system, that has a nearly constant benign cell population with no tumor cells present. It was possible to use the constant states found in section 7.1 as initial values and successfully evolve to this steady state. In terms of the non-dimensional variables many of the initial values were equal to one. However, a set of initial values that lies closer to this steady state is obtained if the values of p, η _B and η _T are chosen instead to be the (constant) solutions of the reduced system, with v = v_e, Y = Y_e fixed that are given in Table 5.

Table 3.

Nondimensionalized variables. Except for length and time scales, the normalizing constants in this table were computed in subsection 7.1.

Quantity	Variable	Dimensionless variable
time	t	τ = t/t
position	x	ŗ=X/L
p53	p(x,t)	p(ŗ,τ)=p(x,t)/p_e
Sp1	s_B(x,t)	s_B(ŗ,τ)=s_B(x,t)/s_e
active Sp1	s_B^a(x,t)	s_B^a(ŗ,τ)=s_B^a(x,t)/s_e
Sp1	s_T(x,t)	s_T(ŗ,τ)=s_T(x,t)/s_e
MMP-1	m(x,t)	M(ŗ,τ)=m(x,t)/m_e
GF	v(x,t)	V(ŗ,τ)=v(x,t)/v_e
collagen	f(x,t)	F(ŗ,τ)=f(x,t)/f_M,F_e=f_e/f_M
Nutrients	Y(x,t)	Y(ŗ,τ)=Y(x,t)/Y_e
benign cell fraction	η_B(x,t)	N_B(ŗ,τ)=η_B(x,t)/N_B
tumor cell fraction	η_N(x,t)	N_T(ŗ,τ)=η_T(x,t)/N_T
cell fraction total	η(x,t)=η_B(x,t)/N_B+η_T	N(ŗ,τ)=N_B(ŗ,τ)+N_T(ŗ,τ)

Table 4.

Numerical values used in simulations.

Equation	Constants	Values	References
	R	0.100(10¹)mm	simulated value
(2.1)	k_p	0.114(10⁰) h^–1	see notes
	K_p	0.500(10³) μM	simulated value
	μ_p	0.140(10¹) h^–1	see notes
(2.6),(2.7)	k_s	0.644(10^–1)h^–1	see notes
	K_s	0.500(10³) μM	simulated value
	$μ_{s_{B}}$	0.180(10¹)h^–1	see notes
	$μ_{s_{T}}$	0.180(10¹)h^–1	see notes
(2.8)	v_e	0.500(10²) (μM)^–1	simulated value
(2.11)	k_v	0.313(10⁰) h^–1	see notes
	K_v	0.500(10³) μM	simulated value
	μ_v	0.208(10¹) h^–1	see notes
	$K_{c a t}^{r}$	0.300(10¹) h^–1	simulated value
	$K_{m}^{r}$	0.500(10³) μM	simulated value
(2.12)	μ_m	0.416(10¹) h^–1	see notes
	n_m	0.800(10¹)	simulated value
(2.13)	T_f	0.180(10²) h	see notes
	K_f	0.500(10³) μM	simulated value
	f_M	0.140(10⁴) μM	see notes
	∈_f	0.100(10⁰)	simulated value
	$K_{c a t}^{f}$	0.170(10²) h^–1	simulated value
	$K_{m}^{f}$	0.700(10³) μM	simulated value
(3.5)	N_B	0.100(10¹³)l^–1	see notes
	D_B	0.360(10^–5) mm²h^–1	see notes
	λ_B	0.120(10^–1) h^–1	see notes
	κ_B	0.500(10³) μM	simulated value
	L_T	0.150(10¹)	simulated value
	μ_B	0.100(10^–3)h^–1	see notes
	λ_tr	0.750(10^–2)h^–1	see notes
	N_T	0.100(10¹³)l^–1	see notes
	D_T	0.360(10^–5)mm²h^–1	see notes
	λ_T	0.120(10^–1) h^–1	see notes
	κ_T	0.500(10³) μM	simulated value
	L_B	0.150(10¹)	simulated value
	μ_T	0.100(10^–3)h^–1	see notes
(4.1)	D_Y	0.500(10^–1)mm²h^–1	simulated value
	ℓ_B	0.200(10^–5) μM	simulated value
	ℓ_T	0.200(10^–5)	simulated value
(5.1)	r_y	0.100(10¹)mmh^–1	simulated value

Table 5.

Derived normalizing constants used in simulations. The choice of Y_e given in the table determines normalized benign cell density to be $η_{B}^{e} / N_{B} = 0.990, f_{e} = 0.139 (10^{4}) μ M .$ and the normalizing constants given in the table.

variable	value
Y_e	0.200(10⁴)μM
p_e	0.645(10^–1)μM
s_e	0.283(10^–1)μM
v_e	0.278(10^–1)μM
m_e	0.127(10⁰)μM

Figure 3 shows the evolution of the benign cell density, starting with the initial values of S_B, S_T, v, m, f and Y chosen as the constants from section 7.1 and the initial values of p, η _B and η _T chosen as the solutions of the reduced system (with all values non-dimensionalized). Spatial profiles at increments of 100 hours are superimposed in this plot. After a fairly rapid initial transient the profiles gradually approach an equilibrium. All other variables in the system were observed to approach their equilibrium profiles much more rapidly than the benign cell density. The equilibrium states of the full system are shown in Figure 4.

Figure 3

Spatial profiles of N_B (r,t) with t = t_i = iΔt, where Δt = 100 hours. As time increases the profiles monotonically approach an equilibrium solution profile.

Figure 4

The state of the system just prior to the start of the loss of function event. All variables have been non-dimensionalized. The system was initialized to its equilibrium configuration. During the first 100 hours of the simulation it remains in an equilibrium state which is shown at the end of this 100 hour period. The loss of function event was then initiated. The functions plotted are: a) growth factor, v (–), protease, m (- ·), p53, p (—); b) total benign Sp1, s_B (—), active benign Sp1, $s_{B}^{a}$ (- ·), tumor Sp1, s_T (–); c) benign cell density, η_B (—), tumor cell density, η_T (–); d) fibronectin, f (—), resources, Y (–).

Results and Conclusions

It is clear that the system of equations and boundary conditions is quite involved. It is perhaps worth taking a paragraph or two to explain what one might expect from the replacement of (2.1) by (2.2). The implementation of the loss of function is through equation (2.2) and is controlled by the loss of function coefficient. Suppose that the level of p53 falls below the critical level. Then there will be a small gain in tumor cell density. This will drive, through the equation for tumor induced Sp1, an increase in growth factor production and hence an increase in MMP-1 production as one sees from the first term in (2.12). This in turn leads to a degradation of matrix collagen/fibronectin as is indicated by the sink term in (2.13). If there are no tumor cells initially present, this small gain in tumor cell density will drive the production of tumor expressed Sp1. The tumor expressed Sp1 is not under the regulation of p53 and hence can induce an increase in the rate of growth factor production. The induced concentration gradients in these two proteins in turn stimulate the chemotactic movement of the two cell types.

As presented in (2.2), the loss of function is modeled by the expression 1 –Ψ (x,t). Suppose that Ψ(x, t) = Iϕ(x)ψ(t) where: i) ϕ(x) = 1 for x ∈ D′ and is zero elsewhere, ii) Ψ (t) = 1 for t₀ < t < t₀ + T and is zero elsewhere (for some t₀), and iii) ℐ ∈ [0,1] is an amplitude factor. Then T is the duration of the loss of function event, and ℐ is the intensity of the loss of function, with ℐ = 0 corresponding to no loss of function and ℐ = 1 being complete loss of function. The domain D′ can also be manipulated, but we will focus more on duration and intensity. For example, suppose D′ = B(r₀) is the ball about zero of a fixed radius $r_{0} < R .$ Then $S = 4 T π r_{0}^{3} ℐ / 3$ is the strength of the loss of function. When $S = 0 = ℐ$ there is no loss of function while when $S = 4 T π_{0}^{3} / 3,$ the loss of function is complete on B(r₀) × [t₀, t₀ + T].

The simulations were done using a spherical domain D, with the origin of the coordinate system at the point where loss of function occurs. Radial symmetry was assumed, so that all functions depended only on r = |x|, the distance to the origin. The domain was non-dimensionalized to the ball bounded by the sphere of radius R = 1. The domain D′ was taken to be a ball with a smaller radius. As with the other step functions appearing in the model description, smoothed versions of the functions ϕ(x) and ѱ (t) defining Ψ(x, t) were used in the simulations. In particular. we chose $ϕ (r) = \frac{1}{2} (1 + \cos (r - r_{0} + δ) π / (2 δ)),$ for |r – r₀| = δ, with ϕ(r) = 1 to the left of this interval and ϕ(r) = 0 to the right. Thus, D′ was a ball of radius r₀ + δ.

Basic Simulations

The model will demonstrate that the loss of p53 function can lead to the uncontrolled growth of a malignant avascular tumor if the loss is sufficiently strong, occurs for a sufficiently long time or involves a sufficiently large number of cells (as measured by the product of cell density times volume). On the other hand if any one of these three factors is sufficiently small, regardless of the intensity of the others, the loss of function will be harmless in the sense that any initial tumor formation will decay. The loss of function is manifested in the model by a decline in the mitosis rate for benign cells together with an increase in the mitosis rate for tumor cells, both rates having the same (maximum) magnitude λ_trη_B where η_B is the local density of benign cells.

In particular, the model demonstrates the existence of a critical intensity, ℐ_c such that if the intensity is below this number, any tumor formed during the time period during which loss of function is experienced will decay once the loss of function is no longer being experienced. On the other hand, if the intensity is larger than this critical intensity, then there is a critical time duration, T_ℐ, such that if the loss of function is experienced for any finite time longer than this critical time, the formed tumor will grow and take over the entire region. (Recall that this region is only one or two mm in diameter and that we are dealing only with solid tumor growth.) If the loss of function occurs only for a time duration smaller than this critical time, the solid tumor will eventually decay.

In this case, the model suggests the existence of an unstable tumor cell region coexisting with the benign cell region which will form if the loss of function is experienced for the critical duration. From a practical point of view, one can view the formation of such a tumor cell region as occurring if the duration time is only close enough to the critical time. For example if the doubling time (when T > T_ℐ) or half life (when (T < T_ℐ) is much larger than the remaining life expectancy of the individual, we can view the tumor to be in a quasi steady state with its surroundings.

Figures 4–6 show the result of a loss of function event that does not result in tumor growth. Figures 7–9 show the result of a loss of function event that does result in tumor growth. Both simulations were started with the system in equilibrium and allowed to evolve, maintaining equilibrium, for 100 hours. The equilibrium state is shown in Figure 4. After a very brief one hour transition the loss of function event begins. The intensity of the event is the same in both cases, ℐ = 0.5; that is a 50% loss of function. The parameters determining ϕ(r) were also the same in both cases: r₀ = 0.1, δ = 0.05. After T hours, another brief one hour transition to a system with no loss of function is made. These times are recorded in the figures, with the one hour transitions ignored. The durations were T = 200 hours in the first case and T = 400 hours in the second case.

Figure 5

The state of the system just after the end of the loss of function event. Functions plotted are the same as in Figure 4. The intensity of the event in the simulation was ℐ = 0.5, that is a 50% loss of function. The duration of the event was T = 200 hours. The state of the system at the end of the event is shown. For subsequent times the system evolves with no loss of function.

Figure 6

The state of the system in Figure 5 at a later time. Functions plotted are the same as in Figure 4. The system has continued to evolve and eventually returns to the original equilibrium state. An immediate time between the end of the loss of function event and full recovery is shown.

Figure 7

The state of the system just after the end of the loss of function event. Functions plotted are the same as in Figure 4. The intensity of the event in the simulation was ℐ = 0.5, that is a 50% loss of function. The duration of the event was T = 400 hours. The state of the system at the end of the event is shown. For subsequent times the system evolves with no loss of function.

Fgure 8

The state of the system in Figure 7 as the tumor cell population expands, causing the density of benign cells to fall dramatically in the region surrounding the location where the loss of function occurred. The spread of the tumor cell population and the decline of benign cell population progresses in the form of a traveling wave. Functions plotted are the same as in Figure 4.

Figure 9

The state of the system shown in Figures 7, 8 at a much later time. Functions plotted are the same as in Figure 4. The tumor cell population has almost completely displaced the benign cell population in the entire region.

Threshold and Critical Values

Given a level of intensity ℐ for a loss of function event, there may be a critical value T_ℐ of the time of duration such that if T > T_ℐ then the event results in tumor growth, while if T < T_ℐ then the event does not result in tumor growth. The simulation results presented in the previous subsection clearly suggest this is the case when ℐ = 0.5. It is also intuitively clear that if the intensity is extremely small then no significant tumor growth will occur, regardless of the duration of the loss of function event. In other words, there should be a threshold of intensity that must be exceeded before significant tumor growth can occur after a sufficiently long loss of function event.

The roles of ℐ and T can be reversed and one can ask, for a given time of duration T, whether or not there is a critical intensity ℐ_c. From this point on view, one would expect a threshold for the time of duration that must be exceeded for significant tumor growth to occur even when the intensity is as large as possible.

Threshold values for both ℐ and T are predicted by the model. In the model, a loss of function event results in a decrease in the production of p53. If the amount being produced does not fall below the critical value p_c then the transfer term, λ_tr H(p_c - p)η_B, appearing in the cell equations (3.5) will remain zero and no change in the dynamics should be expected. (However, there will be a subtle change in the growth factor and resource level that can change production levels slightly.) In order to have the level of p53 fall below the critical value, the loss of function event must have a sufficient intensity, and this results in thresholds for both ℐ and T.

The graph in Figure 10 shows the functional dependence of the critical time T_ℐ of duration on the intensity for loss of function events as predicted by the model. It also shows the threshold value for T at maximum intensity ℐ = 1. The threshold for ℐ is indicated by the vertical asymptote. The curve also depends on the size of the region D′ where the loss of function occurs. For the plotted curve, the same function ϕ(r) described above was used, with the parameters determining ϕ(r) being r₀ = 0.10, δ = 0.05. To determine points on the curve, a bisection technique was used. For example, for a given value of ℐ, two bracketing values T₀ < T₁ of the time of duration were found such that the system returned to equilibrium if the time of duration was T₀, or experience overwhelming tumor growth if the time of duration was T₁. Then the average of these times T_m = (T₀ + T₁)/2 was used as the duration for the loss of function event in a simulation experiment. If the system returned to equilibrium then T₀ was replaced by T_m, otherwise T₁ was replaced by T_m. Then the experiment was repeated.

Figure 10

The curve of critical points that indicates threshold values for tumor growth.

When the domain D′ over which the loss of function takes place increases in size the intensity or duration of the event is diminished. In Figure 11 two curves of critical points are plotted. The lower curve was determined using a domain D′ and function ϕ(r) determine by the parameters r₀ = 0.5, δ = 0.05, while the upper curve corresponds to r₀ = 0.10, δ = 0.05 as in Figure 10. Because mutations in one location are very rare, the domain of a single the loss of function event is expected to have roughly the volume of a single cell or two.

Footnotes

Biochemistry

In this section, using a kinetic approach to the molecular events involved in transcription and translation, we derive equations (2.1), (2.6), (2.7), (2.11), and (2.12). Because the actual events are quite complex, we take an overall approach to them. Our purpose is to give some justification for the simple relationship used in these equations from a kinetic point of view.

The sources of the nucleotides and amino acids involved in the assembly of the mRNA and proteins must ultimately be the blood and tissue.

For our purposes, a cell consists of two regions, the nucleus and the cytoplasm separated by a membrane. The transcription of a gene into RNA and its maturation into mRNA takes place in the former. The mRNA is transferred from the nucleus across the nuclear membrane into the cytoplasm. There the ribosome translates it, via the genetic code, into a protein. This protein may stay associated with the cell or it may be secreted into the extra cellular matrix (ECM) where it may perform several functions. For example, it might be a protease that degrades the ECM or a growth factor that binds to a receptor on the same or another cell (that need not be of the same type as the parent cell) to initiate the expression of other genes via some signal transduction pathway.

Clearly, the sequence of events involved in transcription and translation will involve biochemical mechanisms and processes that may be only incompletely understood quantitatively as well as qualitatively. In order to obtain the “starting” material from the ECM and plasma, one has to transfer these molecules across the cell membrane and, in the case of the nucleotides, across the nuclear membrane. This is typically accomplished by diffusion. There are transmembrane proteins that permit the transfer of such small molecules.

We take the point of view in all of this that the rate limiting step is the translation step in the cytoplasm: That is, the ribosome (Rb), the cytoplasmic mRNA (^cRn) and the tRNA (not shown) work to assemble the protein from the amino acids in the cytoplasm:

(Here n represents the number of molecules of the protein that are assembled during the life time of the mRNA.)

However, if we assume that the second of equations (A.1) is rate limiting and we apply the Law of Mass Action to this chemical equation we should expect that if 7 is in excess,

where [mRNA(P_G)] is the cytoplasmic concentration of mRNA that codes for P_G and where we have included an additional decay term μ[P_G] on the right of this last equation. The Heaviside factor H([T_G](f)) reflects the fact that it is the transcription factor, T_G, that initiates the mRNA(P_G) production. (The factor n is included in the constant K in this last equation.)

Because the available mRNA for each gene is proportional to the number of cells whose genes express it, we can write, for some constants k_G, μ _G

Depending on other underlying processes within the cell, this rate equation may be subject to further modification. For example, one might have two proteins, one of which functions as an enzyme for the degradation of the other. For example, if protein $P_{G_{1}}$ functions as a regulator for protein P_G with K_cat, K_m as kinetic constants, then the equation that replaces (A.2) will be

Remark 3. In (A.2), (A.3) the assumption was that the concentration of cell nutrients [Y] is in excess. This need not be the case, especially in situations that lead to tumor cell necrosis. Therefore, one may need to modify (A.2), (A.3) by replacing k_G by k_g Y/(K_m + Y), which would then reflect the fact that the rate of transcription is also proportional to the level of resources available:

The rate law for nutrient consumption must also be introduced. It will take the form:

where the sum is taken over all genes that are expressed as proteins. The term S_r includes both tissue sources of nutrients and diffusive effects. For example, if the distribution of nutrients is not spatially homogeneous, then

where Δ denotes the Laplace (diffusion) operator (Δ[Y] = div(grad [Y]). The term Y_r (x, y, t) includes any tissue generated nutrients that contribute to the cellular events described above. Boundary conditions for (A.7) will include capillary supplied nutrients that move from the vasculature into the ECM.

Appendix B.

The determination of appropriate biological constants from the literature is notoriously difficult. Not only are in vivo values often different from in vitro values, but very likely the former are unknown. When the constants for a given protein are known, the protein itself may function entirely differently in vivo than in vitro. That is to say, different cell types can respond to the same growth factors in different ways, and the response of a single cell type can be influenced by its environment, including the cells that surround it. Thus, the responses to a growth factor observed in vitro of a particular cell type grown as a mono culture are not necessarily the same as for the same cell type in tissues in vivo where its neighbors are not of the same cell type []. Thus, the task to find good values for the model is daunting, to say the least. We took a pragmatic approach, using available literature values where possible and trial values, when no values were available to us.

References

Adam

J.A.

, Mathematical models of prevascular spheroid development and catastrophe-theoretic description of rapid metastatic growth/tumor remission, Invasion Metastasis, 16(4-5): 247–67, (1996).

Adam

J. A.

, Mathematical models of tumor growth: from empirical description to biological mechanism Adv. Exp. Med. Biol, 537: 287–300, (2003).

Adam

J. A.

, and Maggelakis

S.A.

, Mathematical models of tumor growth. IV. Effects of a necrotic core, Math Biosci., 97(1): 121–36, (1989).

Araujo

R.P.

, and McElwain

D.L.S.

, A history of the study of solid tumour growth: the contribution of mathematical modelling Bulletin of Mathematical Biology, 66(5): 1039–1091, (2004).

Asthagiri

A.R.

, Nelson

CM.

, Horwitz

A.F.

, and Lauffenburger

D.A.

, Quantitative relationship among integrin-ligand binding, adhesion, and signaling via focal adhesion kinase and extracellular signal-regulated kinase 2, J. Biol. Chem., 274: 27119–27, (1999).

Band

, Dalai

, Delmolino

, and Androphy

E. J.

, Enhanced degradation of p53 protein in HPV-t and BPV-1 E6 immortalized human mammary epithelial cells, EMBO J., 12(5): 1847–52, (1993).

Berns

, Cancer: Two in one, Nature, 436–HAugust: 787–789, (2005).

Berthelot

, Muldoon

, rjkowitsch

, Hughes

, and McCarthy

J.E.G.

, Dynamics and porcessivity of 40S ribosome scanning on mRNAin yeast ,Molec. Biol, 51(4): 987–1001, (2004).

Bertsch

, Gurtin

M.E.

, Hilhorst

, and Peletier

L.A.

, On interacting populations that disperse to avoid crowding: The effect of a sedentary colony, J. Math. Biol, 19: 1–12, (1984).

10.

Bowersox

J.C.

, and Sorgente

, Chemotaxis of aortic endothelial cells in response to fibronectin Cancer Res., 4: 2547–2551, (1982).

11.

Burton

TA.

, and Hutson

, Permanence for nonautonomous predatorprey systems, Differential Integral Equations, 4(6): 1269–1280, (1991).

12.

Byrne

H.M.

, and Chaplain

M.A.

, Growth of nonnecrotic tumors in the presence and absence of inhibitors, Math Biosci., 130(2): 151–81, (1995).

13.

Calderon

C. P.

, and Kwembe

T.A.

, Modeling tumor growth Math Biosci., 103(1): 97–114, (1991).

14.

Casciari

J.J.

, Sotirchos

S.V.

, and Sutherland

R.M.

, Mathematical modelling of microenvironment and growth in EMT6/R0 multicellular tumour spheroids. Cell Prolif, 25(1): 1–22, (1992).

15.

Chen

C.Y.

, Byrne

H.M.

, King

J.R.

, The influence of growth-induced stress from the surrounding medium on the development of multicell spheroids, J. Math. Biol, 43(3): 191–220, (2001).

16.

Cizdziel

P.E.

, Hosoi

, Montgomery

J.C.

, Wiseman

R.W.

, and Barrett

J.C.

, Loss of a tumor suppressor gene function is correlated with downregulation of chondrocyte-specific collagen expression in Syrian hamster embrio cells, Mol Carcinog, 4(1): 14–24, (1991).

17.

Chaplain

M.A.

, and Sleeman

B.D.

, A mathematical model for the growth and classification of a solid tumor: a new approach via nonlinear elasticity theory using strain-energy functions, Math. Biosci., 111(2): 169–215, (1992).

18.

Cantrell

R.S.

, Cosner

, and Lou

, Multiple reversals of competitive dominance in ecological reserves via external habitat degradation J. Dynam. Differential Equations, 16(4): 973–1010, (2004).

19.

Cantrell

R. S.

, Cosner

, and Hutson

, Ecological models, permanence and spatial heterogeneity, Rocky Mountain J. Math, 26(1): 1–35, (1996).

20.

Cantrell

R.S.

, Cosner

, and Hutson

, Permanence in ecological systems with spatial heterogeneity Proc. Roy. Soc. Edinburgh Sect. A, 123(3): 533–559, (1993).

21.

Cantrell

R.S.

, and Cosner

, On the uniqueness and stability of positive solutions in the Lotka-Volterra competition model with diffusion, Houston J. Math., 15(3): 341–361, (1989).

22.

Cantrell

R.S.

, and Cosner

, On the effects of spatial heterogeneity on the persistence of interacting species, J. Math. Biol, 37(2): 103–145, (1998).

23.

Eppler

S.M.

, Combs

D.L.

, Henry

T.D.

, Lopez

J.J.

, Ellis

S.G.

, Yi

, Annex

B.H.

, McCluskey

E.R.

, and Zioncheck

T.F.

, A target-mediated model to describe the pharmacokinetics an hemodynamic effects of recombinant human vascular endothelial growth factor in humans, Cli. Pharm. Then, 72: 20–32, (2002).

24.

Fenhalls

, Geyp

, Dent

D. M.

, and Parker

M.I.

, Brest tumour cell induced down-regulation of type I collagen mRNA in fibroblasts, Brit. J. Cancer, 81(7): 1142–1149, (1999).

25.

Filippini

, Griffin

, Uhr

, Eppenberger

, Bonilla

, Cavalli

, and Soldati

, A novel flow cytometric method for the quantification of p53 gene expression, Cytometry, 31: 180–6, (1998).

26.

Gammack

, Byrne

H.M.

, and Lewis

C.E.

, Estimating the selective advantage of mutant p53 tumour cells to repeated rounds of hypoxia, Bull. Math. Biol., 63(1): 135–66, (2001).

27.

Gatenby

R.A.

, and Gawlinski

E.T.

, A reaction-diffusion model of cancer invasion, Cancer Research, 56: 5745–5753, (1996).

28.

Gatenby

R.A.

, and Gawlinski

E.T.

, The glycolytic phenotype in carcinogenesis and tumor invasion: insights through mathematical models Cancer Research, 63: 3847–3854, (2003).

29.

Gatenby

R.A.

Models of tumor-host interaction as competing populations: implications for tumor biology and treatment, J. Theor. Biol., 176(4): 447–55, (1995).

30.

Guiot

, Degiorgis

P.G.

, Delsanto

P.P.

, Gabriele

, and Deisboeck

T.S.

, Does tumor growth follow a “universal law”?, J. Theor. Biol., 225(2: 147–51, (2003).

31.

Hanahan

, and Folkman

, Patterns and emerging mechanisms of the angiogenic switch during tumorigenesis, Cell, 86: 353–364, (1996).

32.

Hermann

M.T.

, Brick

, Teruya-Feldstein

, Herbst

, Nilsson

J.A.

, Cordon-Cardo

, Cleveland

J.L.

, Tansey

W.P.

, and Lowe

S.W.

, Evasion of the p53 tumour surveillance network by tumour-derived MYC mutants, Nature, 463(11): 807–811,(2005).

33.

Hervagault

J.F.

, Ortoleva

P.J.

, Ross

, A plausible model for reversal. of neoplastic transformations in plants based on multiple steady states, Proc. Natl. Acad. Sci. U S A., 88(23): 10797–800, (1988).

34.

Hutson

, and Law

, Permanent coexistence in general models of three interacting species, J. Math. Biol., 21(3): 285–298, (1985).

35.

Jackson

T.L.

, and Byrne

H.M.

, A mathematical model to study the effects of drug resistance and vasculature on the response of solid tumors to chemotherapy Math. Biosci, 164(1): 17–38, (2000).

36.

Jagers

, The proportions of individuals of different kinds in two-type populations. A branching process problem arising in biology, J. Appl. Probability, 6: 249–260, (1969).

37.

Jones

A.F.

, Byrne

H.M.

, Gibson

J.S.

, and Dold

J.W.

, A mathematical model of the stress induced during avascular tumour growth. J. Math Biol., 40(6): 473–99, (2000).

38.

Kelley

, Coexistence in reaction-diffusion systems with mutualist or predatorprey interactions, Methods Appl. Anal., 3(1): 1–11, (1996).

39.

King

, Cancer Biology, 2nd Ed., Prentice Hall, London (2000).

40.

Kuwano

, Ushiro

, Ryuto

, Samoto

, Izumi

, Ito

, Abe

, Nakamura

, Ono

, and Kohno

, Regulation of angiogenesis by growth factors, GANN Monograph on Cancer Research, 42: 113–125, (1994).

41.

Levine

H.A.

, Sleeman

B.D.

, and Nilsen-Hamilton

, Mathematical modeling of the onset of capillary formation initiating angiogenesis, J. Math. Biol., 42(3): 195–238, (2001).

42.

Levine

H.A.

, Sleeman

B.D.

, and Nilsen-Hamilton

, A Mathematical model for the roles of pericytes and macrophages in the onset of angiogenesis: I. The role of protease inhibitors in preventing angiogenesis, Math. Biosci., 168: 77–115, (2000).

43.

Levine

H.A.

, Pamuk

, Sleeman

B.D.

, and Nilsen-Hamilton

, Mathematical Modeling of Capillary Formation and Development in Tumor Angiogenesis: Penetration into the Stroma, Bull. Math. Biol, 63(5): 801–863, (2001).

44.

Levine

H.A.

, Tucker

A.L.

, and Nilsen-Hamilton

, A mathematical model for the role of cell signaling and transduction in the initiation of angiogenesis, Growth Factors, 20: 155–176, (2002).

45.

Mizoguchi

, Ninomiya

, and Yanagida

, Diffusion-induced blowup in a nonlinear parabolic system, J. Dynam. Diff. Eq., 10(4): 619–638, (1998).

46.

Murray

J.D.

, Mathematical Biology, Springer-Verlag, Berlin, (1989).

47.

Niakan

, Steady tumor growth and the spontaneous disappearance of cancer J. Natl. Cancer Inst, 85(1): (1993), 68–9.

48.

Niakan

, Hypothesis on the biochemistry of spontaneous remissions of cancer: coupling of oxidative phosphorylation and the remission of cancer, Cancer Biother Radiopharm, 14(4): 297–8, (1999).

49.

Niakan

, A mechanism of the spontaneous remission and regression of cancer, Cancer Biotherapy Radiopharm, 13(3): 209–10, (1998).

50.

Oren

, Regulation of the p53 tumor suppressor protein, J. Biol. Chem., 274: 36031–36034, (1999).

51.

Pal

, Datta

, and Mukhopadhyay

, Central role of p53 on regulation of vascular permeability factor/vascular endothelial growth factor (VPF/VEGF) expression in mammary carcinoma, Cancer Res., 61(18): 6952–7, (2001).

52.

Painter

K.J.

, and Hillen

, Volume-filling and quorum sensing in models for chemosensitive movement, Can. J. Appl. Math., 10(4): 501–544, (2003).

53.

Patt

, Gries

, Giraldo

, Cervos-Navarro

, Martin

, Janisch

, and Brockmoller

, p53 gene mutations in human astrocytic brain tumors including pilocytic astrocytomas, Hum Pathol, 27(6): 586–9, (1996).

54.

Parikka

, Kainulainen

, Tasanen

, V ‘a'an ‘anen A., Bruckner Tuderman

, and Salo

, Alterations of collagen XVIII during transformation of oral epithelium to dysplasia and carcinoma, J. Histochem. Cytochem., 51(7): 921–929, (2003).

55.

Pavolov

M.Y.

, and Ehernberge

, Rate of translation of natrual mRNAs in an optimized in vitro system, Arch. Biochem. Biophys., 328(1): 9–16, (1996).

56.

Pescarmona

G.P.

, Scalerandi

, Delsanto

P.P.

, and Condat

C.A.

, Non linear model of cancer growth and metastasis: a limiting nutrient as a major determinant of tumor shape and diffusion, Med. Hypotheses, 53(6): 497–503,(1999).

57.

Ravi

, Mookerjee

, Bhujwalla

Z.M.

, Sutter

C.H.

, Artemov

, Zeng

, Dillehay

L.E.

, Madan

, Semenza

G.L.

, and Bedi

, Regulation of tumor angiogenesis by p53-induced degradation of hypoxia inducible factor1 α. Genes Dev, 14(1): 34–44, (2000).

58.

Retsky

, Universal law of tumor growth, J. Theor. Biol., 229(3): 289, (2004).

59.

Tomlinson

I.P.

, and Bodmer

W. F.

, Failure of programmed cell death and differentiation as causes of tumors: some simple mathematical models. Proc. Natl. Acad. Sci. USA., 92(24): 11130–4, (1995).

60.

Tsuzuki

, Fukumura

, Oosthuyse

, Koike

, Carmeliet

, and Jain

R.K.

, Vascular endothelial growth factor (VEGF) modulation by targeting hypoxia-inducible factor-1alpha–¿ hypoxia response element–¿ VEGF cascade differentially regulates vascular response and growth rate in tumors, Cancer Res., 60(22): 6248–52, (2000).

61.

Unemori

E.N.

, Ferrara

, Bauer

E.A.

, and Amento

E.P.

, Vascular endothelial growth factor induces interstitial collagenase expression in human endothelial cells, J. Cell Physiol, 153(3): 557–62, (1992).

62.

Ward

J.P.

, King

J.R.

, Mathematical modelling of avascular-tumour growth. IMA J. Math. Appl. Med. Biol, 14(1): (1997), 39–69.

63.

Weinberger

H.F.

, An example of blowup produced by equal diffusions, J. Differential Equations, 154: 225–237, (1999).

64.

Werbajh

S.E.

, Urtreger

A.J.

, Puricelli

L.I.

, de Lustig

E.S.

, Bal de Kier Joffe

, and Kornblihtt

A.R.

, Downregulation of fibronectin transcription in highly metastatic adenocarcinoma cells, FEBS Let., 440: 277–281, (1998).

65.

Zagzag

, Zhong

, Scalzitti

J.M.

, Laughner

, Simons

J.W.

, and Semenza

G.L.

, Expression of hypoxia-inducible factor 1a in brain tumors: association with angiogenesis, invasion, and progression, Cancer, 88(11): 260–618, (2000).

66.

Schirle

, Keilholz

, Weber

, Gouttefangeas

, Dumrese

, Becker

H.D.

, Stevanovic

, and Rammensee

H.G.

, Identification of tumor-associated MHC class I ligands by a novel T cell-independent approach, Eur. J. Immunol, 30: 2216–25, (2000).

67.

Thoumine

, Kocian

, Kottelat

, and Meister

J.J.

, Short-term binding of fibroblasts to fibronectin: optical tweezers experiments and probabilistic analysis, Eur. Biophys. J. 29: 398–408, (2000).

68.

Wall

F.T.

, Chemical Thermodynamics, 2nd Ed., Freeman, San Francisco and London (1958).

69.

Wodarz

, Krakauer

D.C.

, Genetic instability and the evolution of angiogenic tumor cell lines Oncol. Rep., 8(6): 1195–201, (2001).

70.

Yakovlev

, Boucher

, DiSario

, Modeling insight into spontaneous regression of tumors, Math Biosci., 155(1): 45–60, (1999).

A Mathematical Model for the Onset of Avascular Tumor Growth in Response to the Loss of P53 Function

Abstract

Introduction

Notation and Kinetics

Remark 1

Cell Movement

Mass Balance Equations

Boundary and Initial conditions

Dynamics of the Modified Lotka-Volterra System

Steady State Solutions

Initial Conditions and Normalizing Constants

Remark 2

The boundary value Yb

Evolution to Equilibrium

Results and Conclusions

Basic Simulations

Threshold and Critical Values

Footnotes

Biochemistry

Appendix B.

References

The boundary value Y_b