Selection in spatial stochastic models of cancer: Migration as a key modulator of fitness
 Craig J Thalhauser^{1},
 John S Lowengrub^{1},
 Dwayne Stupack^{2} and
 Natalia L Komarova^{1}Email author
DOI: 10.1186/17456150521
© Thalhauser et al; licensee BioMed Central Ltd. 2010
Received: 1 February 2010
Accepted: 20 April 2010
Published: 20 April 2010
Abstract
Background
We study the selection dynamics in a heterogeneous spatial colony of cells. We use two spatial generalizations of the Moran process, which include cell divisions, death and migration. In the first model, migration is included explicitly as movement to a proximal location. In the second, migration is implicit, through the varied ability of cell types to place their offspring a distance away, in response to another cell's death.
Results
In both models, we find that migration has a direct positive impact on the ability of a single mutant cell to invade a preexisting colony. Thus, a decrease in the growth potential can be compensated by an increase in cell migration. We further find that the neutral ridges (the set of all types with the invasion probability equal to that of the host cells) remain invariant under the increase of system size (for large system sizes), thus making the invasion probability a universal characteristic of the cells selection status. We find that repeated instances of large scale celldeath, such as might arise during therapeutic intervention or host response, strongly select for the migratory phenotype.
Conclusions
These models can help explain the many examples in the biological literature, where genes involved in cell's migratory and invasive machinery are also associated with increased cellular fitness, even though there is no known direct effect of these genes on the cellular reproduction. The models can also help to explain how chemotherapy may provide a selection mechanism for highly invasive phenotypes.
Reviewers
This article was reviewed by Marek Kimmel and Glenn Webb.
Background
Cancerous cells within a tumor compete with one another in a fast paced evolutionary system. At the molecular level, mutations are introduced into the tumoral genome; these mutations may be caused by inherited deficiencies, loss of mismatch repair systems, downregulation of the proofreading checkpoints, and chromosomal instabilities. At the cellular level, these mutations introduce changes in phenotype, some profound but many others subtle. It is the emergence of these mutations, as well as epigenetic events, which generates the incredible flexibility and adaptability of the cancer disease state.
While the effect of certain mutations on the cell's phenotype is reasonably well understood (induction of KRas, loss of p53 and/or Rb, overexpression of matrix metalloproteases [1]), it is conceptually difficult to quantify and analyze the level of genetic heterogeneity within a given tumor. It is even more difficult to measure the forces of natural selection in anything other than broad, descriptive terms.
Speaking in broad evolutionary terms, we would like to understand what cellular characteristics make certain cells more fit than others. If a mutant is introduced in a cell colony, what combinations of the mutant characteristics and "background" characteristics make the mutant cells win the evolutionary competition?
The idea that cancer is an evolutionary process has been applied successfully by many computational biologists, as it allows them to use methods of theoretical population biology and ecology [2–8]. Here we focus on two types of phenotypic changes induced by mutations. The first type involves mutations in genes affecting cell proliferation. Activation of some oncogenes, or inactivation of tumor suppressor genes, change the cells' reproductive capacity, and are thought to be early events in the natural history of many cancers [1]. The second type of genetic change influence the cells' ability to migrate/move. Genes of the second type, while commonly associated with metastases, are also affected in primary tumors [9].
These two types of variation are thought to be implicated in malignant transformations for many (if not all) types of solid tumors. How do the two types of change tradeoff to create a mutant which is "fitter" than the background? Questions of this kind are related to the general theory of fitness landscapes, first introduced by [10]. Fitness is viewed as a surface in a multidimentional space, where the dynamics is assumed to be directed toward local fitness maxima. The global maximum corresponds to the evolutionarily stable strategy [11]. In scenarios where fitness of an individual strategy depends on the current composition of the population (frequencydependent fitness), the formalism of fitness generating functions is used [12].
In this paper we focus on a specific aspect of the general problem of fitness landscapes. Namely, we provide a qualitative framework to study the forces of selection acting within a spatially distributed, stochastic colony of cells, which can vary with regards to the two above mentioned characteristics. The models we construct for this purpose are a spatial generalization of a wellknown Moran process, which was first introduced in [13]. This process has been used recently in cancer modeling (see [14–18]). The first spatial (1D) generalization of the Moran process was described in [19], where we considered the process of onehit and twohit mutant fixation. The simplicity of the (generalized) Moran process enabled us to study analytically, as well as numerically, the role of space in the processes of lossoffunction and gainoffunction mutations, see also [20]. In this spatial Moran process, the cells were allowed to divide in response to a death of a neighboring cell on a 1D grid.
In this paper we construct two spatial generalizations of the Moran process with cell migration. The first model includes explicitly the processes of cell division, death and migration. The second model implicitly describes migration through the varied ability of cell types to place their offspring a distance away, in response to another cell's death. The advantage of the simplified model is its analytical tractability. The two main findings are as follows:

In both models, we find that migration has a direct positive impact on the ability of a single mutant cell to invade a preexisting colony. A decreased fitness due to lesser growth potential may be offset by an increase in cell migration.

