Fifth Mississippi State Conference on Differential Equations and Computational Simulations, Electronic Journal of Differential Equations, Conference 10, 2003, pp. 33–53. ISSN: 1072-6691. http://ejde.math.swt.edu or http://ejde.math.unt.edu ftp ejde.math.swt.edu (login: ftp)

MULTISTAGE EVOLUTIONARY MODEL FOR CARCINOGENESIS MUTATIONS

REZA AHANGAR & XIAO-BIAO LIN

Abstract. We developed a mathematical model for carcinogenesis mutations based on the reaction diffusion, logistic behavior, and interactions between normal, benign, and premalignant mutant cells. We adopted a deterministic view of the multistage evolution of the mutant cell to a tumor with a fast growth rate, and its progress to a malignant stage. In a simple case of this model, the interaction between normal and tumor cells with one or two stages of mutations was analyzed. The stability of the dynamical system and the travelling wave solutions for different stages of evolution of mutant cells were investigated. We observed the effect of variation and natural selection in shaping the carcinogenesis malignant mutation.

1. Formulations of Carcinogenesis Mutations

Introduction. In 1971, Alfred Knudson proposed a theory that united the two forms of retinoblastoma under genetic mutations. He explained that the two mutations would occur, one after another either during embryonic development or shortly after birth, in one of the cells of the retina. This elementary stage of mutation could be inherited. Consequently, a rare somatic mutation is required to trigger the explosive outgrowth of a tumor. The pair of genes inside the thirteenth human chromosome are deactivated during this mutation.

Cell Kinetic Multistage (CKM) cancer risk studied by Bogen (1989) is based on the assumption that cell proliferation follows the exponential growth and geometric model. In this approach, biological evidence indicates that precancerous cells may typically proliferate geometrically. A computer simulation of cell growth governed by stochastic processes, developed by Conolly and Kimbell, assumes that the normal cell is growing exponentially (Kimbell 1993). A stochastic process model for one, two, and three-stage transformation toward malignancy was developed for embryonic and adult mice using the Gompertzian pattern (Mao and Wheldon 1994) and another by Portier (2000) for two-stage.

2000 Mathematics Subject Classification. 92C30, 62P10, 92D10, 92D25.

Key words and phrases. Cancer, mutation, interaction, travelling waves, singular perturbation.

\begin{matrix}c \\ ◯ \end{matrix}2003 Southwest Texas State University.

Published February 28, 2003.

R. Ahangar would like to thank University of Central Arkansas, Conway, AR, for giving me

a visiting position, where most of this work was performed.

X.-B. Lin was partially supported by grant DMS-9973105 from National Science Foundation.

In this paper, we use the evolution principles of variation, natural selection, and reproduction to construct our model. In particular, we shall study how cell interaction affects and causes changes in birth, death, and mutation rates. When selection pressure from the surrounding environment is harsh, cells go through a competition stage for nutrients, space and other resources causing them to alter genetic programming in order to survive under new environmental conditions.

Under environmental pressure, the new DNA program allows cells to exploit biomasses, needed for their growth, from other resources in the body. This internal DNA program change for the cell’s survival is called adaptation. The birth of a new cell with a reprogramming of genetic makeup for adaptation is called a mutation.

In the following steps, we will demonstrate factors and conditions which will affect both normal and tumor cell growth rates either directly or indirectly.

1. Oxygen and nutrients through blood vessels: Folkman and colleagues demonstrated that solid tumors establish their own blood supplies by encouraging the growth of new blood vessels into the tumor tissue to grow beyond the small size of approximately one cubic millimeter (Folkman, 1971, Gimbron et al 1972, and Golrderg 1997). In the absence of a blood supply, the diffusion of oxygen and nutrients across numerous cell layers limits the tumor size.

2. Body’s Immune System: Owen and Sherratt studied the interactions between macrophage, tumor cells, and biochemical regulators (Owen 1997). They showed that tumors contain a high proportion of macrophages, a type of white blood cell which can have a variety of effects on the tumor, leading to a delicate balance between growth promotion and inhibition and the ability of macrophage to kill mutant cells. Macrophages also are able to lyse tumor cells over normal cells.

3. Programmed Death or Cell’s Birth and death process: An article in the “New York Academy of Science” also reviewed the idea of “programmed death” in which the cell uses signals from neighboring cells to either commit suicide or to stay alive (Raff and Wish 1996). “A cell on the verge of becoming cancerous is surrounded by normal cells that undergo apoptosis when damaged. These dying cells leave some space into which the mutated cell can grow by inducing more healthy cells to kill themselves” (Leffel and Brash 1996).

Normal cells are programmed to divide under certain conditions and do not live forever because they are programmed to die. The birth and death of cells is one of the most fascinating features of DNA programming.

The gene’s program can sometimes increase the growth of normal cells (oncogene) and sometimes limit cell growth (tumor suppressor gene). Proto-Oncogenes operate like accelerated pedals in cell proliferation, and tumor suppressor genes work like brakes (Weinberg 1998). Researchers believe that the tumor - suppressor gene or P53 protein normally stops a DNA-damaged cell from reproducing until it has had time to make repairs. Cells that become irreparably damaged rely on their death program for the greater good of the organism. If repairs are made, then P53 allows the cell cycle to continue. But if the damages are too serious to be patched, P53 activates other genes that cause the cell to self-destruct. For example, increasing levels of P53 protein in the cell causes an increase in the process of mutant cell suicide and ultimately a reduction of cancer cells. This cellular proofreading serves to erase genetic mistakes. Undetected genetic error leads the mutation process. The article presented by (Byrne 2001) demonstrates a microscopic model about tumor control by protein P53.

4. Evolution: Evolution is another important factor in cell proliferation and mutation. It is based on the principles of variation and natural selection. Variation involves changing the conditions of the environment and the cell’s internal forces. These forces are imposed on cells to alter DNA programs to adopt new conditions. The mutant cell may survive by extracting nutrients from the surrounding environment. Thus the interactions of mutant cells with their environment play an important role for survival (Davis 2000).

5. Interaction through Signals, Biochemical Regulators and Enzymes: Tumor cells and tumor-associated macrophages both release factors that can affect each other. An introduction to interactions between a cell and its neighboring cells by receiving signals was discussed in the July 1996 issue of “Scientific American”.

