? Author to whom correspondence should be addressed.

E-mail: torquato@matter.princeton.edu
-- T. S. Deisboeck, M.D. is also a$liated with the Depart-
ment of Neurosurgery, University of Munich (Germany).
J. theor. Biol. (2000) 203, 367}382
doi:10.1006/jtbi.2000.2000, available online at http://www.idealibrary.com on
Simulated Brain Tumor Growth Dynamics
Using a Three-Dimensional Cellular Automaton
A. R. KANSAL*, S. TORQUATO*-?, G. R. HARSH IVA,
E. A. CHIOCCA#B** AND T. S. DEISBOECK#B**--
*Department of Chemical Engineering, Department of Chemistry, -Princeton Materials Institute,
Princeton ;niversity, Princeton, NJ 08544, ;.S.A., ADepartment of Neurosurgery, Stanford ;niversity
Medical School, Stanford, CA 94305, ;.S.A., #Neurosurgical Service, BBrain umor Center,
**Molecular Neuro-Oncology aboratory, Massachusetts General Hospital East,
Harvard Medical School, Charlestown, MA 02129, ;.S.A.
(Received on 20 August 1999, Accepted in revised form on 14 January 2000)
We have developed a novel and versatile three-dimensional cellular automaton model of
brain tumor growth. We show that macroscopic tumor behavior can be realistically
modeled using microscopic parameters. Using only four parameters, this model simulates
Gompertzian growth for a tumor growing over nearly three orders of magnitude in radius. It
also predicts the composition and dynamics of the tumor at selected time points in agreement
with medical literature. We also demonstrate the #exibility of the model by showing the
emergence, and eventual dominance, of a second tumor clone with a di!erent genotype. The
model incorporates several important and novel features, both in the rules governing the
model and in the underlying structure of the model. Among these are a new de"nition of how
to model proliferative and non-proliferative cells, an isotropic lattice, and an adaptive grid
lattice.
2000 Academic Press
1. Introduction
The incidence of primary malignant brain tumors
is already 8/100 000 persons per year and is still
increasing. The vast majority (80%) consists of
high-grade malignant neuroepithelial tumors
such as glioblastoma multiforme (GBM) (Fig. 1)
(Annegers et al., 1981; Werner et al., 1995). In
spite of aggressive conventional and advanced
treatments, the prognosis remains uniformly fatal
with a median survival time for patients with
GBM of 8 months (Black, 1991; Whittle, 1996).
The main reason for this grim outcome is not
only the rapid tumor growth but especially the
fact that, long before the neoplasm can be diag-
nosed, it has already grossly invaded the sur-
rounding brain parenchyma, rendering surgical
removal virtually ine!ective. The proposed se-
quence of proliferation, then invasion, followed
by proliferation again suggests that invasive cells
left behind after an operation not only cause
parenchyma destruction, but also eventually
tumor recurrence (Suh & Weiss, 1984; Burger
et al., 1988; Nazzaro & Neuwelt, 1990; Silbergeld
& Chicoine, 1997). Anti-proliferative treatments
fail because of poor transport across the
0022}5193/00/080367#16 $35.00/0 2000 Academic Press
FIG. 1. T1-contrast enhanced brain MRI-scan showing
a right frontal GBM tumor. Perifocal hypointensity is
caused by signi"cant edema formation. The hyper-intense,
white region (ring-enhancement) re#ects an area of extensive
blood}brain/tumor barrier leakage. Since this regional
neovascular setting provides tumor cells with su$cient nu-
trition it contains the highly metabolizing e.g. dividing,
tumor cells (Selker et al., 1982; Burger et al., 1983; Kelly
et al., 1987; Earnest IV et al., 1988). Therefore, this area
corresponds to the outermost concentric shell of highly
proliferating neoplastic cells in our model (see Fig. 6).
blood}brain barrier, acquired treatment resist-
ance, and lack of susceptibility of single invading
cells due to their minimal proliferative activity
(Giese et al., 1996; Schi!er et al., 1997).
The rapid growth and resilience of tumors
make it di$cult to believe that they behave as
random, disorganized and di!use cell masses and
suggests instead that they are emerging, opportun-
istic systems. If this hypothesis holds true, the
growing tumor and not only the single cell
(Kraus & Wolf, 1993) must be investigated and
treated as a self-organizing complex dynamic
system. This cannot be done with currently avail-
able in vitro/in vivo models or common math-
ematical approaches. Therefore, there is a need
for novel computational models to simulate the
mechanistic complexity of solid tumor growth
and invasion, combining a range of disciplines
including medical, engineering and statistical
physics research.
Here we report the "rst step towards a
complex self-organizing tumor model focusing
on volumetric growth of a solid tumor as the
necessary precondition for subsequent and ongo-
ing invasion. We have developed a three-dimen-
sional cellular automaton (CA) model which
describes tumor growth as a function of time. The
algorithm takes into account that this growth
starts out from a few cells, passes through a multi-
cellular tumor spheroid (MTS) stage (Fig. 2) and
proceeds to the macroscopic stages at clinically
designated time-points for a virtual patient:
detectable lesion, diagnosis and death (Suther-
land, 1988; Mueller-Klieser, 1997). This simpli"-
ed growth approach models macroscopic GBM
tumors as an enormous idealized MTS, mathe-
ematically described by a Gompertz function
(Brunton & Wheldon, 1977; Vaidya & Alexandro
Jr., 1982; Norton, 1988; Marusic et al., 1994).
Modeling the ideal tumor at every step as an
oversized spheroid is especially suited for GBM,
since this tumor, like a large MTS, comprises
large areas of central necrosis surrounded by
a rapidly expanding shell of viable cells (Fig. 1).
In accordance with experimental data, the algo-
rithm also implicitly takes into account that in-
vasive cells are continually shed from the tumor
surface. Future work will treat the invasive dy-
namics explicitly.
The simulation contains several features novel
to the simulation of tumor growth:
z The ability of cells to divide is treated in
a new manner. By rede"ning the transition
between dividing and non-dividing cells, it is
possible to generate medically realistic over-
all tumor growth rates using a biologically
reasonable cell-doubling time.
z Previous researchers have used the Voronoi
tessellation in histopathological image anal-
ysis of tumors (Preston Jr. & Siderits, 1992;
Kiss et al., 1995; Haroske et al., 1996). How-
ever, the present work represents the "rst use
of the Voronoi tessellation (described later in
the text) to study the dynamics of tumor
368 A. R. KANSAL E A.
FIG. 2. MTS-gel assay showing a central spheroid with multiple &&chain''-like invasion pathways leading towards the
boundary (magni"cation: ;200).
growth in a cellular automaton. This tessel-
lation is isotropic in space and hence does
not create the arti"cial anisotropies possible
with square or cubic lattices, which have
typically been used in tumor simulations.
z The model uses a varying density of lattice
sites (an adaptive grid lattice). This allows
small tumors to be simulated with greater
accuracy, while still allowing the tumor to
grow to a large size. Using this variable-den-
sity lattice, tumors can be simulated over
nearly three orders of magnitude in radius.
A CA model treats the discrete nature of actual
cells realistically, giving it great adaptability in
treating complex situations (Kau!man, 1984;
Wolfram, 1984). For example, the addition of
a blood vessel or other growth promoting hetero-
geneity could be studied by altering a single para-
meter. Similarly, the e!ect of the surgical removal
of a portion of the tumor could be readily
modeled. Most importantly, a discrete model
readily allows the inclusion of a number of dis-
tinct subpopulations, corresponding to di!erent
cell clones. In addition, this work is intended to
serve as a preliminary step to modeling systemic
invasive growth dynamics, the nature of which is
ideally suited to study using discrete modeling.
Finally, each site in a CA model can be thought
of as a group of actual cells. This interpretation
allows the model to serve as an intuitive comp-
lement to the results obtained from a continuum
model.
The simulation is designed to predict clinically
important criteria such as the fraction of the
tumor which is able to divide (GF), the non-
proliferative (G
"
/G

arrest) and necrotic frac-

tions, as well as the rate of growth (volumetric
doubling time) at given radii. The simulation
results re#ect a test case derived from the medical
literature very well, proving the hypothesis that
macroscopic tumor growth behavior may be
modeled with primarily microscopic data. Cur-
rent limitations and potential implications of this
model for further tumorigenesis research are
discussed. Since an approach claiming to model
malignant tumor complexity must take into ac-
count cell invasion, we will also brie#y describe
TUMOR DYNAMICS VIA CELLULAR AUTOMATON 369
initial theoretical considerations to model tumor
cell invasion.
In summary, our model is characterized by
several biologically important features:
z The model is able to grow the tumor from
a very small size of roughly 1000 real cells
through to a fully developed tumor with 10
cells. This allows a tumor to be grown
from almost any starting point, through to
maturity.
z The thickness of di!erent tumor layers,
i.e. the proliferative rim and the non-prolif-
erative shell, are linked to the overall tumor
radius by a 2/3 power relation. This re#ects
a surface-area-to-volume ratio, which can be
biologically interpreted as nutrients di!using
through a surface.
z Our inclusion of mechanical con"nement
pressure enables us to simulate the physiolo-
gical con"nement by the skull at di!erent
locations within the brain di!erently.
z The discrete nature of the model and the
variable density lattice allow us to control
the inclusion of mutant &&hot spots'' in the
tumor, i.e. consider genetic instability and
emergence of clonal subpopulations. The
variable density lattice will allow us to look
at such an area at a higher resolution.
In the following section (Section 2), we outline
some of the earlier work that has been done in the
"eld of tumor modeling. Section 3 discusses
the details of the procedure for the simulation.
A summary of our results is contained in Section
4. This is followed by a discussion and concluding
remarks regarding our current work, as well as
future work, in Section 5.
2. Previous Work
Signi"cant research has been done in the
modeling of tumors using theoretical models and
computer simulations. Previous modeling has
been done on a range of tumor behaviors, includ-
ing proliferative growth of the tumor core (see
references immediately below), invasive growth
(Tracqui, 1995; Perumpanani et al., 1996) and
immune response (Sherratt & Nowak, 1992).
Here we will focus on works which have explicitly
modeled the proliferative growth of the tumor
core.
Some of the earliest work in modeling of tu-
mors using a three-dimensional cellular automa-
ton on a cubic lattice was carried out by
DuK chting & Vogelsaenger (1985) for very small
tumors. These automaton rules were designed to
re#ect nutritional needs for tumor growth. Other
important factors, such as surrounding cells
and mechanical pressure, however, remained un-
considered.
Qi et al. (1993) considered a two-dimensional
cellular automaton tumor model that reproduced
idealized Gompertz results. However, cells could
only divide if one of their nearest neighbors was
empty. This created an unrealistically small frac-
tion of a tumor which may divide. Furthermore,
dead tumor cells were assumed to simply dissolve
away rather than accumulating into a necrotic
core, as is seen in real tumors. In addition, the
transition from dividing cells to the resting state
was handled in a purely stochastic manner,
rather than the more biologically reasonable
nutrient-based method used by DuK chting and
Vogelsaenger.
Recent work by Smolle & Stettner (1993)
showed that the macroscopic behavior of a
tumor can be a!ected by the presence of growth
factors at the microscopic level and added the
concept of cellular migration to the behavior of
the cells. The work done by Smolle and Stettner,
however, was qualitative and designed to show
the range of behaviors obtainable from a simple
model. Recent work has been carried out to esti-
mate the model parameters needed to generate
a given con"guration (Smolle, 1998). This model,
as well as those designed by Qi et al. and DuK cht-
ing and Vogelsaenger, relied on an underlying
square (or cubic) lattice. While this provides
a simple method of organizing the automata, it
introduces undesirable asymmetries and other
arti"cial lattice e!ects.
Other researchers have attempted to study the
growth patterns of tumors in a more macroscopic
fashion. The work of Wasserman et al. (1996)
used a "nite element analysis technique to de-
scribe the macroscopic behavior of a tumor
based on stresses imposed by various factors.
Wasserman's approach is similar to that of Chap-
lain & Sleeman (1993), in which the authors used
370 A. R. KANSAL E A.
nonlinear elasticity theory to model a tumor. In
that paper, the growth of the tumor was governed
by a strain-energy function.
Another approach that has been taken by
a number of researchers is to create equations
which describe the tumor phenomenologically.
The best known is the Gompertz model, which
describes the volume, <, of a tumor vs. time, t, as
<"<
"
exp

A
B
(1!exp (!Bt))

, (1)
where <
"
is the volume at time t"0 and A and
B are parameters (Steel, 1977). Qualitatively, this
equation gives exponential growth at small times
which then saturates at large times (decelerating
growth). Many tumors, however, contain more
than one clonal population. A second population
might have di!erent division rates and nutri-
tional needs, resulting in competitive e!ects
which cannot be accounted for by the Gompertz
model. The Jansson}Revesz equations model the
interaction of two populations in a competitive
setting (Cruywagen et al., 1995). These equations
are essentially the classical Lotka}Volterra equa-
tions of logistic growth with an added term to
account for the conversion of one species to the
other. Work by Cruywagen et al. has sought to
apply the Jansson}Revesz equations to tumor
growth. In addition, they have included a dif-
fusive term to each equation, to account for pass-
ive cellular motion. These mathematical models
are useful to describe the general size of a tumor
under relatively simple conditions (two popula-
tions, Fickian di!usion), but have yet to be ex-
tended to multiple populations or active cellular
motility. A review of several other mathematical
models is contained in Marusic et al. (1994).
Finally, the growth of the tumor core has been
modeled in a continuum setting using di!erential
equations with explicit spatial dependence by
several researchers. Adam (1986) used an ordi-
nary di!erential equation, which re#ects mass
conservation of tumor cells, coupled with a reac-
tion-di!usion equation, re#ecting the distribu-
tion of nutrients within the tumor. Ward & King
(1997) used nonlinear partial di!erential equa-
tions to generate pro"les for growth of an avascu-
lar tumor based on a nutrient distribution. More
recently, Byrne & Chaplain (1998) have proposed
a model which also includes the e!ects of nutrient
di!usion and incorporates both apoptosis and
necrosis explicitly. Other models, notably those
by Tracqui et al. (1995) and Woodward et al.
(1996), attempt to model both tumor growth and
the impact of treatment.
3. Simulation Procedure
The underlying lattice for our algorithm is the
Delaunay triangulation, which is the dual lattice
of the Voronoi tessellation (Okabe et al., 1992).
The lattice generation is accomplished by "rst
choosing a set of sites, which are simply points
distributed in space according to some random
process. Each cell in the "nal Voronoi lattice will
contain one of these sites. The cell is then de"ned
by the region of space nearer to that particular
site than to any other site. In two spatial dimen-
sions, a Voronoi lattice is a collection of such
polygons, which "ll the plane. This is illustrated
in Fig. 3. In the three-dimensional analog, the
Voronoi cells take the form of polyhedra.
The Delaunay lattice is generated from the
Voronoi network by connecting those sites
whose polyhedra (in three dimensions) share a
common face. This determines which sites are
nearest neighbors of one another. All references
to nearest neighbors in this paper refer to those
determined in this way. A two-dimensional ana-
log of this is also depicted in Fig. 3.
For our algorithm, a prescribed number of
random points in space are generated. The speed
and memory requirements of the subsequent pro-
grams are strongly dependent on the number of
points chosen. In addition, however, the larger
the number of points chosen, the more accurate
the model. As such, the largest number of points
that can be accommodated in computer memory
and run in a reasonable length of time has been
used.
Because a real brain tumor grows over several
orders of magnitude in volume, the lattice was
designed with a variable grid size. In our lattice,
the density of site was allowed to vary continu-
ously with the radius of the tumor. The density of
lattice sites near the center of the lattice was
signi"cantly higher than that at the edge. A high-
er site density corresponds to less real cells per
TUMOR DYNAMICS VIA CELLULAR AUTOMATON 371
FIG. 3. Two-dimensional space tiled into Voronoi cells.
Points represent sites and lines denote boundaries between
cells. (a) and (b) depict a very small section of a lattice.
(a) shows the Voronoi cells, while (b) shows both the
Voronoi cells, along with the Delaunay tessellation. (c) and
(d) show a more representative section of the lattice, with the
variable density of sites evident. (c) shows the entire lattice
section, (d) shows the same section with the darkened cells
representing a tumor.
automaton cell, and so to a higher resolution.
The higher density at the center enables us to
simulate the #at small-time behavior of the Gom-
pertz curve. In the current state of the model, the
innermost automaton cells represent roughly 100
real cells, while the outermost automaton cells
represent roughly 10" real cells. The average dis-
tance between lattice sites was described by the
following relation:
"
1
6
r``, (2)
in which is the average distance between lattice
sites and r is the radial position at which the
density is being measured. This relation restricts
the increase in the number of proliferating cells as
the tumor grows. Note the 2/3 appearing in the
exponent, which is intended to re#ect a surface-
area-to-volume-like relation. This conforms to
the di!usion of nutrients through the surface
of the tumor, which is known to be a crucial
factor in governing the tumor's growth dynamics
(Folkman & Hochberg, 1973).
While an isotropic Voronoi tessellation can be
generated from any list of random points, not all
point sets are equally reasonable. A purely ran-
dom distribution of points (the Poisson distribu-
tion) will have regions in which the density of
points is very high, corresponding to a very small
Voronoi cell, and regions with a very low density,
corresponding to very large cells. While some
variation in the size and shape of cells is impor-
tant to ensure isotropy, it is biologically unreas-
onable to have large variations. To solve this
problem, a technique common in statistical phys-
ics known as the random sequential addition
(RSA) process was used (Cooper, 1988). In this
technique, as random points are generated, they
are tested to ensure that they are not within
a given distance of any other point. This elimin-
ates the possibility of high-density clusters,
though they can be added later to simulate highly
proliferative tumor &&hot spots''. In addition, it is
well known that in the RSA process there is
a maximum fraction of space which can be "lled.
In three-dimensions, this &&jamming limit'' corres-
ponds to an occupied volume fraction of 0.38. By
approaching this density, the possibility of a large
low-density area is also eliminated. This tech-
nique was adapted for our work to allow the
minimum distance between points to vary, giving
a list of points suitable for generating a biolo-
gically reasonable adaptive lattice. The formula
used to do so was
R
Q
"0.146r``, (3)
where R
Q
is the minimum distance between points
at distance r from the lattice center.
The list of random points in space was fed to
a program written by Ernst Mucke called detri
(Mucke, 1997). The detri program generates the
Delaunay triangulation in three spatial dimen-
sions for a given list of points. To test the detri
program, a second Delaunay program, del-tree3,
was also downloaded (Devillers, 1996). Both pro-
grams were run on the same list of points and the
results compared to ensure that they were identi-
cal. Because the detri program provided a more
convenient form of output (listing tetrahedra by
372 A. R. KANSAL E A.
FIG. 4. An idealized tumor. A detailed description of the
di!erent regions is contained in the text.
the points' indices rather than the speci"c coordi-
nates) it was used for the remainder of the work
done.
The detri program was run on an IBM SP2
parallel computer. Due to memory limitations,
the list of sites was divided into sublists. The
program was then run on each sublist, generating
partial lattices. Culling these partial lattices into
a single complete lattice required some care. In
order to avoid missing triangles across the divis-
ions between sublists, a signi"cant amount of
overlap was included in each sublist. A second
challenge, however, are the triangles formed at
the edges of each sublist. This led to extra
triangles being included in the overall lattice.
Because these extra triangles are localized at the
subset boundaries, however, it is possible to elim-
inate them. This is done by generating two di!er-
ent sets of sublists from the same original list and
then comparing the two "nal lattices. Any tri-
angles not found in both lattices are considered
to be e!ects of the division and are discarded.
This produces a Delaunay lattice over the entire
point set. In order to check the procedure fol-
lowed here, the division process was done on
a point set that could be run in its entirety. The
division method gave the same "nal lattice as
when the entire point set was run and so was
considered to be accurate.
Once the lattice is generated, the proliferation
algorithm can be run. This algorithm is designed
to allow a tumor consisting of a few automaton
cells, representing roughly 1000 real cells, to grow
to a full macroscopic size. An idealized model of
a macroscopic tumor is a spherical body consist-
ing of several concentric shells. The inner core,
the gray region in Fig. 4, is composed of necrotic
cells. The necrotic region has a radius R
L
, which
is a function of time, t, and is characterized by its
distance from the proliferation rim,
L
. The next
shell, the cross-hatched region in the "gure, con-
tains cells which are alive but in the G
"
cell-cycle
rest state. This is termed the non-proliferative
region and is de"ned in terms of its distance from
the edge of the tumor,
N
. This thickness is the
maximum distance from the tumor edge with
a high enough nutrient concentration to main-
tain active cellular division. In real MTS tumors,
however, only about one-third of the viable cells
increase the tumor size by proliferation [which is
numerically comparable to the growth fraction in
macroscopic tumors (Hoshino & Wilson, 1975)].
The others are actively dividing, but the new cells
leave the central tumor mass and supposedly
trigger invasion into the surrounding tissue
(Landry et al., 1981; Freyer & Schor, 1989). Fi-
nally, as discussed below, an individual cell can
only divide if free space exists within a certain
distance of it. This distance must also be
N
(as
de"ned above) to properly account for the nutri-
ent gradient basis for the transition of cells be-
tween the actively dividing and G
"
arrested
states. This distance is depicted as a small broken
circle in Fig. 4. Both real tumors and our
simulated tumors, however, are not perfectly
spherical. As such, the values of R
R
and R
L
vary
over the surface of the tumor. The single values
used in the algorithm and listed in our results are
obtained by averaging the radii of all the cells at
the edge of the tumor or of the necrotic region,
respectively, according to the relations:
R
R
"

N
N
G
r
G
N
N
, (4)
R
L
"

N
L
G
r
G
N
L
, (5)
TUMOR DYNAMICS VIA CELLULAR AUTOMATON 373
FIG. 5. Schematic of expansive growth via intercellular
mechanical stress (IMS). In (a) cell A is non-tumorous, B}D
are tumorous and able to divide. (b) Depicts the results of
attempted divisions by cells C and D, in which C is the new
cell created by the successful division of C. A full description
is contained in the text.
where N
N
denotes the number of cells on the edge
of the proliferative region and N
L
denotes the
number of cells on the edge of the necrotic core.
This de"nition of proliferative cells allows cells
that are not on the immediate edge of the tumor
to divide. In order to accommodate these divis-
ions, without allowing discontinuous division, an
algorithm which allows for expansive growth via
intercellular mechanical stress (IMS) is used.
A one-dimensional schematic of IMS growth is
depicted in Fig. 5. In part (a) of the "gure, cell A is
a non-tumorous cell. Non-tumorous cells are
treated as empty in the current model, which
simply means that, when &&"lled'', they are con-
sidered to have been forced into the surrounding
region of indistinguishable non-tumorous cells.
Cells B, C, and D are tumor cells which may
divide. When cell D attempts to divide, it cannot
"nd an empty space within
N
and will turn non-
proliferative. When cell C attempts to divide it
can "nd cell A and so it will divide. Note that,
biologically, cell C actually cannot &&see'' cell A,
but rather senses its location because of the high-
er nutrient level. Division creates a cell C which
"lls the space previously occupied by cell B,
which is in turn forced to the space previously
"lled by cell A. Cell A is regarded as being forced
into the surrounding tissue and disappears (for
now) from our consideration. This is the new
con"guration depicted in part (b) of Fig. 5.
In three spatial dimensions, 1.5 million lattice
sites are used. This has been found to be the
minimum required to give adequate spatial res-
olution over the entire range of tumor radii. Nat-
urally, more sites would be desirable, however,
practical limits dictate that as few sites as possible
be used. The initial tumor is a few automaton
cells, representing roughly 1000 real cells, located
at the center of the lattice.
In summary, the four key quantities R
R
,
N
,
L
,
and p
B
are functions of time calculated within the
model. To "nd them, the simulation utilizes four
microscopic parameters: p
"
, a, b and R
K?V
. These
parameters are linked to the cell-doubling time,
the nutritional needs of growth-arrested cells, the
nutritional needs of dividing cells, and the e!ects
of con"nement pressure, respectively. The quant-
ities are calculated according to the following
algorithm.
z Initial setup: The cells within a "xed initial
radius of the center of the grid are designated
proliferative. All other cells are designated as
non-tumorous.
z Time is discretized and incremented, so that
at each time step:
*Each cell is checked for type: non-tumor-
ous or (apoptotic and) necrotic, non-pro-
liferative or proliferative tumorous cells.
*Non-tumorous cells and tumorous nec-
rotic cells are inert.
*Non-proliferative cells more than a cer-
tain distance,
L
, from the tumor's edge
are turned necrotic. This is designed to
model the e!ects of a nutritional gradi-
ent. The edge of the tumor is taken to be
the nearest non-tumorous cell.