The neutral ridges (the set of all types with the invasion probability equal to that of the host cells) remain invariant under the increase of system size (for large system sizes), thus making the invasion probability a universal characteristic of the cells selection status.
Our work is somewhat different from a large body of recent literature where spatial cancer dynamics is studied by means of cellular automata or agentbased modeling (see the recent reviews [21–29] and the references therein). Rather than adding on many biological processes and subtleties to our model, we focus on understanding how just two forces, proliferation and migration, tradeoff to influence the overall fitness of cells.
Results
Explicit motility increases the invasion probability
The result that an increased migration potential of mutant cells increases their ability to invade, can be explained intuitively. In the case of small mutant motility, mutant cells tend to concentrate in one region, and the expansion can only occur near the boundary of that region. An increase in mutants' motility increases the degree of mixing in the population, such that mutant cells spread throughout the space. In this case, mutant growth in enhanced as it can occur throughout the bulk of the colony.
Identification of neutral ridges in the explicit model
Figure 2 shows that it is possible to find more than one pair of parameters (k_{ B }, λ) which correspond to the same probability of invasion. In other words, we can see that different strategiesthat is, different parameter sets (λ, k_{ B }) relative to a specific background k_{ A }may have the same invasion probability. It is particularly interesting to investigate these strategies for the probability of invasion (see figure 2, right panel) which corresponds to the case where each cell as equal fitness (see Methods section below, as well as prior work, e.g. [15]).
We will say that all the strategies corresponding to the invasion probability P = 1/N are neutral to (have the same fitness as) the host cells A. We call the sets of all such strategies "neutral ridges". The relationship between the notions of invasion probability and fitness is discussed in more detail below. While we were unable to derive a simple law which relates the values of k_{ B }, λ and k_{ A }corresponding to P = 1/N, some features of their relationship are worth noting. Firstly, there appears to be an asymptote as k_{ B }→ ∞, implying for each background k_{ A }there is a minimum λ necessary for equal fitness strategies. That is, no amount of increased motility can compensate for too low a growth rate. Also, we notice that the distance between points on the P = 1/N lines remains approximately constant throughout the parameter regime (for example, the measured distance between the neutral ridges corresponding to k_{ A }= 1.0 (circles) and k_{ A }= 0.1 (squares), is constant for different values of k_{ B }with the accuracy of 2%). This implies that the law relating the strategy parameters λ and k_{ B }is somewhat independent of the background k_{ A }; the background parameter merely shifts the relationship between λ and k_{ B }vertically in the (k_{ B }, λ) plane.
The role of the division radius in the implicit model
As seen in figure 4, increasing the division radius of all cells increases the ability of a superior mutant to invade the system. This is consistent with the previous result that in a spatial Moran process with ν_{ A }= ν_{ B }= 1, the probability of mutant invasion is smaller than that for a spacefree model [19]. Increasing the division radius brings a spatial model closer to a spacefree model. The spacefree model is recovered as soon as the division radius is sufficiently big such that the entire domain is contained inside a circle of radius ν. This of course implies that there is a limit to the extent increasing the division radius can increase the invasion probability of a mutant phenotype. If the division radius already spans the entire space, increasing the radius further will have no effect on the probability of the mutant cell to divide. A cell with such a radius can be thought of as being governed by a nonspatial, bulk tumor growth model, such as a capacity growth or logistic law.
Interestingly, observe that increasing the radius from 1.0 (nearest neighbors) to 2.0 nearly recapitulates the theoretical results for a wellmixed system (invasion probability 0.33), see figure 4. This result is intriguing because, for a two dimensional grid of size 441 total cells (21 by 21), a division radius of 2.0 represents approximately 2.7% of the total area, and yet the results are nearly indistinguishable from a scenario in which all 441 cells can interact with each other.
Level sets of invasion probability
First of all, we observe that increasing r_{ B }and ν_{ B }(separately or together) leads to an increase in the invasion probability. This in itself is not surprising given that larger values of the growth potential and division radius will increase the mutants' probability to divide. Thus, if a cell is free to choose any values of r_{ B }and ν_{ B }, it will lead to an unrestricted growth of both of these parameters, which corresponds to the cells' climbing up the "hill" in the fitness landscape of figure 5, which corresponds to large r_{ B }and ν_{ B }. In reality though this is hardly possible and there are biological and energetic constraints on how often a cell divides, and how far it can travel upon division. Therefore "allowed" strategy trajectories in the (r_{ B }, ν_{ B }) space can be introduced where some external constraints do not allow a simultaneous large increase in both parameters. For example, if we look at trajectories of the type α_{ r }r_{ B }+ α_{ ν }ν_{ B }= const, where α_{ r }and α_{ν} are some weights. We observe from figure 5 that the background condition of r_{ A }= 2.0 and ν_{ A }= 2.5 is able to resist most invasions from mutant strategies which maximize one trait at the cost of another (high ν with low r or vice versa). Thus, in this case the mixed strategy is more fit than those relying heavily on only one trait.
Identification of neutral ridges for implicit motility model
We will use the parameters of figure 5 and proceed according to the following algorithm. By definition of the division radius, small changes in ν_{ B }may or may not lead to a change in the number of sites within division range. Therefore, we start by selecting a set of division radii ν_{ B }, and use the trace of the line of equal invasion probability to provide initial estimates for a value of r_{ B }for each ν_{ B }. We then perform simulations over a small range of r_{ B }centered at the initial estimate and proceed until we have found a candidate value r_{ E }for which the invasion probability after 10 sets of 10000 simulations. We then confirm that the candidate pair (r_{ E }, ν_{ E }) is in fact a strategy of equal invasion probability by performing the reverse experiment; that is, introducing a mutant cell with strategy (2.0, 2.5) into a background population of strategy (r_{ E }, ν_{ E }).
In each case, the correlation coefficient is greater than 0.99, which confirms the observation that the growth potential and the division radius of the cells of constant invasion probability are inversely proportional to each other.
The right panel of figure 6 replots the same data in terms of the relative growth potential, r_{ B }/r_{ A }, and the relative number of neighbors, n_{ A }/n_{ B }. We observe that the three neutral ridges from the left panel align.
We deduce that the neutral ridges satisfy the relationship r_{ A }n_{ A }= r_{ B }n_{ B }.
Once this quantity is calculated, we "flip a coin", and replace an arbitrarily chosen cell with a cell of the given type. In this process, in order for the cells of type B to have the same invasion probability as cells of type A, the probability of cell division for cells of type A should be equal to the probability of division of cells of type B. That is, r_{ A }n_{ A }= r_{ B }n_{ B }. Resolving this equation, we obtain n_{ B }= (n_{ A }r_{ A })/r_{ B }. This is an inverse dependence with the proportionality coefficient c = n_{ A }r_{ A }.
Invasion probability and fitness
So far we have been investigating the probability of invasion of type B starting from one such mutant. This probability is connected to the relative fitness of type B, which we denote as ℱ. Relative fitness is a parameter that is related to the rate of expansion of a phenotype. It is usually defined as the (averaged) frequency of the type in the "next" generation divided by that in the current generation: . A related concept of invasion fitness is defined as the exponential growth rate of the mutant type in a host population. The connection of the relative fitness parameter (ℱ) to the probability of invasion (ρ) has the following useful properties:
(i) If ρ = 1/N (neutral mutant B) then ℱ = 1 (the number of neutral cells of type B stays constant on average, from generation to generation).
(ii) If ℱ → ∞, that is, if the mutants are strongly advantageous, then ρ → 1.
(iii) The invasion probability and fitness are positively correlated (the types that have a higher invasion probability ρ will have a higher fitness and vise versa).
The notions of fitness and probability of invasion are both important in theoretical biology, and both have advantages and disadvantages. Arguably, the notion of invasion probability is more informative in our setting.
The fitness parameter is defined for a particular temporal dynamics. In our simple model, we use a discretetime Markov chain. This is an unrealistic way to treat the mutant dynamics, but it preserves the invasion probabilities. In other words, if one considers the longterm outcomes of a stochastic process, the dynamics becomes unimportant, and can be simplified in the way implemented here. The expansion rate of cells (equivalent to fitness) is dependent on this simplification and thus is affected by this (unrealistic) aspect of our model.
The probability of invasion does not depend on the exact way we choose the timestep for our updates. Thus in this sense it is a better choice of a competitiveness measure for our system. On the other hand, invasion probability depends on the total number of cells, N. The probability to invade a very small constantsize population is higher than that for a large colony. We have investigated how the probability of invasion depends on N, and found the following.
that is, the invasion probability strongly depends on N. Despite this fact, the notion of invasion probability retains a degree of universality even in the case of neutral mutations. Namely, we found numerically that the expressions for the neutral ridges given by functions r_{ B }n_{ B }= c are Nindependent, that is, the proportionality constant c does not change with N (not shown). This is consistent with the derivation of the parabolic dependence of the level sets for the invasion probability.
Ultimately, we are interested in the probability of invasion as long as it is equivalent to the probability of mutants to thrive. When considering cancer, it is irrelevant whether exactly all the host cells in an organ have been replaced by the mutants (which is equivalent to invasion, rigorously speaking). An important notion is whether a mutant colony expands and persists inside an organ for a long time, which, for large values of N and for advantageous mutants, is very close to the probability of invasion. This notion is more universal than the fitness of mutant cells because it is independent of the timeevolution or the system size. Another meaningful concept is neutrality, which is basically a symmetry property: a neutral mutant behaves just like any host cell, and has the same invasion probability as any other cell (1/N). Its fitness is 1 and its expansion rate (invasion fitness) is zero, which is extremely hard to measure. Instead, measuring the invasion probability gives us a useful tool to identify the class of neutral mutants.
In the rest of this section we will discuss the probability of invasion and make some inferences about the fitness of the mutants. As explained above, the two notions are positively correlated, so all the arguments resulting from considering the level sets of the invasion probability hold for the mutant fitness.
The optimal mutant strategy
We can also safely assume that g' ≤ 0, because of the nature of the constraints we consider; simply speaking, an increase in the growth potential must lead to a somewhat reduced division radius, thus making the function g monotonically increasing. The function n_{ B }g(n_{ B }) can take a variety of shapes. Here we restrict ourselves to the cases where it has at most one internal extremum. Multiple minima and maxima could occur, but they are a result of specific biological facts which we cannot specify at the present level of generality. Thus we focus on the simplest types of generic behavior, keeping in mind that other functions can be studied in a similar way. With at most one internal extremum, there can be the following cases:

The function n_{ B }g(n_{ B }) has an internal maximum. Then the function ℱ also has an internal maximum at the corresponding point (ν_{ B }, r_{ B }). This means that the optimal strategy is to find an intermediate value of growth potential and division radius, instead of maximizing them both.

The function n_{ B }g(n_{ B }) is monotonically decreasing; this could be a result of a constraint that small values of ν are not allowed. Then, the maximum of F is achieved for smallest possible values of n_{ B }and the largest possible values of r_{ B }. In other words, the optimal strategy is maximizing the growth potential (within the allowed range) at the cost of the division radius. This situation could arise in the following scenario: suppose that increasing ν is "expensive", that is, the constraint function f (r_{ B }, ν_{ B }) depends stronger on ν_{ B }than it does on r_{ B }. Then, the value of the derivative, g', is relatively large, which shifts the location of the maximum of the function ν_{ B }g(ν_{ B }) to the left. If the dependence of ν_{ B }is sufficiently strong, then this can drive the location of the maximum outside the allowed domain of ν_{ B }. As a result, the optimal strategy will be to find the smallest possible ν_{ B }.

The function n_{ B }g(n_{ B }) is monotonically increasing in the allowed domain. Then, the maximum of F is achieved for largest possible values of n_{ B }and the smallest possible values of r_{ B }. Analogous to the previous argument, such scenarios could arise when increasing ν_{ B }is "cheaper" than increasing r.
To summarize, we find that the growth potential and the division radius are two components of fitness which play a very similar role in the probability of mutant invasion. The fitness landscape has hyperbolic level sets in terms of the growth potential and the number of neighbors. A mixed strategy is optimal unless one of the two fitness components (the growth potential or the number of neighbors) is much more evolutionary "expensive". In the latter case the less expensive characteristic should be maximized at the expense of the other.
Discussion
A central idea of this paper is that the cell with the fastest intrinsic growth rate is not always the most fit. Rather, it is the culmination of ability to divide with increased opportunity to do so. These forces must be balanced, depending on the relative costs of both adaptations for a cell.
We have demonstrated how relatively simple, twocomponent, models can help explore the complicated phenomenon of phenotypic heterogeneity. With the multiphasic nature of fitness presented here, many different strategies have equal or nearly equal fitness (figures 3 and 6). As the number of traits under consideration increases, it is similarly expected that the potential combinations of strategies which lead to equal fitness will increase. This would reflected in the existence of multidimentional sets equivalent to neutral ridges described here.
The heterogeneous nature of tumors provides a primary mechanism for resistance to current therapeutic strategies. Tumors behave as fast evolutionary systems, with many subtly different phenotypes coexisting, cooperating and/or competing within a shared environment. Changes in phenotypes can arise de novo from altered gene expression profiles, rapid cell proliferation outpacing DNA repair mechanisms, loss of repair systems or checkpoints, alterations to chromosomal integrity and epigenetic events. The application of therapy (chemical, radiation) to this system changes the selective forces acting upon the various phenotypes, permitting rigorous selection and socalled "punctuated evolution" within the tumor [30]. These external forces can act to select and advantage a phenotype that might otherwise be unlikely to achieve dominance, or to accelerate the rise of a given phenotype within a population.
Fitness ridges, fitness landscapes and chemotherapy
Assessing the effectiveness of a therapeutic intervention, via mass dieoff and selection, underlies the difficulty in studying, understanding and treating cancer. The plasticity of the disease will alter parameters as the disease is treated. Clearly, a predictive understanding of how mutants penetrate both pretreatment and posttreatment tumors, will be critical for designing new therapy strategies.
Our model can help understand how therapy influences the fitness landscape of the tumor. We test the hypothesis that a higher death rate generates a positive selective pressure on a more mobile phenotype. Such a death rate might describe a form of medical intervention. To test this hypothesis, we allow significantly more than one cell to die per iteration of our process. We begin with a system composed of 10% mutant cells in a background of wild type. At each step of the iteration, each cell in the space has an independent probability to die of P_{ d }. Thus, after the mass death there will be on average (1  P_{ d }) * N cells remaining in the system. The process described previously is then simulated without any new cell death until the space is filled again. At the end of each such iteration we determine the number of type A and type B cells in the system. This process is repeated for a fixed number of mass death events. We also determine whether or not one cell type has completely invaded the system.
Thus, we have demonstrated that application of a toxic therapeutic intervention that results in gross tumor death can significantly promote the invasion of a more motile mutant into the general tumor in fact, our data suggest that this advantages the invasive capacity of a motile tumor relative to the faster proliferating tumor. This somewhat unexpected observation may reflect the clinical situation; many malignant tumors exhibit enhanced expression of genes that promote motility, as discussed below.
In the implicit migration model considered here we found that the the fitness of cells depends on two parameters, the growth potential and division radius. Moreover, we demonstrated that it depends on the product of the two. This is not an unusual finding in evolutionary biology. Another example comes from virus dynamics, where the fitness of viruses (measured in terms of their basic reproductive ratio) can be calculated as a function of various parameters in predatorprey type models. It is typically proportional to the product of two quantities: the infectivity of the virus and the inverse of the death rate of infected cells. The former parameter can be roughly related to the "motility" of the virus, or to how efficiently it infects new cells, and in this sense it is roughly analogous to our division radius. The inverse death rate of infected cells is reflective of our growth potential since a decrease in the cell's lifespan (a decrease in the inverse death rate) leads to a decrease in the total number of progeny.
Several papers have been devoted to studying optimal viral strategies, especially given that the two components of fitness are not independent [31–38]. In spirit, our approach is similar to that. While there has to date been little investigation into this area, it seems reasonable to assume energy availability limits the total growth and migration ability of an invasive cancer cell; that is, since both growth and migration are energyintensive processes, an invasive cell will have to 'budget' its energy to each.
Biological Applications
The model presented in this paper is a very crude approximation of reality. Thus the prediction generated by the model may or may not be relevant for biological applications. However, it is particularly compelling that several genetic alterations have been discovered which are implicated in promoting migratory and/or invasive phenotypes, yet play no direct role in regulating cell division. The in vivo selection for these phenotypes would appear to indicate a strong selection for the fitness conferred. For example, genes that regulate cell interaction with the local extracellular milieu, [39] or those that promote cytoskeleton dynamics [40] have been established to promote cell migration.
The class of small GTPases has been generally divided into families, based on homology and role played in the cell. The Ras family members kras, nras and hras have been linked to cell proliferation, while Rho family members such as Rac, RhoA and RhoC and Rac are small GTP binding proteins the influence cytoskeletal remodeling [41]. These diffierent family members exhibit significant crosstalk, with Ras influencing motility somewhat and Rho proteins sometimes contributing to proliferation. The capacity of numerous cellular proteins found to be elevated in cancer to influence both proliferation and motility, and include a swath of oncogenic proteins such as Src family kinases [42] mitogen activated family kinases [43], and phosphoinnositide 3 kinases [44].
Nonetheless, cytoskeletal effectors which are tightly linked to migration do promote cell invasion and tumor malignancy. Effectors such as RhoC are not oncogenes per se, as they do not transform primary cells, but nonetheless influence cell migration and tumorigeneicity [45, 46]. Moreover, other mediators in the cytoskeleton control system have been shown to impact migration potential and tumorigenicity. The protein ABI1 is an adaptor molecule which, upon RhoA activation, facilitates the formation of multiprotein complexes involved in lamellipodia formation [47]. Overexpression of ABI1 correlates well with migratory and invasion potential in breast cancer cell lines, while suppression of this protein led to the complete loss of migratory ability in once highly invasive cell lines [47].
Coupled to the activity of the cytoskeleton are the cell surface receptors mediating extracellular matrix interaction. Integrins, and their associated proteins, have been shown to be critical for motility and survival. For example, the integrin associated protein focal adhesion kinase (FAK) is a critical determinant of cell migration, yet similar to the small GTPases of the Rho family, the expression of FAK is insufficient to transform primary cells [42, 48]. Moreover, FAK expression is not essential for tumorigenesis, yet FAK expression promotes tumorigenicity and invasion. Accordingly, FAK is found to be enriched in aggressive human tumors. The common theme emerging from these data is that, increased migratory capacity results in a tumor cell with increased fitness. It is not yet clear whether this fitness is sufficient to result in a phenotypic conversion of a tumor in vivo, and it would therefore be interesting to determine whether this occurs in a preclinical model of tumor development. However, our simplified model provides very similar results to those seen in vivo with respect to the importance of migration, and provides an interesting perspective on tumor cell fitness during tumor progression. While the opportunity to divide may often be related to factors that govern the intrinsic proliferative potential of the individual cells, our data suggest that migration can compensate for a lower proliferative rate, and that migration may therefore be an important selection factor during tumor progression in patients. Strikingly, the study also points out the potential for chemotherapy to influence tumor phenotype in an unintended manner.
Conclusions
We have formulated two spatial extensions of the Moran process to study the effects of motility on the ability of a new mutant cell to invade a preexisting tumor. Our results indicate that a mobility phenotype can be dominant and able to invade a background of cells with higher reproduction rates. Further, we have investigated how differing combinations of replication and motility strength blend to form equivalent fitnesses within the tumor system. Finally, we have shown how outside influence on the system inducing mass death within the tumor might unintentionally favor invasion and regrowth by a more aggressively motility phenotype.
Methods
Previous work: the Moran process
The conventional Moran process is formulated as follows. In a population of N cells, each cell is equipped with a nonnegative replication rate. The process is a sequence of updates. At each timestep, one cell is chosen randomly for death and is then replaced by a progeny of another cell (note that in this process, the probability to be picked for death is the same for all cells, regardless of their type). To choose which cell reproduces, one weighs in all cells replication parameters, such that the probability for cell i to reproduce is given by , where r_{ k }is the replication rate of cell k, and the summation is performed over all cells in the population.
In we assume that there are two distinct types in the population, type A with replication parameter r_{ A }and type B with replication parameter r_{ B }, then the stochastic process can be formulated in terms of only one independent random variable, the number of cells of type B. If there are i cells of type B, then after one update the following transitions are possible:

i → i + 1 with probability P_{i → i+1 }= (N  i)/N × r_{ B }i/(r_{ A }(N  i) + r_{ B }i), where the first term is the probability that a cell of type A dies and the second term is the probability for a cell of type B to divide;

i → i  1 with probability P_{i → i1 }= i/N × r_{ A }(N  i)/(r_{ A }(N  i) + r_{ B }i), where the first term is the probability that a cell of type B dies and the second term is the probability for a cell of type A to divide;

i → i with probability 1  P_{i → i+1 } P_{i → i1}.
In the special case of neutral mutants, r_{ A }= r_{ B }, we have ρ = 1/N. This can be obtained from formula (2) by taking the limit r_{ B }→ r_{ A }. Also, this result follows from symmetry considerations: a cell of type B if type B is neutral has the same expansion properties as any of the host cells, and the same probability to invade. Since inevitable one of the N cells will invade, the probability of invasion is 1/N for every cell, including the B cell.
First spatial generalization of the Moran process
In [19] we introduced and analyzed a first spatial generalization of the Moran process. We considered a 1D space, where all the N cells were placed on a regular grid, at locations 1, 2,, N. As before, we assume that the total number of cells does not change. Cells are randomly chosen for death. Each cell death is followed by a cell division of one of its two neighboring cells, which places its daughter cell at the empty slot. Cell death occurs randomly and division is proportional to the relative fitness of the cells.
where r ≡ r_{ B }/r_{ A }. It can be shown that ≤ ρ, with the equality corresponding to neutral mutants, r = 1.
In other words, for nonneutral mutants, the probability to invade is smaller in a spatial model than it is in the spacefree Moran process.
Note that both the spacefree Moran process and the spatial generalization described above depend only on the ratio of the replication parameters r_{ B }/r_{ A }.
The explicit motility model
Our aim in this paper is to study selection processes acting upon cells in spatial settings. The first step toward that goal was made by creating a straightforward spatial generalization of the Moran process described above. However, this model is too simplistic. For example, starting from one mutant cell, a mutant colony can only spread as a solid, connected spot. This is a consequence of the facts that (i) the model is onedimensional, (ii) cells only interact with their nearest neighbors and (iii) cells are not allowed to migrate.
In this paper we introduce two spatial modifications of the Moran process which are higherdimensional, and include cellular motility. This section describes a model in which cell motility is explicitly incorporated. A simpler model, in which motility implicitly described, is discussed in the next section.
Spatial Moran algorithm
 1.