Mutant cells with a proliferative advantage over normal tissue cells produce a generic chemical which regulates macrophage proliferation, influx, activation, and complex formation (Michelson and Leith 1997, Dopazo et al 2001, and Witting 2001).

6. Interaction between a cell and its environment based on the principle of struggle for existence: We will study a macroscopic model of normal and tumor cells in this paper, but the microscopic factors like signal functions which affect the cell division, birth, death, and mutations will not be considered in our model.

In this step, it is assumed that surrounding environmental conditions and the genetic program are in favor of the mutant cell. The genetic program in the mutant cell makes the cell capable of fast exploitation of the space and biomass because of its need to survive and proliferate. The survival of mutant cells with their fast division rate is an indication of proper adaptation, exploitation, and availability of nutrients and space. The capability for a fast-growing rate in the tissue causes changes in the densities of normal and tumor cells at space-time (t, x). Through this conflict and competition in using resources and space, there may be a reduction in normal cell growth rate and an increase in tumor cell density rate.

In the following formulation, we will consider the macroscopic evolution of mutant cells through the stages to the development of a premalignant tumor. Consideration of all changes in DNA and internal forces on the microscopic level that affect cell birth, death, and mutations through the signal functions are beyond the scope of this paper. The effect of migration of mutant cells on tumor growth is investigated by (Pettet et al 2001).

Modeling and Formulation. Assume that cell density depends on the following hypotheses H1: Multistage Mutations: Mutations in DNA cause cancer tumors to evolve through stages with the normal stage as the initial stage and malignancy as the final stage. H2: Diffusion: At every stage of mutation, the rate of density changes by diffusion and the rate of diffusion at every stage remains constant. H3: Cell Proliferation process: Both normal and mutant cells have a limited resource environment. Thus, both behave alike logistically within certain parameters. H4: Interactions Between Stages of Mutant Cells: In all stages, cell density will be affected by quadratic interactions with cells of the current, previous, and next stage. Assume that the effect of interactions within other stages are negligible.

Let us consider Y_{i} as a density of mutant cells of the i-th stage at position (t, x) where i = 0, 1, 2, 3, …, n. We accept all hypotheses, H1 through H4, to develop a multistage model for mutant cell densities Y_{i} at stage i = 0, 1, 2, …, n, where Y_{0} will be the initial stage, Y_{i}(i = 1, 2, 3, …, n − 1) the density of intermediate stages, and Y_{n} represents the density of the final stage. The system of equations for the density function is

\frac{∂Y_{i}}{∂t} = D_{i}ΔY_{i} + a_{i}Y_{i}(1 − \frac{Y_{i}}{K_{i}}) + η_{i}Y_{i}Y_{i−1} − μ_{i+1}Y_{i}Y_{i+1},  (1.1)

for i = 1, 2, …, n−1. In this formulation the constant real numbers D_{i} are diffusion factors, a_{i} and K_{i} are logistic parameters, and η_{i} and μ_{i} are interaction coefficients for the i-th stage. In every stage i, the parameters a_{i} represent the growth rate with unlimited resources. The factor η_{i} is the growth advantage of stage i from the previous stage i − 1, and the parameter μ_{i} is the growth advantage of the next stage i + 1 from the present stage i.

Latest Stage. Note that system (1.1) does not include the normal stage, i = 0, and the final stage, i = n. The positive integer n is the number of mutations leading to the latest stage and is dictated by many factors, including nature, the body, and the type of cancer. The current stage which has not reached malignancy is called the latest stage. Suppose that the mutation is in the latest stage of development 0 ≤ i < n. Obviously the factor μ_{i+1} = 0. To have the next mutation of stage i+1, the cell should go through the variation process (interactions with changes in the environment, genetic alteration, or both) and natural selection.

The equation for the latest stage of mutation has different forms depending on whether it has growth advantage in favorable conditions or growth disadvantage in unfavorable conditions.

1) The Latest Stage of Mutation in a Favorable Environment: In this case, first stage mutant cells act like normal cells and can take advantage of the availability of resources. The equation for the latest stage is in the following form

\frac{∂Y_{i}}{∂t} = D_{i}ΔY_{i} + a_{i}Y_{i}(1 − \frac{Y_{i}}{K_{i}}) + η_{i}Y_{i}Y_{i−1},  (1.2)

for 0 < i < n.

2) The Latest Mutation in a Competition Environment: The density of cells in the latest stage model described in (1.2) does not continue to grow forever. Before mutant cells consume all of the available resources, they will go through a competition process. The body’s immune system causes pressure against the growth of mutant cells and reduces their density.

Mathematically we can describe this latest stage by a competition model:

\frac{∂Y_{i}}{∂t} = D_{i}ΔY_{i} + a_{i}Y_{i}(1 − \frac{Y_{i}}{K_{i}}) − η_{i}Y_{i}Y_{i−1}.  (1.3)

Cells in this stage are using the environmental resources but losing some to other mutant cells causing the latest mutant cells to struggle for survival.

3) The Latest Stage in a Unfavorable Environment: Conditions in (1.3) may change and the mutant cells, which are surrounded by all previous stages of mutant cells, do not have any resources other than the surrounding cells. This premalignant stage occurs when mutant cells are under environmental pressure. In mathematical language, the mutant cells are interacting with normal cells like a predator-prey model and are using the resources of the previous stage at a certain rate η_{i}. Thus, the equation for the latest stage becomes

\frac{∂Y_{i}}{∂t} = D_{i}ΔY_{i} − bY_{i} + η_{i}Y_{i}Y_{i−1}.  (1.4)

The selection pressure at this stage may cause the cells to metastasize for survival and reach the blood vessels.

4) Latest Stage is the Final Stage: The final stage of mutation occurs when cancer cells become malignant and metastasize. The precise biological data needed to specify the number of stages for different kinds of cancers is not known. In melanoma, for example, it takes four mutations for cells to evolve to malignant mutations, and when they reach the blood vessels, it is said to metastasize. We assume that for the malignant mutation in the final stage, the carrying capacity is unlimited; that is K_{n} → ∞, thus equation (1.1) becomes

\frac{∂Y_{n}}{∂t} = D_{n}ΔY_{n} + a_{n}Y_{n} + η_{n}Y_{n}Y_{n−1}.  (1.5)