L
"aR``
R
, (6)
where a, the base necrotic thickness, is
a parameter with units of (length)`.
Note the 2/3 in the exponent, again
indicating a surface-area}volume-type
relation.
374 A. R. KANSAL E A.
TABLE 1
Summary of time-dependent functions and input parameters
for our model
Functions within the model (time dependent)
R
R
Average overall tumor radius

N
Proliferative rim thickness (determines growth fraction)

L
Non-proliferative thickness (determines necrotic fraction)
p
B
Probability of division (varies with time and position)
Parameters (constant inputs to the model)
p
"
Base probability of division, linked to cell-doubling time
a Base necrotic thickness, controlled by nutritional needs
b Base proliferative thickness, controlled by nutritional needs
R
K?V
Maximum tumor extent, controlled by pressure response
*Proliferative cells are checked to see if
they will attempt to divide. This is a ran-
dom process, though the probability of
division, p
B
, is in#uenced by the location
of the dividing cell (r), re#ecting the ef-
fects of mechanical con"nement pressure.
This e!ect requires the use of an addi-
tional parameter, the maximum tumor
extent, R
K?V
. This probability is deter-
mined by the equation:
p
B
"p
"

1!
r
R
K?V

. (7) .
*If a cell attempts to divide, it will search
for su$cient space for the new cell begin-
ning with its nearest neighbors and ex-
panding outwards until either an empty
(non-tumorous) space is found or noth-
ing is found within the proliferation
radius,
N
. The radius searched is
calculated as

N
"bR``
R
, (8)
where b, the base proliferative thickness,
is a parameter with units of (length)`.
*If a cell attempts to divide but cannot
"nd space it is turned into a non-prolif-
erative cell.
(The above two steps constitute the re-
de"nition of the proliferative to non-
proliferative transition that is one of the
most important new features of the
model. They allow a larger number of
cells to divide, since cells no longer need
to be on the outermost surface of the
tumor to divide. In addition, it ensured
that cells which cannot divide are cor-
rectly labeled as such.)
z After a predetermined amount of time has
been stepped through, the volume and radius
of the tumor can be plotted as a function of
time.
z The type of cell contained in each grid can
also be saved at given times.
Table 1 summarizes the important time-depen-
dent functions calculated by the algorithm and
the constant parameters used.
4. Results
The simulation has been compared with avail-
able experimental data for an untreated GBM
tumor from medical literature. The parameters
compared were cell number, growth fraction,
necrotic fraction and volumetric doubling time.
Medically, these data are used to determine a tu-
mor's malignancy and the prognosis for its future
growth. Because it is impossible to determine the
exact time a tumor began growing, the medical
data are listed at "xed radii. The di!erent cell
fractions used were extrapolated from the sphe-
roid level and compared to data published for
TUMOR DYNAMICS VIA CELLULAR AUTOMATON 375
TABLE 2
Comparison of test case data and simulation results (Sim). Note that the time row is
simulation data only and is taken from the start of the simulation not from the
theoretical start of the tumor growth
Spheriod Detect. Diagnosis Death
lesion
Time Sim. Day 69 Day 223 Day 454 Day 560
Radius Data 0.5 mm 5 mm 18.5 mm 25 mm
Sim. 0.5 5 18.5 25
Cell no. Data 10" 10" 5;10" 10
Sim. 7;10` 6;10` 4;10" 9;10"
Growth fraction Data 36% 30% 20% 9%
Sim. 35 30 18 11
Necrotic fraction Data 46% 49% 55% 60%
Sim. 38 53 58 63
Volume-doubling time Data 6 days 45 days 70 days 105 days
Sim. 9 36 68 100
cell fractions at macroscopic stages. Previous
research has shown that the expanding tumor
increases both its cell loss (through necrosis/apo-
ptosis and invasion) and its quiescent cell popula-
tion (G
"
/G

arrested) due to a declining gradient

