\documentclass[reqno]{amsart}
\usepackage{graphicx}     %  for including eps figures.

\AtBeginDocument{{\noindent\small
Fifth Mississippi State Conference on Differential Equations and  
Computational Simulations,
{\em Electronic Journal of Differential Equations},
Conference 10, 2003, pp. 33--53.\newline
ISSN: 1072-6691. http://ejde.math.swt.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.swt.edu  (login: ftp)}
\thanks{\copyright 2003 Southwest Texas State University.}
\thanks{Published February 28, 2003.}
\vspace{9mm}}

\begin{document}
\setcounter{page}{33}

\title[\hfilneg EJDE/Conf/10\hfil Multistage evolutionary model]
{Multistage evolutionary model for carcinogenesis mutations}

\author[Reza Ahangar \& Xiao-Biao Lin\hfil EJDE/Conf/10\hfilneg]
{Reza Ahangar \& Xiao-Biao Lin}

 \address{Reza Ahangar\hfill\break
 Math Department, Kansas Wesleyan University, Salina, KS 67401, USA}
 \email{rahangar@kwu.edu or rahangar@alltel.net}

\address{Xiao-Biao Lin\hfill\break
 Math Department, North Carolina State University,
 Raleigh, NC 27695-8205, USA}
 \email{xblin@math.ncsu.edu}

\date{}
\thanks{R. Ahangar would like to thank University of Central 
Arkansas, Conway, AR, for giving me \hfill\break\indent
a visiting position, where most of 
this work was performed.} 
\thanks{X.-B. Lin was partially supported by grant DMS-9973105
from National Science Foundation} 
\subjclass[2000]{92C30, 62P10, 92D10, 92D25}
\keywords{Cancer, mutation, interaction, travelling waves,
singular perturbation.}

 \begin{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.
\end{abstract}
\maketitle

\newtheorem{theorem}{Theorem}[section]
\numberwithin{equation}{section}

 \section{Formulations of Carcinogenesis Mutations}

 \subsection*{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 \textit{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.

 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 \textit{adaptation}. The birth of a new cell with a reprogramming
of
 genetic makeup  for adaptation is called a \textit{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.

 \noindent1. \textbf{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.

 \noindent2. \textbf{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.

 \noindent3. \textbf{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
 \textit{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).
\textit{ 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.

 \noindent4. \textbf{Evolution:} Evolution is another important factor in
 cell proliferation and mutation. It is based on the \textit{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).

 \noindent5. \textbf{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).

 \noindent6. \textbf{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 \textit{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).

 \subsection*{Modeling and Formulation}

 Assume that  cell density depends on the following hypotheses

 \noindent\textbf{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.

 \noindent\textbf{H2: Diffusion:} At every stage of mutation, the rate of
 density changes by diffusion and the rate of
 diffusion at every stage remains constant.

 \noindent\textbf{H3: Cell Proliferation process:} Both
 normal and mutant cells have a limited resource environment. Thus, both
 behave alike logistically within certain parameters.

 \noindent\textbf{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,\dots,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,\dots,n$, where Y$_{0}$ will be the initial stage, $%
 Y_{i}(i=1,2,3,\dots,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
 \begin{equation}
 \frac{\partial Y_{i}}{\partial t}=D_{i}\Delta
Y_{i}+a_{i}Y_{i}(1-\frac{Y_{i}%
 }{K_{i}})+\eta_{i}Y_{i}Y_{i-1}-\mu_{i+1}Y_{i}Y_{i+1},  \tag{1.1}
 \end{equation}
 for $i=1,2,\dots,n-1$.
 In this formulation the constant real numbers $D_{i}$ are diffusion
factors,
 $a_{i}$ and $K_{i}$ are logistic parameters, and $\eta_{i}$ and $\mu_{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 $\eta_{i}$ is the \textit{growth advantage of stage} i from the
 previous stage $i-1$, and the parameter $\mu_{i}$ is the \textit{growth
 advantage of the next stage} $i+1$ from the present stage $i$.

 \subsection*{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\leq $ $i<n$. Obviously the factor $\mu _{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. \smallskip

 \noindent\textbf{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
 \begin{equation}
 \frac{\partial Y_{i}}{\partial t}=D_{i}\Delta
Y_{i}+a_{i}Y_{i}(1-\frac{Y_{i}%
 }{K_{i}})+\eta_{i}Y_{i}Y_{i-1},  \tag{1.2}
 \end{equation}
 for $0<i<n$. \smallskip

 \noindent 2) \textbf{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:
 \begin{equation}
 \frac{\partial Y_{i}}{\partial t}=D_{i}\Delta
Y_{i}+a_{i}Y_{i}(1-\frac{Y_{i}%
 }{K_{i}})-\eta _{i}Y_{i}Y_{i-1}.  \tag{1.3}
 \end{equation}
 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. \smallskip

 \noindent\textbf{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 $\eta_{i}$.  Thus, the equation for the latest stage becomes
 \begin{equation}
 \frac{\partial Y_{i}}{\partial t}=D_{i}\Delta
 Y_{i}-bY_{i}+\eta_{i}Y_{i}Y_{i-1}.  \tag{1.4}
 \end{equation}
 The selection pressure at this stage may cause the cells to metastasize
 for survival and reach the blood vessels. \smallskip

 \noindent\textbf{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}\rightarrow\infty$, thus equation (1.1) becomes
 \begin{equation}
 \frac{\partial Y_{n}}{\partial t}=D_{n}\Delta
 Y_{n}+a_{n}Y_{n}+\eta_{n}Y_{n}Y_{n-1}.  \tag{1.5}
 \end{equation}
 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
\textit{intermediate
 stage or Benign}. The full blown developed malignant mutation can be
 described in the following system of \textit{the n-stage model}, i.e.
 initial, benign, and malignant,
 \begin{equation}
 \begin{gathered}
 \frac{\partial u}{\partial t}=D_{0}\Delta u+a_{0}u(1-\frac{u}{K_{0}})-\mu
 _{1}uY_{1}, \\
 \frac{\partial Y_{i}}{\partial t}=D_{i}\Delta
Y_{i}+a_{i}Y_{i}(1-\frac{Y_{i}%
 }{K_{i}})+\eta_{i}Y_{i}Y_{i-1}-\mu_{i+1}Y_{i}Y_{i+1}, \\
 \frac{\partial M}{\partial t}=D_{n}\Delta M+a_{n}M+\eta_{n}MY_{n-1},
 \end{gathered} \tag{1.6}
 \end{equation}
 where $i=1,2,3,\dots,n-1$. Coefficients $\eta_{i}$ and $\mu_{i}$, for  $
 i=0,1,\dots,n$, are constants and are called \textit{interaction
coefficients}.

 \subsection*{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),
\quad\mbox{in }I\times\Omega
$$
The boundary conditions of system (1.8) will be
 B.C:
$$
\frac{\partial u}{\partial n}=0,\frac{\partial Y_{i}}{\partial n}
 =0,(i=1,2,\dots,n-1),\quad t>0,x\in\partial\Omega.
$$
The boundary conditions are to be interpreted as
\textit{no flux}
 conditions. This means that the mutation is not in the metastasized stage
 and there is no migration of cells across $\partial\Omega$. For all of
these
 stages, $\Omega$ is considered the only habitat for cancer cells.

 \subsection*{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{equation}
 \begin{gathered}
 \frac{\partial u}{\partial t}=D_{0}\Delta u+a_{0}u(1-\frac{u}{K_{0}})-\mu
 _{1}uv, \\
 \frac{\partial v}{\partial t}=D_{1}\Delta v+a_{1}v(1-\frac{v}{K_{1}})+\eta
 _{1}uv-\mu_{2}vw, \\
 \frac{\partial w}{\partial t}=D_{2}\Delta w-a_{2}w+\eta_{2}vw. %
 \end{gathered} \tag{1.7}
 \end{equation}
 The initial conditions are  $u(0)=u_{0},v(0)=v_{0},$and $%
 w(0)=w_{0}$ with no flux on the boundary.


 \section{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.

 \subsection*{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{equation}
 \begin{gathered}
 \frac{\partial u}{\partial t}=D_{1}\Delta u+au(1-\frac{u}{K_{1}})-r_{1}uv,\\
 \frac{\partial v}{\partial t}=D_{2}\Delta v-ev+r_{2}uv,
 \end{gathered}  \tag{2.1}
 \end{equation}
 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,
\begin{gather*}
 U=\frac{u}{K_{1}},\quad V=\frac{r_{1}}{a},t'=at,\quad
 x'=x(D_{2}/a)^{-1/2}=(\frac{a}{D_{2}})^{1/2}x \\
D=\frac{D_{1}}{D_{2}},\quad\gamma=\frac{r_{2}K_{1}}{a},\quad\delta
=\frac{ e}{r_{2}K_{1}},
 \end{gather*}
we obtain
 \begin{equation}
 \begin{gathered}
 U_{t}=DU_{xx}+U(1-U-V), \\
 V_{t}=V_{xx}+\gamma(U-\delta)V.
 \end{gathered}  \tag{2.2}
 \end{equation}

 The necessary condition for the survival of mutant cells in this stage is
$ 0<\delta=\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=\dfrac{D_{1}}{D_{2}}$
is  very small and negligible. Thus,
 \[
 U(x,t)=w(x+ct),\quad 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 $\quad$%
 \begin{equation}
 \begin{gathered}
 cw'=Dw''+w(1-w-z), \\
 cz'=z''+\gamma z(w-\delta).
 \end{gathered}  \tag{2.3}
 \end{equation}

 The second order nonlinear system can be substituted by a first order
system  of ODE with respect to the variable $s$
 \begin{equation}
 \begin{gathered}
 w'=(1/c)w(1-w-z), \\
 z'=y, \\
 y'=cy-\gamma z(w-\delta).
 \end{gathered}  \tag{2.4}
 \end{equation}
 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{equation}
 \begin{gathered}
 w(-\infty)=1,\quad w(\infty)=\delta, \\
 z(-\infty)=0,z(\infty)=1-\delta
 \end{gathered} \tag{2.5}
 \end{equation}
 are called type I solutions. Non-negative solutions of (2.5) satisfying
 \begin{equation}
 \begin{gathered}
 w(-\infty)=0,\quad w(\infty)=\delta, \\
 z(-\infty)=0,z(\infty)=1-\delta
 \end{gathered}\tag{2.6}
 \end{equation}
 are called type II solutions. The following conclusion for $D=0$ can be
 shown (Dunbar 1983).

\begin{theorem} \label{thm2.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\gamma(1-\delta)}$. The travelling wave front
solution  to (2.2) satisfying (2.5) exists when
$c\geq \sqrt{4\gamma (1-\delta)}$.
\end{theorem}

 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 ($\delta,1-\delta)$ 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}$.

 \section{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 \textit{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.

 \noindent\textbf{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.

 \noindent\textbf{C2: Analysis of Multistage with Premalignant Stage:} We
 want to study the evolution of a two-stage mutation model under
 several different environmental conditions.

 \noindent\textbf{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.

 \noindent\textbf{C4: Stability:} We study the instability of the
 latest-stage mutant cells which may cause another mutation.

 \noindent\textbf{C5:} We will measure the migration factors to see the
 abnormality of the mutant cells in transition to a metastasized stage
 mutation. \smallskip

 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 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:
$$
 \mbox{Normal Cell } \longrightarrow \mbox{ Benign }
  \longrightarrow \mbox{ Premalignant }
\longrightarrow \mbox{ 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). \smallskip

 \noindent\textbf{Remark:} The stochastic simulation model when $n=4$ has
 been established by Sherman and Portier (1996).

 \subsection*{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\leq x\leq40$.

 Simulation of the travelling wave solution $U(t)$ over the parameter $D$ is
 depicted in Fig.1 and $V(t)$ over the parameter $\delta$ 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 $\delta$,
 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 $\delta$.

 
 \begin{figure}
 \begin{center}
\includegraphics[width=0.7\textwidth]{fig1.eps}
 \end{center}
 \caption{Simulation of the travelling wave solution U(t) over the parameter
$ 0\leq D\leq3$ measuring the migration factor.}
 \label{F1}
 \end{figure}

 \begin{figure}
 \begin{center}
\includegraphics[width=0.7\textwidth]{fig2.eps}
 \end{center}
 \caption{Simulation of the travelling wave solution V(t) over the
 parameter $0\leq\protect\delta\leq4$. This parameter measures how much the
 environment favors mutant cells.}
 \label{F2}
 \end{figure}

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

\noindent{\small \begin{verbatim}
# 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
#
\end{verbatim}}

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


 \section{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

 \noindent i) \textbf{Changing variables:}
 \begin{equation}
 U  =\frac{u}{K_{0}},V=\frac{v}{K_{1}},W=w.  \tag{4.1}
 \end{equation}
 As a result of this change,
$ u_{xx}=K_{0}U_{xx},\;v_{xx}=K_{1}V_{xx},\;w_{xx}=W_{xx}$.

 \noindent ii) \textbf{Changing the coordinate system
 $(t,x)\rightarrow (\tau,\xi)$:}
 \begin{equation}
 \tau=a_{2}t,\text{ }\xi=\sqrt{a_{2/D_{2}}}x.  \tag{4.2}
 \end{equation}
 According to this transformation
$(U_{t},V_{t},W_{t})=a_{2}(U_{\tau},V_{\tau },W_{\tau})$,
 \[
(u_{xx},v_{xx},w_{xx})=(\frac{a_{2}}{D_{2}})(U_{\xi\xi},V_{\xi\xi},W_{\xi\xi
 }).
 \]
 System (1.7) will be in the form
 \begin{equation}
  \begin{gathered}
 U_{\tau}=\frac{D_{0}}{D_{2}}U_{\xi\xi}+\frac{a_{0}}{a_{2}}U(1-U)-\frac{\mu
 _{1}K_{1}}{a_{2}}UV¸ \\
V_{\tau}=\frac{D_{1}}{D_{2}}V_{\xi\xi}+\frac{a_{1}}{a_{2}}V(1-V)+\frac{\eta
 _{1}K_{0}}{a_{2}}UV-\frac{\mu_{2}}{a_{2}}VW, \\
 W_{\tau}=W_{\xi\xi}-W+\frac{\eta_{2}K_{1}}{a_{2}}VW.
 \end{gathered}\tag{4.3}
 \end{equation}

 \noindent iii) \textbf{Rename the Coefficients:} If the constant
 coefficients are not unusually small or large, redefine them by
 \begin{equation}
 \begin{gathered}
\frac{\mu_{1}K_{1}}{a_{2}}=c,\frac{\eta_{1}K_{0}}{a_{2}}=m,\frac{\mu_{2}}{%
 a_{2}}=n,\frac{\eta_{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{gathered}  \tag{4.4}
 \end{equation}
 For convenience, label our new coordinate system $(t,x)$ so that system
 (4.3) becomes
 \begin{equation}
 \begin{gathered}
 U_{t}=d_{1}U_{xx}+aU(1-U)-cUV, \\
 V_{t}=d_{2}V_{xx}+bV(1-V)+mUV-nVW, \\
 W_{t}=W_{xx}-W+eVW.
 \end{gathered} \tag{4.5}
 \end{equation}

 \noindent iv) \textbf{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).

 \noindent The ratio $\frac{a_{1}}{a_{2}}$ is not unusually large or
small,
 but $\frac{a_{0}}{a_{2}}$ is denoted by $\frac{1}{\varepsilon_{1}}$ (where
$\epsilon_{1}\ll 1,$ implies  smaller values for $a_{2})$.

 \noindent 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
 \begin{equation}
\frac{D_{0}}{D_{2}}=\varepsilon_{2},\;\frac{D_{1}}{D_{2}}=D,\;\frac{a_{0}}{%
 a_{2}}=1/\varepsilon_{1},\;\frac{a_{1}}{a_{2}}=b.  \tag{4.6}
 \end{equation}
 For simplicity but without loosing generality, we will replace
 $ \varepsilon_{i}=\varepsilon$ \ (for $i=1,2$). Accepting all of these
changes  in (i), system (4.3) becomes
 \begin{equation}
 \begin{gathered}
 U_{t}=\varepsilon
U_{xx}+\frac{1}{\varepsilon}U(1-U)-\frac{1}{\varepsilon }%
 \frac{\mu_{1}K_{1}}{a_{0}}UV, \\
 V_{t}=DV_{xx}+bV(1-V)+mUV-nVW, \\
 W_{t}=W_{xx}-W+eVW.
 \end{gathered} \tag{4.7}
 \end{equation}
 Renaming the value $\frac{\mu_{1}K_{1}}{a_{0}}=\gamma,$ system (4.7) leads
 to the system
 \begin{gather}
 \varepsilon U_{t}=\varepsilon^{2}U_{xx}+U(1-U-\gamma V),  \tag{4.8-a}\\[4pt]
 \begin{gathered}
 V_{t}=DV_{xx}+bV(1-V)+mUV-nVW, \\
 W_{t}=W_{xx}-W+eVW.
 \end{gathered}  \tag{4.8-b}
 \end{gather}
 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).

 \noindent- $U,\, V$, and $W$ represent the densities of three stages of
 mutant cells (after rescaling).

 \noindent $W$ is the latest stage and is considered the premalignant
stage.

 \noindent  $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
$eVW$
 represents the harsh interactions with surrounding cells.

 \noindent The parameter $\varepsilon$ in the first equation measures, on
 one hand, the relative movement of $\frac{D_{0}}{D_{2}}=\varepsilon$ ($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/\varepsilon$ getting very large).

 \noindent We use the \textit{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.

 \section{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, $\gamma=m$. System (4.8) will be in the following form
 \begin{gather}
 \quad\varepsilon U_{t}=\varepsilon^{2}U_{xx}+U(1-U-mV),  \tag{5.1-a}\\[4pt]
 \begin{gathered}
 V_{t}=DV_{xx}+bV(1-V)+mUV-nVW, \\
 W_{t}=W_{xx}-W+eVW, \\
 U(0)=U_{0},V(0)=V_{0},W(0)=W_{0}.
 \end{gathered}\tag{5.1-b}
 \end{gather}
 For sufficiently small $\varepsilon,$ 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\neq0$.

 \begin{theorem}[Existence of Traveling Wave Solution] The singularly
 perturbed system (5.1) has a travelling wave solution.
 \end{theorem}

 \begin{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
 \textbf{Fisher-Kolomogorov} model
 \begin{equation}
 U_{t}=\varepsilon U_{xx}+\frac{1}{\varepsilon}U(1-U)  \tag{5.2}
 \end{equation}
 is decoupled from the other two equations and has a travelling wave $U(s)$.
 Dunbar (1983) proved that if $D=0$ or $0<D\leq1,$ travelling wave solutions
 exist for the last two equations  in (5.1-b)
 \begin{equation}
 \begin{gathered}
 V_{t}=DV_{xx}+bV(1-V)-nVW, \\
 W_{t}=W_{xx}-W+eVW.
 \end{gathered} \tag{5.3}
 \end{equation}
 At certain parameter values the wave speeds for $U(s)$ and
  $(V(s),W(s))$ agree. Assume that $\varepsilon=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.
 \end{proof}

 \subsection*{Discussions About Singular Perturbation Method}

 When $\varepsilon\rightarrow0$ 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:

 \noindent I) For $U=0$,
 \begin{equation}
 \begin{gathered}
 V_{t}=DV_{xx}+bV(1-V)-nVW, \\
 W_{t}=W_{xx}-W+eVW.
 \end{gathered} \tag{5.4}
 \end{equation}
 II) For $U=1-V$,
 \begin{equation}
 \begin{gathered}
 V_{t}=DV_{xx}+(b+m)V(1-V)-nVW, \\
 W_{t}=W_{xx}-W+eVW.
 \end{gathered}\tag{5.5}
 \end{equation}
 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,
 \begin{equation}
 U(t,x)=\begin{cases} 0, & \text{for }x<x_{0},\\
 1-V, & \text{for }x>x_{0}, \end{cases}  \tag{5.6}
 \end{equation}
 where at the point $x=x_{0}$, the value of $\widetilde{V}=V(x_{0})$, and
 there is a jump for the value of $U$ from $U=0$ to
 $\widetilde{U}=1- \widetilde{V}$.

 To study the solution of the first equation, after rescaling
 \[
 \xi=(x-x_{0})/\varepsilon,\tau=t/\varepsilon,
 \]
 the internal layer will satisfy the ``\textit{slow variable}'' system
 \begin{equation}
 U_{\tau}=U_{\xi\xi}+U(1-U-m\widetilde{V}).  \tag{5.7}
 \end{equation}
 This equation has a travelling wave front solution connecting $U=0$ to
 $\widetilde{U}=1-V(x_{0})$ with the same speed of the last two equations.

 \section{Concluding Remarks}

 The following is a brief review of our results.

 \noindent\textbf{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 \textit{premalignant stage}.

 \noindent\textbf{2- Mathematical Representations for Stages:} A
 deterministic and macroscopic model representing the evolution of
 carcinogenesis mutation using parabolic PDEs was developed.

 \noindent\textbf{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 \textit{cells by selection
pressures}.

 \noindent\textbf{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$\quad$%
 \begin{equation}
 \begin{gathered}
 \frac{\partial u}{\partial t}=D_{1}\Delta u+au(1-\frac{u}{K_{1}})-r_{1}uv,\\
 \frac{\partial v}{\partial t}=D_{2}\Delta v-ev+r_{2}uv.%
 \end{gathered}  \tag{6.1}
 \end{equation}
When $D=\frac{D_{1}}{D_{2}}$ is very small, we have the following statement.

 \begin{theorem} \label{thm6.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\gamma(1-\delta)}$.
 The travelling wave front solution to (2.2) satisfying (2.5) exists when
$ c\geq\sqrt{4\gamma(1-\delta)}$.
\end{theorem}

 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
($ \delta,1-\delta)$ 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}}\cong0$. The simulation of
the
 travelling wave solution to system (2.2) has been demonstrated in Fig.1
 and Fig.2.

 \noindent\textbf{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{equation}
 \begin{gathered}
 \varepsilon U_{t}=\varepsilon^{2}U_{xx}+U(1-U-mV), \\
 V_{t}=DV_{xx}+bV(1-V)+mUV-nVW, \\
 W_{t}=W_{xx}-W+eVW.
 \end{gathered}  \tag{6.2}
 \end{equation}
is the singular perturbation form. We used the idea of
\textit{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.

 \noindent\textbf{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, $\epsilon$
 measures the mobility of a premalignant cell with respect to the normal
 cell. $\epsilon$ also measures the selection pressure on the premalignant
 mutant cell. By making $\epsilon$ small, we increase the moving capability
 and selection pressures on the premalignant cell. Simulation over the
 variation of $\epsilon$ 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 $\epsilon$ 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).

 
 \subsection*{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:

 \noindent\textbf{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.

 \noindent\textbf{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.

 \noindent\textbf{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.

 \noindent\textbf{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.

 \begin{figure}
 \begin{center}
\includegraphics[width=0.7\textwidth]{fig3.eps}
 \end{center}
 \caption{Simulation of the evolution of normal cell density over the
 variation of epsilon $0\leq k1=\varepsilon\leq1$.}
 \label{F3}
 \end{figure}


 \begin{figure}
 \begin{center}
\includegraphics[width=0.7\textwidth]{fig4.eps}
 \end{center}
 \caption{Simulation of the travelling wave solution $V(t)$ over the
 parameter $0\leq\delta\leq4$. This parameter measures how much the
environment favors mutant cells.}
 \label{F4}
 \end{figure}

 \begin{figure}
 \begin{center}
\includegraphics[width=0.7\textwidth]{fig5.eps}
 \end{center}
 \caption{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\leq
 u_{0}\leq1$ over the parameters  $0\leq D\leq3$ and epsilon $0\leq
k1\leq1$.}
 \label{F5}
 \end{figure}

 \begin{figure}
 \begin{center}
\includegraphics[width=0.7\textwidth]{fig6.eps}
 \end{center}
 \caption{Simulation of the evolution of mutant cell density $U(t)$ over
 parameter epsilon $0\leq k1\leq1$.}
 \label{F6}
 \end{figure}

 \begin{figure}[ht]
 \begin{center}
\includegraphics[width=0.7\textwidth]{fig7.eps}
 \end{center}
 \caption{Simulation of the evolution of the singularly perturbed solution
$W(t)$, premalignant cell density over the parameters $1\leq e\leq4$ and
$0 \leq D\leq3$.}
 \label{F7}
 \end{figure}

 \begin{figure}
 \begin{center}
\includegraphics[width=0.7\textwidth]{fig8.eps}
 \end{center}
 \caption{Simulation of the solution $W(t)$ over the parameters
 $0\leq e\leq4 $ and $0\leq D\leq3$.}
 \label{F8}
 \end{figure}

 \begin{figure}
 \begin{center}
\includegraphics[width=0.7\textwidth]{fig9.eps}
 \end{center}
 \caption{Simulation of the singularly perturbed solution
 $W(t):0\leq x\leq12$
 over the parameters $1\leq e\leq4$ and 0$\leq k1\leq1$.}
 \label{F9}
 \end{figure}

\clearpage

 \begin{thebibliography}{99}

 \bibitem{Bogen 89} Bogen,  ``Cell Proliferation Kinetics and Multistage
 Cancer Risk Model'', Journal of the National Cancer Institute, Vol. 81,
No.  4, Feb. 15, 1989.

 \bibitem{Byrne 1995} 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.

 \bibitem{Byrne 2001} 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.

 \bibitem{Davis 2000} Davis, Brian K.; ``Darwinian Aspects of Molecular
 Evolution at Sublinear Propagation Rates'', Bulletin of Mathematical
Biology
 (2000) 62, 675-694.

 \bibitem{Dopazo 2001} 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.

 \bibitem{Dunbar 83} Dunbar, Steve R.; ``Traveling Wave Solutions of
Diffusive
 Lotka-Volterra Equations'', J. Math. Biology 1983, 17:11-32.

 \bibitem{Dunbar 96} 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.

 \bibitem{Edy and Johnson 89} Maitland, A. Edey and Donald C. Johnson;
 ``Blueprints Solving the Mystery of Evolution'', Penguin Books, 1989.

 \bibitem{Fitzgibbon 1984} Fitzgibbon, W. E.; ``Partial Differential
 Equations and Dynamical Systems'', Research Notes in Mathematics 101,
Pitman
 Advanced Publishing Program 1984.

 \bibitem{Frank 93}  Frank-Kamenetskii, Maxim D.; ``Unraveling DNA''
Translated
 from Russian by Lev Liapin, Published by VCH, 1993.

 \bibitem{Goldberg 97} Goldberg, I. D. and Rosen, E. M.; ``Regulation of
 Angiogenesis'' Birkhauser Verlag, 1997.

 \bibitem{Ken-On 1994} Ken-On, Y., ``Fisher-Type Property of Traveling
Waves
 for Competition-Diffusion Equations'', Research Report CDSNS94-159,
GA-Tech,
 1994.

 \bibitem{Kimbell 93} 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).

 \bibitem{Kirshner 96} Kirschner, Denise; ``Using Mathematics to Understand
 HIV Immune Dynamics'', AMS Notices, Feb, 1996, Vol 43, No.2.

 \bibitem{Leffell 96} Leffell, David J. and Brash, Douglas E; ``Sunlight and
 Skin Cancer'', Scientific America, July 1996.

 \bibitem{Mao 96} 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.

 \bibitem{Moolgavkar 86} Moolgavkar Suresh H., ``Carcinogenesis Modeling:
 >From Molecular Biology to Epidemiology'', Ann. Rev. Public Health 1986,
 7.151-69.

 \bibitem{Michelson and Leith 97} 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.

 \bibitem{Owen and Sherratt 97} 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.

 \bibitem{Pettet 2001} 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.

 \bibitem{Portier-Sherman 2000} 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

 \bibitem{Portier -Smith 2000} Portier, C. J., Smith, M. V.; ``Incorporating
 Observability Thresholds of Tumors into the Two-Stage Carcinogenesis
 Model'', Mathematical Biosciences, Elsevier, 163, (2000) 75-89.

 \bibitem{Raff 96} Raff, Martin C.; ``Death Wish'' The Sciences Published by
 New York Academy of Sciences. August 1996.

 \bibitem{Spouge -Shrager 96} John L. Spouge, Richard Shrager, and D. S.
 Dimitrov,\quad``HIV-1 Infection Kinetics in Tissue Cultures'',
Mathematical
 Biosciences, 138:1-22, Elsevier, March 1996.

 \bibitem{Schinazi 95} Schinazi, Rinaldo; ``On an interaction particle
 modeling an epidemic'', by  Mathematical Biology, Spring Verlag 1996, 34:
 915-925.

 \bibitem{Sherman and Portier 96} Sherman, Claire D. and Portier, Christopher
 J.; ``Stochastic Simulation of a Multistage Model of Carcinogenesis'',
 Mathematical Biosciences, volume 134, Number 1, May 1996.

 \bibitem{Smitalova -Sujan 91} Smitalova, Kristina and Sujan Stefan; ``A
 Mathematical Treatment of Dynamical Models in\ Biological Science,'' By
 Translated from Czechoslovakia, by B. Sleeman, Ellis Harwood 1991.

 \bibitem{Travis-Wiarda 1997} Wiarda, D.,and Travis, C.C.; ``Determinibality
 of Model Parameters in a Two-Stage Deterministic Cancer Model''
Mathematical
 Biology, 1997.

 \bibitem{Weinberg 96} Weinberg, Robert; ``How Cancer Arises\dots.'',
 Scientific American, September 1996.

 \bibitem{Weinberg 1998} Weinberg, Robert; ``One Renegade Cell, How Cancer
 Begins'', The Science Masters Series Publisher, 1998.

 \bibitem{Witting 2001} Witting, Lars; ``Population Cycles Caused by
Selection
 by Dependent Competitive Interactions'', Bulletin of Mathematical Biology,
 (2000), 62, 1109-1136.

 \bibitem{Wu 96} Wu, Jianhong; ``Theory and Applications of Partial
 Functional Differential Equations,'' Applied Mathematical Sciences 119,
 Springer 1996.

 \bibitem{Yanagida 1994} Yanagida, E. and Ei, S. I.; ``Dynamics in
 Competition-Diffusion Systems'', SIAM J. Appl. Math., vol. 54, No. 5, pp.
 1355-1373, October 1994.

 \end{thebibliography}

 \end{document}