When the evolution of the mutations of normal cells reaches this stage, the system will blow up, but it is not of our interest of study.

For i = 0, equation (1.1) turns out to be the equation for the initial or normal stage and will be denoted by Y_{0} = u and the final stage Y_{n} = M. We also can postulate that the malignant tumor is growing exponentially rather than logistically. When the i^{th} stage of evolution of the normal tumor is not the initial or final stage, we call it the intermediate stage or Benign. The full blown developed malignant mutation can be described in the following system of the n-stage model, i.e. initial, benign, and malignant,

\begin{matrix}\frac{∂u}{∂t} = D_{0}Δu + a_{0}u(1 − \frac{u}{K_{0}} ) − μ_{1}uY_{1}, \\ \frac{∂Y_{i}}{∂t} = D_{i}ΔY_{i} + a_{i}Y_{i}(1 − \frac{Y_{i}}{K_{i}}) + η_{i}Y_{i}Y_{i−1} − μ_{i+1}Y_{i}Y_{i+1}, \\ \frac{∂M}{∂t} = D_{n}ΔM + a_{n}M + η_{n}MY_{n−1}, \end{matrix}  (1.6)

where i = 1, 2, 3, …, n − 1. Coefficients η_{i} and μ_{i}, for i = 0, 1, …, n, are constants and are called interaction coefficients.

Initial and Boundary Conditions. The initial conditions for systems (1.1)–(1.5) are I.C:

u(x, 0) = u_{0}(x), Y_{i}(x, 0) = Y_{i0}(x), M(x, 0) = M_{0}(x), in I × Ω

The boundary conditions of system (1.8) will be B.C:

\frac{∂u}{∂n} = 0, \frac{∂Y_{i}}{∂n} = 0, (i = 1, 2, …, n − 1), t > 0, x ∈ ∂Ω.

The boundary conditions are to be interpreted as no flux conditions. This means that the mutation is not in the metastasized stage and there is no migration of cells across ∂Ω. For all of these stages, Ω is considered the only habitat for cancer cells.

Formulation of Two-Stage Mutation. For simplicity we consider a two-stage mutation model: benign and the latest stage of premalignancy. When all conditions are in favor of the mutant cell, the latest stage may lead to malignant mutation. However, we are interested in studying the latest stage under selection pressure.

To demonstrate the mathematical form of the selection pressure, we adopt equation (1.4) for the latest stage. Thus the two-stage model with premalignancy as latest stage will be in the following form

\begin{matrix}\frac{∂u}{∂t} = D_{0}Δu + a_{0}u(1 − \frac{u}{K_{0}} ) − μ_{1}uv, \\ \frac{∂v}{∂t} = D_{1}Δv + a_{1}v(1 − \frac{v}{K_{1}} ) + η_{1}uv − μ_{2}vw, \\ \frac{∂w}{∂t} = D_{2}Δw − a_{2}w + η_{2}vw. \end{matrix}  (1.7)

The initial conditions are u(0) = u_{0}, v(0) = v_{0},and w(0) = w_{0} with no flux on the boundary.

2. One-Stage Mutation Interacting System

Epidemiologists succeeded in converting normal cells to cancer cells by introducing single oncogenes into them. A single oncogene was enough to create a malignant cancer cell in one hit (Weinberg 1996). We accept the one-stage model for the purpose of its simplicity of evolution in developing the cancer cell. This type of transformation of normal cell to a full blown cancer cell is possible experimentally in the laboratory by using chemical agents (Weinberg 1998). However, cancer formation is a complex process involving a long sequence of steps, rather than a simple one-hit event that converts a fully normal cell into a highly malignant one in a single step. When we eliminate all intermediate stages, the mathematical form of (1.7) include the initial, the intermediate, and the latest stage of evolution toward cancer cells.

We will present a single hit mutation for the one-stage model. Michelson and Leith (1997) studied this type of interaction between tumor cells with no diffusion factors of angiogenesis in tumor growth control. By certain assumptions which are dependent on the nature of the tumor on one side and the genetic program with environmental conditions on the other side, this model will be reduced to a particular case of tumor growth for which enormous work has been done during the past few decades in mathematical biology. For example, in the one-stage model when there is no diffusion factor, systems (1.3) and (1.4) will be reduced to Lotka-Volterra equations or to the logistic case where it will be the Pearl-Verhulst model.

One-Stage Carcinogenesis Mutation in an Unfavorable Environment. Our mathematical form should be able to explain the principles of evolution: variation, natural selection, and survival of the fittest. In the following model, we assume that mutant cells go through the selection pressure in an unfavorable environment. The more aggressive mutant cells are able to exploit the environment and the resources of cells of previous stages and have a better chance to survive. Since we are eliminating the intermediate stages, the environment and normal cells provide resources for the latest stage.

Assume that the initial stage cells are the only resource for the growth of the mutant cells. Then the quadratic interaction behaves as a predator-prey model. In this case, the body’s immune system is very strong and is harshly active against the growth of the tumor cell, that is, the density of the mutant cell in the absence of nutrients is negatively proportional to its size.

When the environment is not in favor of the mutant cell, then according to one of the principles of evolution, the mutant cell will face selection pressure. If it cannot adapt to a new condition, it will not survive. The mutant cells that can use the resources of normal cells have a better chance of survival. With these conditions imposed on the mutant cells, the model will be

\begin{matrix}\frac{∂u}{∂t} = D_{1}Δu + au(1 − \frac{u}{K_{1}} ) − r_{1}uv, \\ \frac{∂v}{∂t} = D_{2}Δv − ev + r_{2}uv, \end{matrix}  (2.1)

where e, r_{1}, and r_{2} are positive constant real numbers. Dunbar (1983) studied this model. See also Smitalova and Sujan (1991). By rescaling (2.1) in one dimension, that is,

U = \frac{u}{K_{1}} , V = \frac{r_{1}}{a} , t′ = at, x′ = x(D_{2}/a)^{−1/2} = (\frac{a}{D_{2}} )^{1/2}x

D = \frac{D_{1}}{D_{2}} , γ = \frac{r_{2}K_{1}}{a} , δ = \frac{e}{r_{2}K_{1}} ,

we obtain