a background cell divides;
 2.
an invading cell divides;
 3.
a background cell migrates;
 4.
an invading cell migrates.
where we denoted K = n_{ A }(r_{ A }+ m_{ A }) + n_{ B }(r_{ B }+ m_{ B }).
In our simulations, we use a square grid with periodic boundary conditions. At time 0 we introduce a single type B mutant at a random position in the space grid which is otherwise filled with type A cells. The process described above is repeated until one species, A or B, has been eliminated from the environment.
Model parameters: the growth and migration potentials
where is the ratio of migration over division potentials, is the ratio of invader to background growth potentials, and = n_{ A }(1 + k_{ A }) + n_{ B }λ (1 + k_{ B }).
Note that increasing the migration potential in this model is similar to decreasing the typical timescale of migration compared to that of cell division. It amounts to an increase in the average distance traveled by a cell in each iteration. An assumption of the model is that the death rate of cells is low compared to the division rate. In other words, we typically do not expect to have multiple cell deaths within a migration radius of a given cell. This is why we can complete each iteration (from a cell death event to a division event) independently.
The implicit motility model
The advantage of the model described in the previous section is the explicit representation of motility within the cell population. A disadvantage is that the relationships for determining equality amongst strategies were found to be too complex to study analytically. Thus, we designed a simpler model of motility, in which cell movement is implicitly tied to the survival of new cell progeny. In this section we present the formulation and the analysis of this simpler model.
Model Scheme
Following the standard Moran process, at each timestep, one cell is selected randomly to die, to be immediately replaced by the progeny of one of its neighboring cells. As before, the difference between the two phenotypes A and B is reflected in their ability to reproduce. Suppose that cells of type A have growth potential r_{ A }and division radius ν_{ A }, and cells of type B have growth potential r_{ B }and division radius ν_{ B }(the distance over which a cell can place its progeny); the biological meaning of these parameters is discussed in the next subsection. The probability for the empty slot to be filled by either the background (A) type or mutant (B) type depends on the division radii and growth potential of each.
The probability for the slot to be filled by a progeny of a cell of type A is P_{ B }= 1  P_{ A }.
Note that the first spatial Moran process described above is a special case of the implicit model which corresponds to a 1D space and ν_{ A }= ν_{ B }= 1.
Model parameters: the growth potential and the division radius
In the process described above, each cell is characterized by two parameters, r_{ A }(or r_{ B }) and ν_{ A }(or ν_{ B }). We call the former the growth potential and the latter division radius. Growth potential is a measure of, given that an empty space has opened sufficiently close to a cell, how likely it is that that cell will divide and fill the space. The corresponding parameter, r, can take any positive values. It is not in any sense a probability to divide, it instead measures how often the cell can divide compared to other phenotypes in the colony. Thus, the results will only depend on the ratios of growth potentials.
The division radius is a measure of what "sufficiently close" in the prior definition represents. From a technical standpoint, a division radius of 1.0 in a 2D square grid corresponds to nearest neighbor interactions, while a radius of 1.5 corresponds to next nearest neighbors.
This implicit model only tracks net division and migration events. Further, the model depends only on two parameters: the relative growth potential r_{ B }/r_{ A }and the relative division radius ν_{ B }/ν_{ A }(this follows from equation (4). Thus, it is a simpler to analyze than the previous model that accounts for migration explicitly. Note that as currently implemented, all cells within the division radius have equal growth potential. It is straightforward to generalize the model such that the growth potential depends on the distance from the empty space.
It is important to note here that the implicit motility model is not intended to be a mechanistic model of growth and evolution of a spatially defined tumor. It is instead an abstraction of a more mechanistic system, the explicit motility model, made under a simplifying assumption, which allows a more analytic treatment of the results. The fundamental assumption allowing us to construct the implicit motility model is that the timescale of migration is much faster than the timescale of tumor cell death in a fully formed tumor. Thus, the division radius might be considered to be the average distance a daughter cell migrates away from a parent cell before the next death event occurs.
Reviewers' Comments
Glenn Webb: There is a continuum modeling approach related to the issues raised in this paper given in the paper. Transforming growth factor TGFbeta is known to inhibit cell proliferation, and also increase cell motility and decrease cellcell adhesion [49]. A version of the FisherKolmogorov equation was used to quantify the simultaneous effects of TGFbeta to increase the tendency of individual cells and cell clusters to move randomly and to decrease overall population growth. Accompanying experiments demonstrated that TGFbeta increased the percentage of mobile cells in an in vitro cell population in a dosedependent manner, consistent with model simulations. Have the authors tried to develop continuum models to quantify migratory impact on invasive cell population dynamics?
Marek Kimmel: The paper introduces very interesting spatial extensions of the Moran model with resulting conclusions concerning invasion by a mutant depending on cell motility and division radius. However, mathematical and simulation results provided, concern only aggregate characteristics such as invasion probability, without discussing the spatial patterns emerging. These latter might be interesting as demonstrated by a class of spatial models, which explore the consequences of linking the model of spatial growth of precancerous cells with diffusion of the growth factors [50]. The picture emerging from modeling indicates that production of growth factors by cells may lead to diffusiondriven instability, which in turn may lead either to decay of both population, or to emergence of local growth foci represented by spikelike solutions. Interesting dependencies arise when two mutualistic cell populations are considered. It would be probably interesting to check if any of the models in the current paper may lead to similar dynamics?
Declarations
Acknowledgements
This work was supported by the NIH/NSF joint initiative on Mathematical Biology through NIH grant GM75309 and NIH/NIGMS grant 1P50GM07651601A1. N.K. also gratefully acknowledges support from a Sloan Foundation Fellowship. J.L. also gratefully acknowledges partial support from the National Science Foundation Division of Mathematical Sciences. D.S. also gratefully support from the NCI through grant NCI R01CA107263. As a onetime exception to the publishing policy of Biology Direct the articles in this series are being published with two reviewers.
Authors’ Affiliations
References
 Hanahan D, Weinberg R: The hallmarks of cancer. CELL. 2000, 100 (1): 5770. 10.1016/S00928674(00)816839.PubMedView ArticleGoogle Scholar
 Nowell P: The clonal evolution of tumor cell populations. Science. 1976, 194: 2328. 10.1126/science.959840.PubMedView ArticleGoogle Scholar
 Breivik J, Gaudernack G: Carcinogenesis and natural selection: a new perspective to the genetics and epigenetics of colorectal cancer. Adv Cancer Res. 1999, 76: 187212. full_text.PubMedView ArticleGoogle Scholar
 Gatenby R, Maini P: Mathematical oncology: cancer summed up. Nature. 2003, 421: 32110.1038/421321a.PubMedView ArticleGoogle Scholar
 Nowak M, Sigmund K: Evolutionary dynamics of biological games. Science. 2004, 303: 793799. 10.1126/science.1093411.PubMedView ArticleGoogle Scholar
 Wodarz D, Komarova N: Computational biology of cancer: lecture notes and mathematical modeling. 2005, World ScientificView ArticleGoogle Scholar
 Vineis P, Berwick M: The population dynamics of cancer: a Darwinian perspective. Int J Epidemiol. 2006, 35: 11511159. 10.1093/ije/dyl185.PubMedView ArticleGoogle Scholar
 Merlo L, Pepper J, Reid B, Maley C: Cancer as an evolutionary and ecological process. Nat Rev Cancer. 2006, 6: 924935. 10.1038/nrc2013.PubMedView ArticleGoogle Scholar
 RinkerSchaeffer CW, O'Keefe JP, Welch DR, Theodorescu D: Metastasis suppressor proteins: Discovery, molecular mechanisms, and clinical application. Clinical cancer research. 2006, 12 (13): 38823889. 10.1158/10780432.CCR061014.PubMedPubMed CentralView ArticleGoogle Scholar
 Wright S: The roles of mutation, inbreeding, crossbreeding, and selection in evolution. Proceedings of the Sixth International Congress on Genetics. 1932, 355366.Google Scholar
 Smith JM: Evolution and the Theory of Games. 1982, Cambridge University Press. CambridgeView ArticleGoogle Scholar
 Vincent TL, Brown JS: Evolutionary game theory, natural selection, and Darwinian dynamics. 2005, Cambridge University Press. NYView ArticleGoogle Scholar
 Moran P: The Statistical Processes of Evolutionary Theory. 1962, Oxford: ClarendonGoogle Scholar
 Nowak MA, Komarova NL, Sengupta A, Jallepalli PV, Shih IM, Vogelstein B, Lengauer C: The role of chromosomal instability in tumor initiation. Proc Natl Acad Sci USA. 2002, 99 (25): 1622616231. 10.1073/pnas.202617399.PubMedPubMed CentralView ArticleGoogle Scholar
 Komarova NL, Sengupta A, Nowak MA: Mutationselection networks of cancer initiation: tumor suppressor genes and chromosomal instability. J Theor Biol. 2003, 223 (4): 433450. 10.1016/S00225193(03)001206.PubMedView ArticleGoogle Scholar
 Nowak MA, Michor F, Komarova NL, Iwasa Y: Evolutionary dynamics of tumor suppressor gene inactivation. Proc Natl Acad Sci USA. 2004, 101 (29): 1063510638. 10.1073/pnas.0400747101.PubMedPubMed CentralView ArticleGoogle Scholar
 Michor F, Iwasa Y, Rajagopalan H, Lengauer C, Nowak MA: Linear model of colon cancer initiation. Cell Cycle. 2004, 3 (3): 358362.PubMedView ArticleGoogle Scholar
 Iwasa Y, Michor F, Nowak MA: Stochastic tunnels in evolutionary dynamics. Genetics. 2004, 166 (3): 15711579. 10.1534/genetics.166.3.1571.PubMedPubMed CentralView ArticleGoogle Scholar
 Komarova NL: Spatial stochastic models for cancer initiation and progression. Bull Math Biol. 2006, 68 (7): 15731599. 10.1007/s1153800590468.PubMedView ArticleGoogle Scholar
 Komarova N: Loss and gainoffunction mutations in cancer: massaction, spatial and hierarchical models. Jour Stat Phys. 2007, 128: 413446. 10.1007/s1095500692380.View ArticleGoogle Scholar
 Deutsch A, Dormann S: Cellular automaton modeling of biological pattern formation. 2005, Birkhauser. BostonGoogle Scholar
 Byrne H, Alarcón T, Owen M, Webb S, Maini P: Modeling Aspects of Cancer Dynamics: A Review. Phi Trans R Soc A. 2006, 364: 15631578. 10.1098/rsta.2006.1786.View ArticleGoogle Scholar
 Fasano A, Bertuzzi A, Gandolfi A: Complex systems in biomedicine. Milan: Springer 2006 chap. Mathematical modelling of tumour growth and treatment, 71108.
 Galle J, Aust G, Schaller G, Beyer T, Drasdo D: Individual cellbased models of the spatial temporal organization of multicellular systemsAchievements and limitations. Cytometry. 2006, 69A: 704710. 10.1002/cyto.a.20287.View ArticleGoogle Scholar
 Drasdo D, Höhme S: On the role of physics in the growth and pattern of multicellular systems: What we learn from individualcell based models?. J Stat Phys. 2007, 128: 287345. 10.1007/s109550079289x.View ArticleGoogle Scholar
 Anderson A, Chaplain M, Rejniak K, Fozard J: Singlecell based models in biology and medicine. Math Med Biol. 2008, 25: 185186. 10.1093/imammb/dqn008.View ArticleGoogle Scholar
 Deisboeck T, Zhang L, Yoon J, Costa J: In silico cancer modeling: is it ready for prime time?. 2008, 6: 3442.Google Scholar
 Anderson A, Quaranta V: Integrative mathematical oncology. Nature Reviews Cancer. 2008, 8: 227244. 10.1038/nrc2329.PubMedView ArticleGoogle Scholar
 Quaranta V, Rejniak K, Gerlee P, Anderson A: Invasion emerges from cancer cell adaptation to competitive microenvironments: Quantitative predictions from multiscale mathematical models. Sem Cancer Biol. 2008, 18: 338348. 10.1016/j.semcancer.2008.03.018.View ArticleGoogle Scholar
 Gould S: Tempo and mode in the macroevolutionary reconstruction of Darwinism. Proc Nat Acad Sci USA. 1994, 91: 67646771. 10.1073/pnas.91.15.6764.PubMedPubMed CentralView ArticleGoogle Scholar
 Anderson RM, May RM: Coevolution of hosts and parasites. Parasitology. 1982, 85 (Pt 2): 411426. 10.1017/S0031182000055360.PubMedView ArticleGoogle Scholar
 Nowak MA, May RM: Superinfection and the evolution of parasite virulence. Proc Biol Sci. 1994, 255 (1342): 8189. 10.1098/rspb.1994.0012.PubMedView ArticleGoogle Scholar
 Levin BR: The evolution and maintenance of virulence in microparasites. Emerg Infect Dis. 1996, 2 (2): 93102. 10.3201/eid0202.960203.PubMedPubMed CentralView ArticleGoogle Scholar
 Ebert D, Herre EA: The evolution of parasitic diseases. Parasitol Today. 1996, 12 (3): 96101. 10.1016/01694758(96)806685.PubMedView ArticleGoogle Scholar
 Frank SA: Models of parasite virulence. Q Rev Biol. 1996, 71: 3778. 10.1086/419267.PubMedView ArticleGoogle Scholar
 Ebert D, Mangin KL: The influence of host demography on the evolution of virulence of a microsporidian gut parasite. Evolution. 1997, 51: 18281837. 10.2307/2411005.View ArticleGoogle Scholar
 Sasaki A, Boots M: Parasite evolution and extinctions. Ecology Letters. 2003, 6 (3): 17610.1046/j.14610248.2003.00426.x.View ArticleGoogle Scholar
 Boots M, Hudson PJ, Sasaki A: Large shifts in pathogen virulence relate to host population structure. Science. 2004, 303 (5659): 842844. 10.1126/science.1088542.PubMedView ArticleGoogle Scholar
 Streuli C, Akhtar N: Signal cooperation between integrins and other receptor systems. Biochem J. 2009, 418: 491506. 10.1042/BJ20081948.PubMedView ArticleGoogle Scholar
 Vogel V, Sheetz M: Cell fate regulation by coupling mechanical cycles to biochemical signalling pathways. Biochem J. 2009, 418: 491506. 10.1042/BJ20081948.View ArticleGoogle Scholar
 Bustelo X, Sauzeau V, Berenjeno I: GTPbinding proteins of the Rho/Rac family: regulation, effectors and functions in vivo. Bioessays. 2007, 29: 356370. 10.1002/bies.20558.PubMedPubMed CentralView ArticleGoogle Scholar
 Burnton V, Frame M: Src and focal adhesion kinase as therapeutic targets in cancer. Curr Opin Pharmacol. 2008, 8: 427432. 10.1016/j.coph.2008.06.012.View ArticleGoogle Scholar
 Boutros T, Chevet E, Metrakos P: Mitogenactivated protein (MAP) kinase/MAP kinase phosphatase regulation: roles in cell growth, death and cancer. Pharmacol Rev. 2008, 60: 261310. 10.1124/pr.107.00106.PubMedView ArticleGoogle Scholar
 Yuan T, Cantley L: PI3K pathway alterations in cancer: variations on a theme. Oncogene. 2008, 27: 54975510. 10.1038/onc.2008.245.PubMedPubMed CentralView ArticleGoogle Scholar
 Clark E, Golub T, Lander E, Hynes R: Genomic analysis of metastasis reveals an essential role for RhoC. Nature. 2000, 406: 532535. 10.1038/35020106.PubMedView ArticleGoogle Scholar
 Stoletov K, Motel V, Lester R, Gonias S, Klemke R: Highresolution imaging of the dynamic tumor cellvascular interface in transparent zebrafish. Nature. 2007, 104: 1740617411.Google Scholar
 Larkins TL, Nowell M, Singh S, Sanford GL: Inhibition of cyclooxygenase2 decreases breast cancer cell motility, invasion and matrix metalloproteinase expression. BMC CANCER. 2006, 6: 10.1186/147124076181.Google Scholar
 Lim S, Mikolon D, Stupack D, Schlaepfer D: FERM control of FAK function: Implications for cancer therapy. Cell Cycle. 2008, 7: 23062314.PubMedPubMed CentralView ArticleGoogle Scholar
 Wang SE, Hinow P, Bryce N, Weaver AM, Estrada L, Arteaga CL, Webb GF: A mathematical model quantifies proliferation and motility effects of TGF on cancer cells. Computational and mathematical methods in medicine. 2009, 10 (1): 7183. 10.1080/17486700802171993.PubMedPubMed CentralView ArticleGoogle Scholar
 MarciniakCzochra A, Kimmel M: Modelling of early lung cancer progression: Influence of growth factor production and cooperation between partially transformed cells. Mathmatical models and methods in applied sciences. 2007, 17 (1, Suppl S): 16931719. 10.1142/S0218202507002443.View ArticleGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.