of nutritional elements towards the center of the
rapidly growing avascular mass (Folkman &
Hochberg, 1973; Durand, 1976; Landry et al., 1981;
Freyer & Sutherland, 1986; Mueller-Klieser et al.,
1986; Rotin et al., 1986; Freyer & Schor, 1989). At
advanced tumor stages, volumetric growth slows
down mostly due to the declining growth fraction
(GF), an increasing cell loss and G
"
-fraction
(Turner & Weiss, 1980; Bauer et al., 1982; Freyer
& Sutherland, 1985) caused by a lack of nutrition
as well as an increasing con"nement pressure.
Summarized in Table 2 is the comparison
between simulation results and data (experi-
mental, as well as clinical) taken from the medical
literature. For macroscopic glial tumors (the
latter three time points) cell numbers, tumor
volumes, and volumetric doubling time measure-
ments, including growth and necrotic fractions
and survival times, were taken from Hoshino
& Wilson (1975); Hoshino & Wilson (1979),
Yamashita & Kuwabara (1983), Yoshii et al.
(1986), Alvord (1995), Blankenberg et al. (1995),
Pierallini et al. (1996) and Burgess et al. (1997).
For the microscopic spheroid level, doubling
times as well as viable rim diameter, rim cell
fractions, necrotic fractions and cell shedding
data were taken from Haji-Karim & Carlsson
(1978), Carlsson et al. (1983), Carlsson & Acker
(1988), Freyer & Schor (1989) and Landry et al.
(1981). Since it is virtually impossible to measure
total cell numbers in macroscopic tumor (the
latter three time points), we have used the vol-
ume-doubling times to estimate the number of
cells in the tumors, including at the microscopic
(MTS) level. This method, however, leads to
a rather high cell number at the MTS level [ex-
perimental values range from 10" to 10` cells
(Freyer, 1988, 1998)]. In a few cases, glioma MTS
data were not available and thus data from other
well-characterized cell lines were used after care-
ful evaluation. On the whole, the simulation data
reproduce the test case very well. The virtual
patient would die roughly 11 months after the
tumor radius reached 5 mm and 3.5 months after
the expected time of diagnosis. The fatal tumor
volume is about 65 cm`.
These data were created using a tumor which
was grown from an initial radius of 0.1 mm. The
following parameter set was used:
p
"
"0.192, a"0.42 mm`, b"0.11 mm`,
R
K?V
"37.5 mm.
This value of p
"
corresponds to a cell-doubling
time of 4.0 days, which is reasonable for high-
grade glial tumors (Hoshino & Wilson, 1979;
376 A. R. KANSAL E A.
FIG. 6. The development of the cross-central section of
a tumor in time. (a) corresponds to the tumor spheroid stage,
(b) to the "rst detectable lesion, (c) to diagnosis and (d) to
death. The dark-gray outer region is comprised of proliferat-
ing cells, the light-gray region is non-proliferative cells and
the black region is necrotic cells. The scales are in milli-
meters.
FIG. 7. Cross sections of a fully developed tumor
(radius"25 mm). (a) The central slice. (b) and (c) are taken
10 mm from the center. (d) Taken 20 mm from the center, on
the same side as (b). The dark-gray region is proliferating
cells, the light-gray region non-proliferative cells and the
black region necrotic cells. The scales are in millimeters.
Pertuiset et al., 1985). The a and b parameters
have been chosen to give a growth history that
quantitatively "ts the test case. As discussed be-
low, the speci"cation of these parameters corres-
ponds to the speci"cation of a clonal strain. This
is manifested in qualitative behavior that is inde-
pendent of the choice of the a and b parameters,
but quantative behavior that is strongly a!ected
by them. The R
K?V
parameter was similarly
chosen to match the test case history. In this case,
however, the "t is relatively insensitive to the
value of R
K?V
, as long as the parameter is some-
what larger than the fatal radius in the test case.
Indeed, the "t is relatively insensitive to the exact
form of eqn (7) in general.
Since a three-dimensional CA image is di$cult
to visualize, cross sections of the tumors are
shown instead. The growth of the tumor can be
followed graphically over time in Fig. 6. It depicts
the central cross section of the tumor, with the
convention that necrotic cells are labeled with
black, non-proliferative tumorous cells with light
gray and proliferative tumor cells with dark gray.
Note that the plotted points are not intended
to depict the exact shape of the Voronoi cells,
but rather just their positions. Figure 7 depicts
several cross sections taken at di!erent positions
from within a single fully developed tumor. As
expected in this idealized case, the tumor is essen-
tially spherical, within a small degree of random-
ness. The high degree of spherical symmetry
ensures that the central cross section is a
representative view. The volume and radius of
the developing tumor are shown vs. time in Fig. 8.
Note that the virtual patient dies while the
untreated tumor is in the rapid growth phase.
While these results con"rm that the algorithm
is able to reproduce a very idealized case, they do
not illustrate its ability to model more complex
situations. As mentioned in the Introduction, one
such complexity is the inclusion of multiple dis-
tinct tumor clones. An interesting question we
have begun to address is that of the emergence of
a clonal population from a small mutated hot
spot in the original cell population. To address
this, a second strain is de"ned by the parameter
set:
p
"
"0.384, a"0.42 mm`, b"0.11 mm`,
R
K?V
"37.5 mm.
TUMOR DYNAMICS VIA CELLULAR AUTOMATON 377
FIG. 8. Plots of the radius and volume of the tumor
versus time. The lines correspond to simulation predictions,
using the "rst parameter set given in the text. The plotted
points re#ect the test case derived from the medical litera-
ture. A quantitative comparison of the simulation with the
test case is given in Table 2.
The second strain is a more rapidly dividing
clone than the "rst, parental strain de"ned above.
The increase in p
"
gives a decreased cell-doubling
time (1.7 vs. 4 days in the primary strain), chosen
to re#ect the degree of heterogeneity between
di!erent clonal strains. This example corres-
ponds to the emergence of a more rapidly divid-
ing clone from the more slowly proliferating
malignant parental strain and thereby models the
&&progression pathway'' of GBM tumors from less
malignant precursors (Lang et al., 1994). This
simulation begins in the same way as the pre-
vious one, with a very small tumor, composed
entirely of cells of the primary strain. This tumor
is then allowed to grow until it reaches a
prede"ned overall radius of 3.8 mm. A mutated
hot spot is then introduced as a single, randomly
selected, CA cell (corresponding to roughly 10`
real cells, i.e. the size of a spheroid) changing from
the primary strain to the secondary strain. This
represents the appearance of the second genotype
in roughly 0.01% of the viable tumor cells at that
time, assuming monoclonal expansion. At this
time we do not distinguish between random
mutational activity (intrinsic) or environmentally
caused mutational stresses. Figure 9 depicts the
development of the secondary strain (shown in
blue) in the cross-sectional slice nearest the ap-
pearance of the mutation, which is displaced from
the tumor center by roughly 2.5 mm in this case.
To aid in distinguishing between strains, cells of
the primary clonal strain (dark gray in Fig. 6)
have been colored red, while the non-prolif-
erative cells (light gray in previous "gures) have
been colored yellow. A bulge near the initial
position of mutation appears soon after the
mutation occurs. This corresponds to a more
rapidly advancing tumor rim due to a higher
local proliferation rate. This higher proliferation
rate also leads to additional cells being shed per
unit of time, which in turn leads to a higher
invasive potential towards this area. Gradually,
the secondary strain overtakes the primary strain
and comes to dominate the entire tumor. This
example shows the spatial evolution of a second-
ary tumor population. From such a simulation, it
is possible to consider the e!ects of mutations on
the shape and location of a tumor. In this simula-
tion, the secondary strain overtakes the primary
strain su$ciently rapidly so that the center of mass
of the entire tumor shifts considerably from its
original location. In addition, the bulge created by
the more rapidly advancing front increases the
surface-area-to-volume ratio of the tumor, there-
by increasing the growth fraction and causing the
tumor growth to accelerate more than would be
the case if the tumor remained spherical.
Currently, we are exploring the survival prob-
ability of a secondary strain as a function of the
number of cells that have mutated and the com-
petitive advantage it enjoys as well as the posi-
tion within the tumor at which the second strain
arises. Given too small an advantage, the second
strain will be quickly overwhelmed by the much
more common primary strain and will likely dis-
appear. Preliminary results show that for a fairly
small mutation, like the one considered here,
a large competitive advantage in doubling time is
required (of the order of a factor of two), for the
secondary population to survive. In biomedical
terms, the tumor is &&searching'' for the best suited
heterogeneic pattern, i.e. the genotypes and
phenotypes (considering gene expression changes,
along with mutations) that are optimal relative to
the tumor's growth rate. The current model
focuses on the importance of growth-promoting
genetic combinations. From a biological per-
spective, the model incorporates selection
pressures leading to the dominance of the clonal
378 A. R. KANSAL E A.
FIG. 9. Cross-section of a tumor, showing the emergence and eventual dominance of a more rapidly growing secondary
strain. The red region corresponds to proliferating cells of the primary strain (darker gray in Figs 6 and 7), the blue region to
those of the secondary strain, while the yellow (light gray in previous "gures) and black regions correspond to non-
proliferative and necrotic cells of either strain, respectively. The cross-sections are taken 2.5 mm from the tumor's central
plane. (a) Depicts the tumor roughly 1 month after the initial mutated CA cell was introduced in the simulation, (b) 5 months
after the mutation, (c) 10 months after the mutation, and (d) 20 months after the mutation. Note that the tumor origin is
located at (0, 0), but because of the mutation, the center of mass of the tumor shifts substantially from this location, before
being forced back by the boundary conditions.
population best suited to the tumor's environ-
ment. Realistically, environmental in#uences
show spatial heterogeneity and as such will not
necessarily lead to the dominance scenario depic-
ted here, but rather may lead to a tumor with
multiple, coexisting tumor clones (i.e. the hall-
mark of GBM tumors), each best suited to their
local environment.
5. Discussion and Conclusions
Substantial progress has been made in the vari-
ous specialized areas of cancer research. Yet the
complexity of the disease on both the single cell
level as well as the multicellular tumor stage has
led to the "rst attempts to describe tumors as
complex, dynamic, self-organizing biosystems,
rather than merely focusing on single features
(Bagley et al., 1989; Schwab & Pienta, 1996;
Co!ey, 1998; Waliszewski et al., 1998). To begin
to understand the complexity of the proposed
system, novel simulations must be developed, in-
corporating concepts from many scienti"c areas
such as cancer research, statistical mechanics,
applied mathematics and nonlinear dynamical
systems.
To our knowledge, this is the "rst three-dimen-
sional cellular automaton model of solid tumor
growth, which realistically models the macro-
scopic behavior of a malignant tumor as a func-
tion of time using predominantly microscopic
parameters. This four-parameter model predicts
the composition and dynamics of malignant
brain tumor proliferation at selected clinically
relevant time points in agreement with experi-
mental and clinical data. From a modeling per-
spective, this is also the "rst use of the Voronoi
tessellation to study tumor growth in a cellular
automaton. The use of a variable density lattice
enables the simulation of tumor growth over
nearly three orders of magnitude in radius.
Finally, the ability of internal cells to divide rep-
resents a physiologically more realistic situation
and changes the proliferative dynamics from pre-
vious models.
In addition, the discrete nature of our model
enables us to directly simulate more complex
physiological situations with only minor alter-
ations. Examples include regional di!erences in
structural con"nement and the in#uence of envir-
onmental factors, such as blood vessels. The im-
pact of surgical procedures or other treatments
can also be modeled. Further, we can currently
include the hallmark of GBM tumors*hetero-
geneity*by using tumorous subpopulations
with di!erent growth behaviors (Giangaspero
& Burger, 1983; Burger & Kleihues, 1989; Paulus
& Pfei!er, 1989). These tumor subpopulations
represent a clear step forward from previous
&&monoclonal'' (single population) models. Envir-
onmental stresses are known to exert a strong
in#uence on tumor growth (Helmlinger et al.,
1997). These stresses can both cause genetic
variations, and so increase the number of sub-
populations, and select for the most "t tumor
subpopulation. It is this diversity in structure and
function which enables tumors to react locally
and globally. It is important to mention that,
according to this concept, systemic tumoral
growth increases both intrinsic and extrinsic
stress and therefore leads to a higher selection
pressure and mutational probability. Since this
may advance tumor progression, it clearly argues
against Foulds' rule III (that tumor growth and
progression are independent) (Foulds, 1954;
Rubin, 1994). In the next iteration of the model,
we will simulate the conditions which create re-
gional genetic instability and study the e!ect of
speci"c mutations in tumor cells on the macro-
scopic growth of tumors. Such information could
allow tumor biopsies, which are currently region-
ally limited, to be put into the context of an
evolving system. This will lead to important in-
formation about intratumoral competition and
cooperation and how tumors adapt to maintain
their "tness in the environment they are present-
ed (Greenspan, 1976; Axelrod & Hamilton, 1981;
Gatenby, 1996).
Another important step on our way to a
complex dynamic tumor model is the speci"c
inclusion of the other key feature of tumors,
namely single cell invasion. The current model
only includes such an e!ect implicitly. Active cell
motility is a crucial feature, not only because of
its local disruptive capacity, but also because of
its signi"cance to treatment. If a solid tumor is
removed, these invasive cells, which would be left
behind, can eventually cause recurrence of the
tumor. We are currently working to theoretically
assess the factors which may drive some of the
TUMOR DYNAMICS VIA CELLULAR AUTOMATON 379
structural elements within the invasive network,
in accordance with a proposed deterministic
chaos principle (Habib et al., unpublished data).
The spatial distribution of these factors will also
depend on the heterogeneous intercellular space.
While this space is not included in the current work,
it will be properly accounted for in further work.
In conclusion, such complex interdisciplinary
models are likely to lead to insights within cancer
research and into complex biosystems in general.
Such insights hold promise for increasing our
understanding of tumors as self-organizing sys-
tems, which in turn may have signi"cant impact
both on cancer research and on clinical practice.
By no means does this model claim such com-
pleteness. It is instead a "rst, but promising, step
towards advanced modeling techniques treating
tumors as complex dynamic systems.
S.T. gratefully acknowledges the Guggenheim
Foundation for his Guggenheimfellowship to conduct
this work. This work has also been supported in part
by grants CA84509 and CA69246 from the National
Institutes of Health. Calculations were carried out on
an IBM SP2, which was kindly provided by the IBM
Corporation (equipment grant to Princeton Univer-
sity for the Harvard}Princeton Tumor Modeling Pro-
ject). The technical assistance provided by Drs Kirk E.
Jordan and Gyan V. Bhanot of the IBM T.J. Watson
Labs is gratefully acknowledged. The authors would
also like to thank Drs Stuart A. Kau!man of the Santa
Fe Institute for Complex Science, Jerome B. Posner
of Memorial Sloan Kettering Cancer Center, Michael
E. Berens of the Barrows Neurological Institute and
H.J. Reulen of the University of Munich Neuro-
surgical Department for valuable discussions.
REFERENCES
ADAM, J. A. (1986). A simpli"ed mathematical model of
tumour growth. Math. Biosci. 81, 229}242.
ALVORD, E. C. (1995). Patterns of growth of gliomas. Am. J.
Neuroradiol. 16, 1013}1017.
ANNEGERS, J. F., SCHOENBERG, B. D., OKAZAKI, H. &
KURLAND, L. T. (1981). Epidemiologic study of primary
intracranial neoplasms. Arch. Neurol. 38, 217}219.
AXELROD, R. & HAMILTON, W. D. (1981). The evolution of
cooperation. Science 211, 1390}1396.
BAGLEY, R. J., FARMER, J. D., KAUFFMAN, S. A., PACKARD,
N. H., PERELSON, A. S. & STADNYK, I. M. (1989).
Modeling adaptive biological systems. Biosystems 23,
113}138.
BAUER, K. D., KENG, P. C. & SUTHERLAND, R. M. (1982).
Isolation of quiescent cells from multicellular tumor sphe-
roids using centrifugal elutriation. Cancer Res. 42, 72}78.
BLACK, P. M. (1991). Brain tumors, Part 2. N. Engl. J. Med.
324, 1555}1564.
BLANKENBERG, F. G., TEGLITZ, R. L., ELLIS, W., SALAMAT,
M. S., MIN, B. H., HALL, L., BOOTHROYD, D. B., JOHN-
STONE, I. M. & ENZMANN, D. R. (1995). The in#uence of
volumetric doubling time, DNA ploidy, and histological
grade on the survival of patients with intracranial as-
trocytomas. Am. J. Neuroradiol. 16, 1001}1012.
BRUNTON, G. F. & WHELDON, T. E. (1977). Prediction of the
complete growth pattern of human multiple myeloma
from restricted intial measurements. Cell issue Kinet. 10,
591}594.
BURGER, P. C. & KLEIHUES, P. (1989). Cytologic composi-
tion of the untreated glioblastoma with implications for
evaluation of needle biopsies. Cancer 63, 2014}2023.
BURGER, P. C., HEINZ, E. R., SHIBATA, T. & KLEIHUES, P.
(1988). Topographic anatomy and CT correlations in
the untreated glioblastoma multiforme. J. Neurosurg. 68,
698}704.
BURGER, P. C., SCHOLD, P. J., SMITH Jr., K. R., ODOM,
G. L., CRAFTS, D. C. & GIANGASPERO, F. (1983).
Computerized tomographic and pathologic studies of
the untreated, quiescent, and recurrent glioblastoma multi-
forme. J. Neurosurg. 58, 159}169.
BURGESS, P. K., KULESA, P. M., MURRAY, J. D. & ALVORD,
E. C. (1997). The interaction of growth rates and di!usion
coe$cients in a three-dimensional mathematical model of
gliomas. J. Neuropathol. Exp. Neurol. 56, 704}713.
BYRNE, H. M. & CHAPLAIN, M. A. J. (1998). Necrosis and
apoptosis: distinct cell loss mechanisms in a mathematical
model of avascular tumour growth. J. heor. Med. 1,
223}235.
CARLSSON, J. & ACKER, H. (1988). Relations between pH,
oxygen partial pressure and growth in cultured cell sphe-
roids. Int. J. Cancer 42, 715}720.
CARLSSON, J., NILSSON, K., WESTERMARK, B., PONTEN,
J., SUNDSTROM, C., LARSSON, E., BERGH, J., PAHLMAN, S.,
BUSCH, C. & COLLINS, V. P. (1983). Formation and growth
of multicellular spheroids of human origin. Int. J. Cancer
31, 523}533.
CHAPLAIN, M. A. J. & SLEEMAN, B. D. (1993). Modelling the
growth of solid tumours and incorporating a method for
their classi"cation using nonlinear elasticity theory.
J. Math. Biol. 31, 431}473.
COFFEY, D. S. (1998). Self-organization, complexity
and chaos: the new biology for medicine. Nat. Med. 4,
882}885.
COOPER, D. W. (1988). Random sequential packing simula-
tions in 3 dimensions for spheres. Phys. Rev. A 38, 522}524.
CRUYWAGEN, G. C., WOODWARD, D. E., TRACQUI, P.,
BARTOO, G. T., MURRAY, J. D. & ALVORD, E. C. (1995).
The modelling of di!usive tumours. J. Biol. Sys. 3,
937}945.
DEVILLERS, O. (1996). http://www.inria.fr:/prisme/logiciel/
del-tree.html
DUG CHTING, W. & VOGELSAENGER, T. (1985). Recent pro-
gress in modelling and simulation of three-dimensional
tumor growth and treatment. Biosystems 18, 79}91.
DURAND, R. E. (1976). Cell cycle kinetics in an in vitro
tumor model. Cell issue Kinet. 9, 403}412.
EARNEST IV, F., KELLY, P. J., SCHEITHAUER, B. W., KALL,
B. A., CASCINO, T. L., EHMAN, R. L., FORBES, G. S. &
AXLEY, P. L. (1988). Cerebral astrocytomas: his-
topathological correlation of MR and CT contrast en-
hancement with stereotactic biopsy. Radiology 166,
823}827.
380 A. R. KANSAL E A.
FOLKMAN, J. & HOCHBERG, M. (1973). Self-regulation
of growth in three dimensions, J. Exp. Med. 138,
745}753.
FOULDS, L. (1954). The experimental study of tumor pro-
gression: a review. Cancer Res. 14, 327}339.
FREYER, J. P. (1988). Role of necrosis in regulating the
growth saturation of multicellular spheroids. Cancer Res.
48, 2432}2439.
FREYER, J. P. (1998). Decreased mitochondrial function in
quiescent cells isolated from multicellular tumor spheroids.
J. Cell. Physiol. 176, 138}149.
FREYER, J. P. & SCHOR, P. L. (1989). Regrowth kinetics of
cells from di!erent regions of multicellular spheroids of
four cell lines. J. Cell. Physiol. 138, 384}392.
FREYER, J. P. & SUTHERLAND, R. M. (1985). A reduction in
the in situ rates of oxygen and glucose consumption of cells
in the EMT6/Ro spheroids during growth. J. Cell. Physiol.
124, 516}524.
FREYER, J. P. & SUTHERLAND, R. M. (1986). Regulation
of growth saturation and development of necrosis in
EMT6/Ro multicellular spheroids by the glucose and oxy-
gen supply. Cancer Res. 46, 3504}3512.
GATENBY, R. A. (1996). Applications of competition theory
to tumour growth: implications for tumour biology and
treatment. Eur. J. Cancer 32A, 722}726.
GIANGASPERO, F. & BURGER, P. C. (1983). Correlations
between cytologic composition and biologic behavior in
the glioblastoma multiforme: a postmorterm study of 50
cases. Cancer 52, 2320}2333.
GIESE, A., LOO, M. A., TRAN, N., HASKETT, D., COONS,
S. W. & BERENS, M. E. (1996). Dichotomy of astro-
cytoma migration and proliferation. Int. J. Cancer 67,
275}282.
GREENSPAN, H. P. (1976). On the growth and stability of
cell cultures and solid tumors. J. theor. Biol. 56,
229}242.
HAJI-KARIM, M. & CARLSSON, J. (1978). Proliferation and
viability in cellular spheroids of human origin. Cancer Res.
38, 1457}1464.
HAROSKE, G., DIMMER, V., STEINDORF, D., SCHILLING, U.,
THEISSIG, F. & KUNZE, K. D. (1996). Cellular sociology of
proliferating tumor cells in invasive ductal breast cancer.
Analyt. Quant. Cytol. Histol. 18, 191}198.
HELMLINGER, G., NETTI, P. A., LICHTENBELD, H. C.,
MELDER, R. J. & JAIN, R. K. (1997). Solid stress inhibits the
growth of multicellular tumor spheroids. Nat. Biotechnol.
15, 778}783.
HOSHINO, T. & WILSON, C. B. (1975). Review of basic
concepts of cell kinetics as applied to brain tumors.
J. Neurosurg. 42, 123}131.
HOSHINO, T. & WILSON, C. B. (1979). Cell kinetic analyses
of human malignant brain tumors (gliomas). Cancer 44,
956}962.
KAUFFMAN, S. A. (1984). Emergent properties in random
complex automata. Physica 10D, 145}156.
KELLY, P. J., DUMAS-DAVENPORT, C., KISPERT, D. B.,
KALL, B. A., SCHEITHAUER, B. W. & ILLIG, J. J.
(1987). Imaging-based sterotactic serial biopsies in un-
treated intracranial glial neoplasms. J. Neurosurg. 66,
865}874.
KISS, R., CAMBY, I., SALMON, I., VAN HAM, P., BROTCHI,
J. & PASTEELS, J. (1995). Relationship between DNA
ploidy level and tumor sociology behavior in 12 nervous
cell lines. Cytometry 20, 118}126.
KRAUS, M. & WOLF, B. (1993). Emergence of self-organiza-
tion in tumor cells: relevance for diagnosis and therapy.
umor Biol. 14, 338}353.
LANDRY, J., FREYER, J. P. & SUTHERLAND, R. M. (1981).
Shedding of mitotic cells from the surface of multicell
spheroids during growth. J. Cell. Physiol. 106, 23}32.
LANG, F. F., MILLER, D. C., KOSLOV, M. & NEWCOMB,
E. W. (1994). Pathways leading to glioblastoma-multi-
forme: a molecular analysis of genetic alteration in 65. J.
Neurosurg. 81, 427}436.
MARUSIC, M., BAJZER, Z., FREYER, J. P. & VUC-PAVLOVIC,
S. (1994). Analysis of growth of multicellular tumour sphe-
roids by mathematical models. Cell Prolif. 27, 73}94.
MUCKE, E. (1997). http://www.geom.umn.edu/locate/cglist/
GeomDir.
MUELLER-KLIESER, W. (1997). Three-dimensional cell
cultures: from molecular mechanisms to clinical applica-
tions. Am. J. Physiol. 273, C1109}C1123.
MUELLER-KLIESER, W., FREYER, J. P. & SUTHERLAND,
R. M. (1986). In#uence of glucose and oxygen supply
conditions on the oxygenation of multicellular spheroids.
Br. J. Cancer 53, 345}353.
NAZZARO, J. M. & NEUWELT, E. A. (1990). The role of
surgery in the management of supratentorial intermediate
and high-grade astrocytomas in adults. J. Neurosurg. 73,
331}344.
NORTON, L. (1988). A Gompertzian model of human breast
cancer growth. Cancer Res. 48, 7067}7071.
OKABE, A., BOOTS, B. & SUGIHARA, K. (1992). Spatial
esselations. New York: Wiley.
PAULUS, W. & PFEIFFER, J. (1989). Intratumoral histologic
heterogeneity of gliomas: a quantitative study. Cancer 64,
442}447.
PERTUISET, B., DOUGHERTY, D., CROMEYER, C., HOSHINO,
T., BERGER, M. & ROSENBLUM, M. L. (1985). Stem cell
studies of human malignant brain tumors. Part 2.: Prolifer-
ation kinetics of brain-tumor cells in vitro in early-passage
cultures. J. Neurosurg. 63, 426}432.
PERUMPANANI, A., SHERRATT, J. A., NORBURY, J. & BYRNE,
H. M. (1996). Biological inferences from a mathematical
model for malignant invasion. Invasion Metastasis 16,
209}221.
PIERALLINI, A., BONAMIMI, M., OSTI, M. F., PANTANO, P.,
PALMEGGIANI, F., SANTORO, A., ENRICI, R. M. & BOZZAO,
L. (1996). Supratentorial glioblastoma: neuroradiological
"ndings and survival after surgery and radiotherapy.
Neuroradiology 38, S26}S30.
PRESTON JR., K. & SIDERITS, R. (1992). New techniques for
three-dimensional data analysis in histopathology. Analyt.
Quant. Cytol. Histol. 14, 398}406.
QI, A. S., ZHENG, X., DU, C. Y. & AN, B. S. (1993). A cellular
automaton model of cancerous growth. J. theor. Biol. 161,
1}12.
ROTIN, D., ROBINSON, B. & TANNOCK, I. F. (1986).
In#uence of hypoxia and an acidic environment on the
metabolism and viability of cultured cells: potential
implications for cell death in tumors. Cancer Res. 46,
2821}2826.
RUBIN, H. (1994). Experimental control of neoplastic pro-
gression in cell populations: Foulds' rules revisited. Proc.
Nat. Acad. Sci. ;.S.A. 91, 6619}6623.
SCHIFFER, D., CAVALLA, P., DUTTO, A. & BORSOTTI, L.
(1997). Cell proliferation and invasion in malignant
gliomas. Anticancer Res. 17, 61}70.
TUMOR DYNAMICS VIA CELLULAR AUTOMATON 381
SCHWAB, E. D. & PIENTA, K. J. (1996). Cancer as a complex
adaptive system. Med. Hypotheses 47, 235}241.
SELKER, R. G., MENDELOW, H., WALKER, M., SHEPTAK,
P. E. & PHILIPS, J. G. (1982). Pathological correlation of
CT ring in recurrent, previously treated gliomas. Surg.
Neurol. 17, 251}254.
SHERRATT, J. A. & NOWAK, M. A. (1992). Oncogenes,
anti-oncogenes and the immune response to cancer:
a mathematical model. Proc. Roy. Soc. ondon [B] 248,
261}271.
SILBERGELD, D. L. & CHICOINE, M. R. (1997). Isolation and
characterization of human malignant glioma cells from
histologically normal brain. J. Neurosurg. 86, 525}531.
SMOLLE, J. (1998). Cellular automaton simulation of tumour
growth*equivocal relationships between simulation
parameters and morphologic pattern features. Analyt. Cell.
Pathol. 17, 71}82.
SMOLLE, J. & STETTNER, H. (1993). Computer simulation of
tumour cell invasion by a stochastic growth model.
J. theor. Biol. 160, 63}72.
STEEL, G. G. (1977). Growth Kinetics of umors. Oxford:
Clarendon Press.
SUH, O. & WEISS, L. (1984). The development of a technique
for the morphometric analysis of invasion in cancer.
J. theor. Biol. 107, 547}561.
SUTHERLAND, R. M. (1988). Cell and environment interac-
tions in tumor microregions: the multicell spheroid model.
Science 240, 177}184.
TRACQUI, P. (1995). From passive di!usion to active cellular
migration in mathematical models of tumor invasion. Acta
Biotheor. 43, 443}464.
TRACQUI, P., CRUYWAGEN, G. C., WOODWARD, D. E.,
BARTOO, G. T., MURRAY, J. D. & ALVORD, E. C. (1995).
A mathematical model of glioma growth: the e!ect of
chemotherapy on spatio-temporal growth. Cell Prolif. 28,
17}31.
TURNER, G. A. & WEISS, L. (1980). Some e!ects of products
from necrotic regions of tumours on the in vitro migration
of cancer and peritoneal exudate cells. Int. J. Cancer 26,
247}254.
VAIDYA, V. G. & ALEXANDRO Jr., F. J. (1982). Evaluation
of some mathematical models for tumor growth. Int. J.
Biomed. Comput. 13, 19}35.
WALISZEWSKI, P., MOLSKI, M. & KONARSKI, J. (1998).
On the holistic approach in cellular and cancer biology:
nonlinearity, complexity and quasi-determinism of the
dynamic cellular network. J. Surg. Oncol. 68, 70}78.
WARD, J. P. & KING, J. R. (1997). Mathematical modelling
of avascular-tumour growth. IMA J. Math. Appl. Med.
Biol. 14, 39}69.
WASSERMAN, R., ACHARYA, R., SIBATA, C. & SHIN, K. H.
(1996). A patient-speci"c in vivo tumor model. Math. Bio-
sci. 136, 111}140.
WERNER, M. H., PHUPHANICH, S. & LYMAN, G. H. (1995).
The increasing incidence of malignant gliomas and pri-
mary central nervous system lymphoma in the elderly.
Cancer 76, 1634}1642.
WHITTLE, I. R. (1996). Management of primary malignant
brain tumours. J. Neurol. Neurosurg. Psychiatry 60, 2}5.
WOLFRAM, S. (1984). Cellular automata as models of com-
plexity. Nature 311, 419}424.
WOODWARD, D. E., COOK, J., TRACQUI, P., CRUYWAGEN,
G. C., MURRAY, J. D. & ALVORD, E. C. (1996). A math-
ematical model of glioma growth: the e!ect of extent of
surgical resection. Cell Prolif. 29, 269}328.
YAMASHITA, T. & KUWABARA, T. (1983). Estimation of rate
of growth of malignant brain tumors by computed tomo-
graphy scanning. Surg. Neurol. 20, 464}470.
YOSHII, Y., MAKI, Y., TSUBOI, K., TOMONO, Y., NAKAGAWA,
K. & HOSHINO, T. (1986). Estimation of growth fraction
with bromodeoxyuridine in human central nervous system
tumors. J. Neurosurg. 65, 659}663.
382 A. R. KANSAL E A.