\begin{matrix}U_{t} = DU_{xx} + U(1 − U − V), \\ V_{t} = V_{xx} + γ(U − δ)V. \end{matrix}  (2.2)

The necessary condition for the survival of mutant cells in this stage is 0 < δ = \frac{e}{r_{2}K_{1}} < 1. This means that the tumor growth factor e cannot be very large, but its interaction capability against normal cells r_{2} and the carrying capacity of normal cells will help the tumor to survive. The travelling wave exists in (2.2) when

D = \frac{D_{1}}{D_{2}} is very small and negligible. Thus,

U(x, t) = w(x + ct), V (x, t) = z(x + ct)

where the wave with a single variable s = x + ct is moving in a stationary system with the speed c. Substituting in (2.2) yields

\begin{matrix}cw′ = Dw′′ + w(1 − w − z), \\ cz′ = z′′ + γz(w − δ). \end{matrix}  (2.3)

The second order nonlinear system can be substituted by a first order system of ODE with respect to the variable s

\begin{matrix}w′ = (1/c)w(1 − w − z), \\ z′ = y, \\ y′ = cy − γz(w − δ). \end{matrix}  (2.4)

For system (2.4) there are two pairs of critical points whose connections will be the solution to the system. The non-negative solutions satisfying the condition

\begin{matrix}w(−∞) = 1, w(∞) = δ, \\ z(−∞) = 0, z(∞) = 1 − δ \end{matrix}  (2.5)

are called type I solutions. Non-negative solutions of (2.5) satisfying

\begin{matrix}w(−∞) = 0, w(∞) = δ, \\ z(−∞) = 0, z(∞) = 1 − δ \end{matrix}  (2.6)

are called type II solutions. The following conclusion for D = 0 can be shown (Dunbar 1983).

Theorem 2.1. The travelling wave front solutions (w(s), z(s)) of system (2.2) satisfying condition (2.6) (type II) exist if 0 < c < \sqrt{4γ(1 − δ)}. The travelling wave front solution to (2.2) satisfying (2.5) exists when c ≥ \sqrt{4γ(1 − δ)}.

The reaction term of system (2.2) has three equilibrium points. The equilibrium (0, 0) represents meaningless extinction of normal and mutant cells. The point (1,0) represents the extinction of mutant cells which is desirable for every cancer patient. But for the sake of studying the evolution, behavior, and survival of mutant cells, this critical point is not interesting. Both of these points are unstable.

The equilibrium (δ, 1 − δ) represents the stable coexistence of both normal and mutant cells. At this stage in the development of cell mutation which has reached the stage of metastasizing, the migration factor D_{1} is very small relative to D_{2}.

3. The Latest Stage Mutation in Unfavorable Conditions

In this section, we will study the case where mutations occur for the second time in unfavorable conditions.

The premalignant stage is the latest stage to reach the malignant stage if the mutant cells survive selection pressure. Environmental conditions are the types of interactions which may or may not be favorable for the latest stage of mutant cells. When we treat the one-stage tumor development, if changes in the environment are not favorable for the life of the latest stage mutant cell, then it causes a force of selection pressure on the mutant cells. The fittest are those cells which can develop the ability to change for adaptation. These changes in the cell’s DNA will lead to a new mutation.

In this paper, we consider the following five factors, C1 through C5, to develop, analyze, and study the premalignant mutations. C1: Mutant cells are using the surrounding cells of previous stages as survival resources and have not yet reached the final malignant stage. We are considering all interactions of the i^{th} stage with the previous stage (i-1) and the latest stage i+1. The effect of all other interactions on the cell density will be negligible. C2: Analysis of Multistage with Premalignant Stage: We want to study the evolution of a two-stage mutation model under several different environmental conditions. C3: Evolution: We are considering the selection pressure in the latest stage, meaning that the surrounding cells are the only resources for the latest-stage mutant cells. C4: Stability: We study the instability of the latest-stage mutant cells which may cause another mutation. C5: We will measure the migration factors to see the abnormality of the mutant cells in transition to a metastasized stage mutation.

For simplicity, we drop all subscripts and use u, v,and w to represent the densities of the three stages Y_{i−1}, Y_{i}, and Y_{i+1} respectively. We use M for the density of the Figure 1. Simulation of the travelling wave solution U(t) over the parameter 0 ≤ D ≤ 3 measuring the migration factor.

malignant stage. Based on different conditions which govern the latest stage of mutant cells, the mathematical formulations will be different.

For a three-stage mutation model which could lead to malignant cells, the evolution of the stages is of the following form:

Normal Cell ⟶ Benign ⟶ Premalignant ⟶ Malignant .

For example, in skin cancer, normal cells will undergo the following stages: initial stage (no mutation, but cells are susceptible to mutate), rapid growth (divided rate mutation), invasion mutation, angiogenesis, and metastasis (malignant stage). Remark: The stochastic simulation model when n = 4 has been established by Sherman and Portier (1996).

Numerical Solutions. The following algorithm will provide the numerical solutions to the system of singularly perturbed equations with two-stage mutations. All parameters related to growth, birth, death, mutation, and interaction rates are from Conolly and Kimbell (1994). We are using the Winpp computational tools to simulate the travelling wave solution to the system of PDEs for 0 ≤ x ≤ 40.

Simulation of the travelling wave solution U(t) over the parameter D is depicted in Fig.1 and V (t) over the parameter δ is depicted in Fig.2. In Figure 1, the parameter D represents the migration capability of mutant cell with respect to normal cells. When D is close to one, the mutant cell is metastasized. The smallness of D may lead to the high capability of the travelling of the mutant cells for adaptation and survival. The simulation of V (t), in Fig.2 over the parameter δ, can be interpreted as the harshness of the environment against mutant cells which may lead to their adaptation and survival. This model demonstrates the evolution of the mutant cell in a simple form of macroscopic one-stage mutation without considering the cell’s genetic changes and its relation to the parameters D and δ.

Figure 2. Simulation of the travelling wave solution V(t) over the parameter 0 ≤ δ ≤ 4. This parameter measures how much the environment favors mutant cells.

In the next section, we will show that this coexistence may last years before it is disturbed by another premalignant mutation.

# Program Using Winpp or Xtc
# discretization of the singular perturbation system of
# reaction diffusion equations
# USING TWO-STAGE CARCINOGENESIS MUTATIONS
# with no flux boundary conditions
# Notation: k1 is used for epsilon and c is replacement for the
# parameter n in the system of PDE.
#
par h=1,b=.7,c=.2,m=.3,e=.4,D=2.5,k1=.0008
u0'=u0*(1-u0-m*v0)/k1+(k1)*(u1-u0)/(h^2)
v0'=v0*(b-m*u0-b*v0-c*w0)+D*(v1-v0)/h^2
w0'=-w0*(1-e*w0)+(w1-w0)/(h^2)
u[1..39]'=u[j]*(1-u[j]-m*v[j])/k1 +(k1)*(u[j-1]-2*u[j]+u[j+1])/(h^2)
v[1..39]'=v[j]*(b-m*u[j]-b*v[j]-c*w[j])+D*(v[j-1]-2*v[j]+v[j+1])/(h^2)
w[1..39]'=-w[j]*(1-e*w[j])+(w[j-1]-2*w[j]+w[j+1])/(h^2)
u40'=u40*(1-u40-m*v40)/k1+(k1)*(u39-u40)/(h^2)
v40'=v40*(b-m*u40-b*v40-c*w40)+D*(v39-v40)/h^2
w40'=-w40*(1-e*w40)+(w39-w40)/(h^2)
init u0=.1,v0=.03,w0=.4
@ total=40,dt=.2,meth=qualrk
@ xhi=40,yhi=1,ylo=0,yp=5,zhi=40
done
#

For detailed information about this program is available at the web site http://www.math.pitt.edu/pub/bardware To download this program, use ftp://ftp.pitt.edu/pub/bardware

4. The Latest Premalignant Stage in Unfavorable Conditions

To study the possibility of the final malignant mutation, we study the stability or instability of the premalignant stage of mutation. Accepting hypothesis H1-H4 and conditions C1-C5, we will pursue the following for further development of the two-stage model of system (1.7) by

i) Changing variables:

U = \frac{u}{K_{0}} , V = \frac{v}{K_{1}} , W = w.  (4.1)

As a result of this change, u_{xx} = K_{0}U_{xx}, v_{xx} = K_{1}V_{xx}, w_{xx} = W_{xx}. ii) Changing the coordinate system (t, x) → (τ, ξ):

τ = a_{2}t, ξ = \sqrt{a_{2/D_{2}}}x.  (4.2)

According to this transformation (U_{t}, V_{t}, W_{t}) = a_{2}(U_{τ} , V_{τ} , W_{τ}),

(u_{xx}, v_{xx}, w_{xx}) = (\frac{a_{2}}{D_{2}} )(U_{ξξ}, V_{ξξ}, W_{ξξ}).

System (1.7) will be in the form

\begin{matrix}U_{τ} = \frac{D_{0}}{D_{2}} U_{ξξ} + \frac{a_{0}}{a_{2}} U(1 − U) − \frac{μ_{1}K_{1}}{a_{2}} UV \\ V_{τ} = \frac{D_{1}}{D_{2}} V_{ξξ} + \frac{a_{1}}{a_{2}} V (1 − V) + \frac{η_{1}K_{0}}{a_{2}} UV − \frac{μ_{2}}{a_{2}} V W, \\ W_{τ} = W_{ξξ} − W + \frac{η_{2}K_{1}}{a_{2}} V W. \end{matrix}  (4.3)

iii) Rename the Coefficients: If the constant coefficients are not unusually small or large, redefine them by

\begin{matrix}\frac{μ_{1}K_{1}}{a_{2}} = c, \frac{η_{1}K_{0}}{a_{2}} = m, \frac{μ_{2}}{a_{2}} = n, \frac{η_{2}K_{1}}{a_{2}} = e, \\ \frac{D_{0}}{D_{2}} = d_{1}, \frac{D_{1}}{D_{2}} = d_{2}, \frac{a_{0}}{a_{2}} = a, \frac{a_{1}}{a_{2}} = b. \end{matrix}  (4.4)

For convenience, label our new coordinate system (t, x) so that system (4.3) becomes

\begin{matrix}U_{t} = d_{1}U_{xx} + aU(1 − U) − cUV, \\ V_{t} = d_{2}V_{xx} + bV (1 − V) + mUV − nV W, \\ W_{t} = W_{xx} − W + eV W. \end{matrix}  (4.5)

iv) Introduce the Perturbation Parameter: To study the carcinogenesis mutation using model (4.3), we impose a few more conditions on the parameters in (4.4). The ratio \frac{a_{1}}{a_{2}} is not unusually large or small, but \frac{a_{0}}{a_{2}} is denoted by \frac{1}{ε_{1}} (where ϵ_{1} ≪ 1, implies smaller values for a_{2}). The migration rate D_{2} is not very large with respect to the migration rate of benign cells, but is large in comparison to the diffusion factor D_{0}. Thus, the second line of the parameters in (4.4) become

\frac{D_{0}}{D_{2}} = ε_{2}, \frac{D_{1}}{D_{2}} = D, \frac{a_{0}}{a_{2}} = 1/ε_{1}, \frac{a_{1}}{a_{2}} = b.  (4.6)

For simplicity but without loosing generality, we will replace ε_{i} = ε (for i = 1, 2). Accepting all of these changes in (i), system (4.3) becomes

\begin{matrix}U_{t} = εU_{xx} + \frac{1}{ε}U(1 − U) − \frac{1}{ε} \frac{μ_{1}K_{1}}{a_{0}} UV, \\ V_{t} = DV_{xx} + bV (1 − V) + mUV − nV W, \\ W_{t} = W_{xx} − W + eV W. \end{matrix}  (4.7)

Renaming the value \frac{μ_{1}K_{1}}{a_{0}} = γ, system (4.7) leads to the system

εU_{t} = ε^{2}U_{xx} + U(1 − U − γV),  (4.8-a)
\begin{matrix}V_{t} = DV_{xx} + bV (1 − V) + mUV − nV W, \\ W_{t} = W_{xx} − W + eV W. \end{matrix}  (4.8-b)

This is a system of singularly perturbed equations where U(t, x) in (4.8-a) is the “fast variable”, whereas V (t, x) and W (t, x) in (4.8-b) are the “slow variables”.

For further clarification, let us review the justifications behind our assumptions and goals in the obtaining of three singularly perturbed system of equations (4.8) from the nonlinear multistage system of (1.4).

- U, V, and W represent the densities of three stages of mutant cells (after rescaling). W is the latest stage and is considered the premalignant stage. W cells are the mutations inside the V cells environment and are under selection pressure for adaptation and survival. The negative term −W represents the evolutionary selection pressures, and the last term eV W represents the harsh interactions with surrounding cells. The parameter ε in the first equation measures, on one hand, the relative movement of \frac{D_{0}}{D_{2}} = ε (W with respect to U) and, on the other hand, the effect of the selection pressure on the reduction of the growth rate a_{2} with respect a_{0} ( \frac{a_{0}}{a_{2}} = 1/ε getting very large). We use the singular perturbation method to study the stability of the system.

In the next section, we would like to present the analysis of the singularly perturbed solution of the dynamical systems and the numerical methods to approximate the travelling wave solutions.

5. Singularly Perturbed Solution to the Dynamical System

The geometrical interpretation of the approximated solution using a dynamical system approach is the existence of travelling wave solutions for a class of parabolic PDE’s. The desired trajectory can be constructed as a transverse intersection of a pair of invariant manifolds. The actual solution (if it exists) is near the singular solution which must pass close to the slow manifolds.

To resolve problem (4.8) we would like to accept the condition that mutant cells take advantage of all space and nutrients against normal cells, that is, γ = m.

System (4.8) will be in the following form

εU_{t} = ε^{2}U_{xx} + U(1 − U − mV),  (5.1-a)

\begin{matrix}V_{t} = DV_{xx} + bV (1 − V) + mUV − nV W, \\ W_{t} = W_{xx} − W + eV W, \\ U(0) = U_{0}, V (0) = V_{0}, W (0) = W_{0}. \end{matrix}  (5.1-b)

For sufficiently small ε, the travelling wave solution of system (5.1) is an internal layer solution. Since the wave speed is dependent on the shape of initial distribution, Dunbar Steve (1983) conjectured that the stable travelling wave solutions of the system with D = 0 are the singular limit solutions of stable travelling wave solutions with D ≠ 0.

Theorem 5.1 (Existence of Traveling Wave Solution). The singularly perturbed system (5.1) has a travelling wave solution.

Proof. The transversality of the intersection shows the connections between critical points of three dimensional phase space. First, assume in system (5.1) that the parameter m = 0. The first equation which is the Fisher-Kolomogorov model

U_{t} = εU_{xx} + \frac{1}{ε}U(1 − U)  (5.2)

is decoupled from the other two equations and has a travelling wave U(s). Dunbar (1983) proved that if D = 0 or 0 < D ≤ 1, travelling wave solutions exist for the last two equations in (5.1-b)

\begin{matrix}V_{t} = DV_{xx} + bV (1 − V) − nV W, \\ W_{t} = W_{xx} − W + eV W. \end{matrix}  (5.3)

At certain parameter values the wave speeds for U(s) and (V (s), W (s)) agree. Assume that ε = 1 and m is nonzero but very small. Using the transverse intersection of unstable manifold of one equilibrium with the stable manifold of another equilibrium, we find that the travelling wave solution persists for small and nonzero m. This proves that the travelling wave solution for systems (5.1-a) and (5.1-b) exists. □

Discussions About Singular Perturbation Method. When ε → 0 and m = 1 the first equation in (5.1) becomes algebraic, and the manifold of critical points is defined by

U(1 − U − V) = 0.

The center manifolds or slow manifolds is produced by U = 0 or U = 1 − V. We can use these center manifolds to approximate the solutions. These so-called singular solutions are the union of solutions to the “slow equations” and the “fast equations”. Substituting these values of U into the other two equations, we have: I) For U = 0,

\begin{matrix}V_{t} = DV_{xx} + bV (1 − V) − nV W, \\ W_{t} = W_{xx} − W + eV W. \end{matrix}  (5.4)

II) For U = 1 − V,

\begin{matrix}V_{t} = DV_{xx} + (b + m)V (1 − V) − nV W, \\ W_{t} = W_{xx} − W + eV W. \end{matrix}  (5.5)

According to Dunbar, both of these systems have travelling front wave solutions. That is, there are two sets of solutions (0, V (s), W (s)) and (1 − V (s), V (s), W (s)) satisfying (5.4) and (5.5) respectively. There is a point x = x_{0} such that this solution for U jumps from zero to 1 − V, that is,

U(t, x) = \begin{cases}0, & for x < x_{0}, \\ 1 − V, & for x > x_{0}, \end{cases}  (5.6)

where at the point x = x_{0}, the value of Ṽ = V (x_{0}), and there is a jump for the value of U from U = 0 to Ũ = 1 − Ṽ.

To study the solution of the first equation, after rescaling

ξ = (x − x_{0})/ε, τ = t/ε,

the internal layer will satisfy the “slow variable” system

U_{τ} = U_{ξξ} + U(1 − U − mṼ).  (5.7)

This equation has a travelling wave front solution connecting U = 0 to Ũ = 1 − V (x_{0}) with the same speed of the last two equations.

6. Concluding Remarks

The following is a brief review of our results. 1- Variation and Natural Selection: The tumor formed by a one-stage mutation may approach a stable equilibrium state. By one of the evolution principles of “variation”, the environment, and consequently, the DNA algorithm will change. This equilibrium state persists until new conditions force the mutant cells to develop an adaptation procedure to survive. Some conditions are not favorable for mutant cells, but those that can adapt to the new conditions will survive. These changes are inevitable and may cause a sequence of mutations.

If the cell is in the latest stage of its mutation, and has not reached the final malignant stage under unfavorable conditions such as lack of oxygen and nutrients, it is in the premalignant stage. 2- Mathematical Representations for Stages: A deterministic and macroscopic model representing the evolution of carcinogenesis mutation using parabolic PDEs was developed. 3- Stability and instability: In order to study the chance of malignant mutation, we have studied the instability of premalignant mutations. This will lead us to understand the chance of another mutation toward the malignant stage. Particularly, our interest has been to study the evolution of premalignant cells when unfavorable conditions are imposed on the mutant cells by selection pressures. 4- One-stage with Diffusion factor in Unfavorable Conditions: If the environment is not in favor of mutant cells, the following system can represent the normal and tumor densities and their interactions

\begin{matrix}\frac{∂u}{∂t} = D_{1}Δu + au(1 − \frac{u}{K_{1}} ) − r_{1}uv, \\ \frac{∂v}{∂t} = D_{2}Δv − ev + r_{2}uv. \end{matrix}  (6.1)

When D = \frac{D_{1}}{D_{2}} is very small, we have the following statement.

Theorem 6.1. The travelling wave front solutions (w(s),z(s)) of system (2.2) satisfying condition (2.6) (type II) exists if 0 < c < \sqrt{4γ(1 − δ)}. The travelling wave front solution to (2.2) satisfying (2.5) exists when c ≥ \sqrt{4γ(1 − δ)}.

The reaction term of system (2.2) has three equilibrium points. The equilibria (0, 0) represents the meaningless extinction of normal and mutant cells. The point (1, 0) represents another uninteresting stage of the extinction of mutant cells since we want to study the evolution of mutant cells. Both of these points are unstable. The stable equilibrium (δ, 1 − δ) represents coexistence of both normal and mutant cells. At this stage, cell mutation has reached metastasis. The migration factor D_{1} is very small relative to D_{2}. This implies D = \frac{D_{1}}{D_{2}} ≅ 0. The simulation of the travelling wave solution to system (2.2) has been demonstrated in Fig.1 and Fig.2. 5- Two-Stage Model: For the two-stage model, the normal cell will undergo two mutations. In the second stage, the mutant cells are not in a favorable environment. We used singular perturbation to prove the existence of a stable travelling wave front. Consequently, the survival of the mutant cell under selection pressure may lead to another mutation. The following system

\begin{matrix}εU_{t} = ε^{2}U_{xx} + U(1 − U − mV), \\ V_{t} = DV_{xx} + bV (1 − V) + mUV − nV W, \\ W_{t} = W_{xx} − W + eV W. \end{matrix}  (6.2)

is the singular perturbation form. We used the idea of transverse intersection of unstable manifold of one equilibrium with the stable manifold of another equilibrium preserving the travelling wave solution. This idea can be extended to a higher dimension or to a model with more than two-stage mutations. 6- Final Malignant Stage: There will not be any coexistent state when the malignant cells metastasize which means that there will not be any stable equilibrium point to the dynamical system (5.2 ) when a single malignant cell is an aggressive survivor. This is the time when the malignant cells aggressively metastasize and use the relatively unlimited resources of the body.

The numerical method developed in the previous section can be applied to system (5.1) using Winpp. All figures in this paper were produced by this algorithm and are approximate solutions to the singularly perturbed system (5.1). Simulations of the travelling wave solutions (U(s), V (s), W (s)) over the variations of indicated parameters and initial data are demonstrated later in this paper (see Figures 3–9). In these figures, ϵ measures the mobility of a premalignant cell with respect to the normal cell. ϵ also measures the selection pressure on the premalignant mutant cell. By making ϵ small, we increase the moving capability and selection pressures on the premalignant cell. Simulation over the variation of ϵ can be observed in Fig. 3 and Fig. 6 for producing a travelling wave solution of U when x ranges from 0 to 40. Simulation over the variations of two parameters D and ϵ and generation of the solution of V is depicted in Fig. 5. The parameter e in system (6.2) measures the aggressive nature of the premalignant cell in its struggle for existence. Simulation over this parameter demonstrates the survival of the premalignant mutation in Fig. 8 and Fig. 9. The simulation of the solution of W when x = 0 to x = 12 over the parameters e and D shows the travelling front and coexistence of the singularly perturbed system of (6.2) (see Figures 7-9).

Figure 3. Simulation of the evolution of normal cell density over the variation of epsilon 0 ≤ k1 = ε ≤ 1.

Future Development. In this paper we investigated the possibilities of another mutation for the premalignant stage through instability caused by evolutionary selection pressure using the travelling wave solution of the singular perturbation system. The following problems are open for further development:

1. Although extensive work has been done on multistage stochastic modelling, the deterministic system of (1.4) can be applied for a particular cancer with more than two stages. Notice that the paper mentioned above by Portier (2000) is a probabilistic development of a two-stage carcinogenesis model.

2. The quadratic interaction used in our model may be consistent with certain kinds of cancer. Testing the behavior of a variety of interaction functions for different kinds of cancer will help to study the behavior of tumor growth.

3. In this macroscopic model, we assumed all constant parameters are independent from all internal actions and reactions including the effect of continuous changes of environment on DNA leading to the mutation. A microscopic formulation which considers the cell’s timing clock for birth, death, and mutation is required.

4. Developing the stochastic model by imposing a random perturbation on the deterministic model caused by random change in the environment or gene will be very interesting.

Figure 4. Simulation of the travelling wave solution V (t) over the parameter 0 ≤ δ ≤ 4. This parameter measures how much the environment favors mutant cells.

Figure 5. Simulation of the singularly perturbed solution V (t) to equations (5.2) and (5.3) with the parameters h = 1, c = 0.2, b = 0.7, e = 0.4, m = 0.3, 0 ≤ u_{0} ≤ 1 over the parameters 0 ≤ D ≤ 3 and epsilon 0 ≤ k1 ≤ 1.

Figure 6. Simulation of the evolution of mutant cell density U(t) over parameter epsilon 0 ≤ k1 ≤ 1.

Figure 7. Simulation of the evolution of the singularly perturbed solution W (t), premalignant cell density over the parameters 1 ≤ e ≤ 4 and 0 ≤ D ≤ 3.

Figure 8. Simulation of the solution W (t) over the parameters 0 ≤ e ≤ 4 and 0 ≤ D ≤ 3.

Figure 9. Simulation of the singularly perturbed solution W (t) : 0 ≤ x ≤ 12 over the parameters 1 ≤ e ≤ 4 and 0≤ k1 ≤ 1.

References

[1] Bogen, “Cell Proliferation Kinetics and Multistage Cancer Risk Model”, Journal of the National Cancer Institute, Vol. 81, No. 4, Feb. 15, 1989.

[2] Byrne, H. M.; Chaplain ,M. A. J.; “Mathematical Models for Tumor Angiogenesis: Numerical Simulations and Nonlinear Wave Solutions”, Bulletin of Mathematical Biology, Vol. 57, No. 3, pp. 461-486, 1995.

[3] Byrne, H. M. and Gammack, D.; “ Estimating the Selective Advantage of Mutant P53 Tumor Cells to Repeated Rounds of Hypoxia”, Bulletin of Mathematical Biology (2001) 63, 135-166.

[4] Davis, Brian K.; “Darwinian Aspects of Molecular Evolution at Sublinear Propagation Rates”, Bulletin of Mathematical Biology (2000) 62, 675-694.

[5] Depazo, H.; Gordon, M. B.; Perazzo, R. and Risau-Gusman, S., “A Model for the Interaction of Learning and Evolution”, Bulletin of Mathematical Biology, (2001) 63, 117-134.

[6] Dunbar, Steve R.; “Traveling Wave Solutions of Diffusive Lotka-Volterra Equations”, J. Math. Biology 1983, 17:11-32.

[7] Dunbar, Steve R.; “Traveling Waves in Diffusive Predator-Prey Equations: Periodic Orbits and Point-To-Periodic Heteroclinic Orbits”. SIAM J. Appl. Math, vol. 46, No. 6, December 1996.

[8] Maitland, A. Edey and Donald C. Johnson; “Blueprints Solving the Mystery of Evolution”, Penguin Books, 1989.

[9] Fitzgibbon, W. E.; “Partial Differential Equations and Dynamical Systems”, Research Notes in Mathematics 101, Pitman Advanced Publishing Program 1984.

[10] Frank-Kamenetskii, Maxim D.; “Unraveling DNA” Translated from Russian by Lev Liapin, Published by VCH, 1993.

[11] Goldberg, I. D. and Rosen, E. M.; “Regulation of Angiogenesis” Birkhauser Verlag, 1997.

[12] Ken-On, Y., “Fisher-Type Property of Traveling Waves for Competition-Diffusion Equations”, Research Report CDSNS94-159, GA-Tech, 1994.

[13] Conolly, Rory B. and Kimbell, Julia S.; “Computer Simulation of Cell Growth Governed by Stochastic Processes: Application to Clonal Growth Cancer Models”, Toxicology and Applied Pharmacology 124, 284-295 (1994).

[14] Kirschner, Denise; “Using Mathematics to Understand HIV Immune Dynamics”, AMS Notices, Feb, 1996, Vol 43, No.2.

[15] Leffell, David J. and Brash, Douglas E; “Sunlight and Skin Cancer”, Scientific America, July 1996.

[16] Mao Jian-Hua and Wheldon Tom E., “A Stochastic Model for Multistage Tumorigenesis in Developing Adult Mice”, Mathematical Biosciences, vol 129, Number 2, pp95-110, October 1996.

[17] Moolgavkar Suresh H., “Carcinogenesis Modeling: ¿From Molecular Biology to Epidemiology”, Ann. Rev. Public Health 1986, 7.151-69.

[18] Michelson, Seth and John T. Leith; “Positive Feedback and Angiogenesis in Tumor Growth Control”, Elsevier Science Inc. Bulletin of Mathematical Biology, Society for Mathematical Biology, Vol. 59, No.2, pp. 233-254, 1997.

[19] Owen, Markus R. and Sherratt, Jonathan A.; “Pattern formation and Spatiotemporal Irregularity in a Model for Macrophage-Tumour Interactions”, Jan. 1997, Nonlinear System Laboratory, Mathematics Institute, University of Warwick, Coventry CV4 7AL, U.K.

[20] Pettet, G. J.; Please C. P.; Tindall, and McElwain D.L.S.; “ The Migration of Cells in Multicell Tumor Spheroids”, Bulletin of Mathematical Biology (2001) 63, 231-257.

[21] Sherman, Claire D. and Portier Christopher J.; “Calculation of Cumulative Distribution Function of the Time to a Small Observable Tumor”, Bulletin of Mathematical Biology (2000) 62, 229-240

[22] Portier, C. J., Smith, M. V.; “Incorporating Observability Thresholds of Tumors into the Two-Stage Carcinogenesis Model”, Mathematical Biosciences, Elsevier, 163, (2000) 75-89.

[23] Raff, Martin C.; “Death Wish” The Sciences Published by New York Academy of Sciences. August 1996.

[24] John L. Spouge, Richard Shrager, and D. S. Dimitrov, “HIV-1 Infection Kinetics in Tissue Cultures”, Mathematical Biosciences, 138:1-22, Elsevier, March 1996.

[25] Schinazi, Rinaldo; “On an interaction particle modeling an epidemic”, by Mathematical Biology, Spring Verlag 1996, 34: 915-925.

[26] Sherman, Claire D. and Portier, Christopher J.; “Stochastic Simulation of a Multistage Model of Carcinogenesis”, Mathematical Biosciences, volume 134, Number 1, May 1996.

[27] Smitalova, Kristina and Sujan Stefan; “A Mathematical Treatment of Dynamical Models in Biological Science,” By Translated from Czechoslovakia, by B. Sleeman, Ellis Harwood 1991.

[28] Wiarda, D.,and Travis, C.C.; “Determinibality of Model Parameters in a Two-Stage Deterministic Cancer Model” Mathematical Biology, 1997.

[29] Weinberg, Robert; “How Cancer Arises. …”, Scientific American, September 1996.

[30] Weinberg, Robert; “One Renegade Cell, How Cancer Begins”, The Science Masters Series Publisher, 1998.

[31] Witting, Lars; “Population Cycles Caused by Selection by Dependent Competitive Interactions”, Bulletin of Mathematical Biology, (2000), 62, 1109-1136.

[32] Wu, Jianhong; “Theory and Applications of Partial Functional Differential Equations,” Applied Mathematical Sciences 119, Springer 1996.

[33] Yanagida, E. and Ei, S. I.; “Dynamics in Competition-Diffusion Systems”, SIAM J. Appl. Math., vol. 54, No. 5, pp. 1355-1373, October 1994.

Reza Ahangar Math Department, Kansas Wesleyan University, Salina, KS 67401, USA

   E-mail address: rahangar@kwu.edu or rahangar@alltel.net

Xiao-Biao Lin Math Department, North Carolina State University, Raleigh, NC 27695-8205, USA

   E-mail address: xblin@math.ncsu.edu