Spatial measures of genetic heterogeneity during carcinogenesis
SSpatial measures of genetic heterogeneity duringcarcinogenesis
Kathleen Storey , Marc D. Ryser , Kevin Leder , and Jasmine Foo School of Mathematics, University of Minnesota, Minneapolis, MN, USA Department of Mathematics, Duke University, Durham, NC, USA Industrial and Systems Engineering, Minneapolis, MN, USA, [email protected] School of Mathematics, Minneapolis, MN, USA, [email protected]
October 21, 2018
Abstract
In this work we explore the temporal dynamics of spatial heterogeneityduring the process of tumorigenesis from healthy tissue. We utilize a spatialstochastic process model of mutation accumulation and clonal expansion ina structured tissue to describe this process. Under a two-step tumorigenesismodel, we first derive estimates of a non-spatial measure of diversity: Simp-son’s Index, which is the probability that two individuals sampled at randomfrom the population are identical, in the premalignant population. We nextanalyze two new measures of spatial population heterogeneity. In particularwe study the typical length scale of genetic heterogeneity during the carcino-genesis process and estimate the extent of a surrounding premalignant clonegiven a clinical observation of a premalignant point biopsy. This evolutionaryframework contributes to a growing literature focused on developing a betterunderstanding of the spatial population dynamics of cancer initiation and pro-gression. Although initially motivated by understanding questions in cancer,these results can be applied more generally to help understand the dynam-ics of heterogeneity and diversity in a variety of spatially structured, evolvingpopulations.
Carcinogenesis, the transformation from healthy tissue to invasive cancer, is alengthy and complex process driven by a variety of factors including hereditarypredisposition [1], exposure to environmental factors [2] and a changing microenvi-ronment in the affected organ [3]. Irrespective of the driving factors, most cancersare characterized by the progressive accumulation of genetic alterations in a smallgroup of founder cells. These alterations are either deleterious or neutral (passengermutation), and some can confer a fitness advantage to the affected cell (driver mu-tation) by increasing the reproductive rate or inhibiting cell-regulatory mechanisms1 a r X i v : . [ q - b i o . P E ] O c t α -, β - and γ -diversity to denote the averagespecies richness at the single habitat level ( α ), the diversity between habitats ( β ),and total species richness ( γ ) [10]. These measures are useful to quantify largescale organismal diversity in an ecological setting with spatial variation betweenwell-defined habitats. However in the present work we are interested in developing2ew measures of diversity to explore more specifically the intrinsic length scales ofgenetic heterogeneity driven by clonal expansion dynamics in a spatially structuredtissue population.There have been other mathematical modeling efforts on the topic of heterogene-ity during cancer initiation and expansion. In particular, previous work by Iwasaand Michor explored the Simpson’s Index in a Moran process of tumorigenesis [11].This study focused on understanding the impact of neutral and advantageous mu-tations in a non-spatial, homogeneously mixed population setting. The work byDurrett and et. al. [12] developed formulas for Simpson’s Index and other hetero-geneity measures in a multitype branching process model of cancer evolution. Morerecently, Dhawan and colleagues [13] developed a computational platform for thecomparison of alternative spatial heterogeneity measures as potential biomarkers fortumor progression. Finally, within the broader context of spatial tumor growth, ourwork adds to a vast body of literature, e.g. [14, 15, 16, 17, 18, 7].The outline of this paper is as follows: In Section 2 we introduce a cell-basedstochastic evolutionary model of spatial carcinogenesis, as well as a mesoscopic ap-proximation to this model that was analyzed in [8]. In Section 3 we first analyzethe non-spatial Simpson’s Index for this spatially-structured population. Then, inSection 4 we formulate and analyze two clinically relevant spatial measures of het-erogeneity and study their dependence on cancer-specific parameters. Finally, wesummarize and discuss our findings in Section 5. We introduce a spatial evolutionary model that describes the dynamic transitionfrom physiological homeostasis to onset of invasive cancer. In between, the tissueundergoes a sequence of genetic changes that manifest themselves at the phenotypiclevel in the form of increased proliferation rates, and hence a fitness advantage ofmutant cells over normal cells. It is important to note that in many cancers, thereis a succinct lack of a clearly defined genetic sequence [19]. On the other hand,the morphological changes from normal tissue to dysplasia, carcinoma in situ andinvasive cancer is common in carcinomas, which account for over 80% of all cancers.Therefore, one might prefer to interpret mutations as phenotypic transitions ratherthan genetic aberrations. With this interpretation in mind, we are going to introducea linear 3-stage model, where type 0 cells represent normal tissue, type 1 cells arepre-malignant (dysplasia/CIS) and type 2 cells are malignant cancer cells. Notethat the model can be extended to a setting with more than two mutations, eitherto represent a more refined phenotypic progression, or to account for select cancer-specific genetic events.To render this model spatial, we introduce a cell-based stochastic model on theinteger lattice Z d ∩ [ − L/ , L/ d , where L >
0, and equip this domain with periodicboundary conditions. On this lattice we have three different types of cells, labeledas type 0, type 1 and type 2. For i ∈ { , , } a type i cell reproduces at rate (1 + s ) i ,and when the cell reproduces it replaces one of its 2 d neighboring cells at random.3n addition, we assume that for i ∈ { , } , a type i cell mutates to type i + 1 atrate u i +1 . Initially our entire lattice is occupied by type 0 cells which representnormal cells without any oncogenic mutations. Tumor initiation is defined as thebirth of the first type-2 cell that does not go extinct. In the biological applicationwe are interested in (somatic cells in the body) L is generally at least 10 while s , u and u are quite small. Therefore we will, unless stated otherwise, restrict ouranalysis to the regime L (cid:28) u (cid:28) u (cid:28)
1, and s (cid:28)
1. Before we can discussthe specific conditions imposed on the model parameters, we need to review thedynamic properties of the model.In [8] we established that the arrival of type-1 mutants that are successful (i.e.whose progeny does not go extinct) can be described as a Poisson process with rate u s/ (1 + s ). Here, u is the mutation rate to type-1 and s/ (1 + s ) is the survivalprobability of each type-1 mutant. We also characterized the radial expansion rateof type-1 families as a function of the selective advantage s in each dimension. Inparticular, it was established in [20, 21] that each successful type-1 family has anasymptotic shape D which is a unit ball in an unknown norm, and grows linearly intime. Let e the unit vector in the first axis, and let c d ( s ) be the linear expansion rateof the radius of this ball: D ∩ { ze : z ∈ R } = [ − c d ( s ) , c d ( s )] . Then, we establishedthat as s → c d ( s ) ∼ s d = 1 (cid:112) πs/ log(1 /s ) d = 2 √ β d s d ≥ , where β d is the probability that two d dimensional simple random walks started at0 and e = (1 , , . . . ,
0) never collide.Based on this result, we then introduced a mesoscopic approximation to themodel. Here, the growth of successful mutant families is deterministic while thearrival of these families follows a non-homogeneous Poisson process. To ensure thatthis mesocopic model accurately recapitulates the dynamics of the cell-based model,we will make the following assumptions on the relationships between parameters inthe model: ( A u (cid:28) /(cid:96) ( s ) ( d +2) / (1)( A (cid:18) c d u s (cid:19) d/ ( d +1) (cid:28) N ( A
2) (
N u s ) d +1 ( c dd u s ) − → c ∈ [0 , ∞ ) (2)( A u (cid:28) /(cid:96) ( s )These assumptions generally hold for the parameter ranges appropriate for our bio-logical application of carcinogenesis, see [8] for a details. In addition, we will focuson dimensions d = 1 and d = 2 since most epithelial tissues can be viewed as oneor two dimensional structures, e.g., the cells lining a mammary duct ( d = 1), the4rypts in the colon ( d = 2), or the stratified squamous epithelia of bladder, thecervix and the skin ( d = 2).In the simplified mesoscopic model we consider the cells to live on a spatialcontinuum D = [ − L, L ] d . The state-space of the system is given by a set-valuedfunction χ t , which characterizes the regions of D occupied by type-1 cells at time t .Mutations to type-1 cells occur as a Poisson process at rate u s in the set χ ct = D \ χ t ,i.e. in regions where type-0 cells reside. Each newly created type-1 mutation initiatesan expanding ball whose radius grows linearly at rate c d . Denoting the Euclideanball by B x,r = { y : || y − x || ≤ r } , then after k mutations at the space-time points { ( x , t ) , . . . , ( x k , t k ) } , we have χ t = k (cid:91) i =1 B x i ,c d ( t − t i ) . Thus, the state of the system at any time t is the union of balls occupied by ex-panding mutant type-1 families. In [8] we proved that under assumption (A2) wecan neglect the possibility that a second mutation arises from a type-1 family thatdies out eventually. Therefore, we model successful type-2 mutations as Poissonarrivals into the space occupied by type-1 cells, χ t , with rate u s . Recall that in ourtwo-step cancer initiation model, the type-2 mutant represents a malignant cancercell.We define the cancer initiation time σ as the time when the first successfultype 2 cell is born. Then, σ is a random variable with complimentary cumulativedistribution function given by P ( σ > t ) = E exp (cid:18) − u s (cid:90) t | χ t | dt (cid:19) , where | χ t | is the area of type-1 cells at time t . Simpson’s Index, a traditional non-spatial measure of heterogeneity, is defined as theprobability that two individuals, sampled at random from a population, are geneti-cally identical. More precisely, if there are N types of individuals in a population,the Simpson’s Index is defined as R = N (cid:88) i =1 (cid:18) Y i Y (cid:19) , (3)where Y i is the number of individuals of the i -th type, and Y is the size of the entirepopulation. Although this measure is usually used to characterize well-mixed popu-lations, we investigate here how it evolves over time within the spatially structuredpopulation described by the mesoscopic spatial model from Section 2. In the cancersetting, one question of interest is to determine the degree of heterogeneity of the5remalignant cell population. In our mesoscopic model, suppose there are N t type-1clones present at time t >
0. We then extend definition (3) as the time-dependentquantity R ( t ) = N t (cid:88) i =1 (cid:18) Y ,i ( t ) Y ( t ) (cid:19) , (4)where N t is a Poisson random variable with parameter N u st , Y ,i denotes thevolume of the type-1 subclone originating from the i th type 1 mutation, and Y ( t ) isthe total volume of all the type-1 families present at time t , i.e. Y ( t ) = (cid:80) N t k =1 Y ,i ( t ).From conditions (A0)-(A3) and Theorem 4 of [8] we know that overlaps betweendistinct type-1 clones occur with negligible probability by time σ ; we thereforeignore this possibility in the computation of Simpson’s Index.Building on the theory of size-biased permutations, it is possible to characterizethe distribution of R ( t ) as follows. Proposition 3.1.
The conditional expectation of Simpson’s Index for the spatialmesoscopic model is E [ R ( t ) | N t = n ] = n E (cid:34)(cid:18) S S n (cid:19) (cid:35) , (5) where S n := B + . . . + B n with B i are i.i.d. Beta ( d , random variables. The proof of this result is found in Appendix A.2.
Proposition 3.2.
The conditional variance of the Simpson’s Index is bounded asfollows [ E ( R ( t ) | N t = n )] ≤ n (cid:90) ∞ (cid:90) ∞ (cid:16) rx (cid:17) ν ( r ) ν n − ( x − r ) dx dr, (6) where ν k is the probability density function of S k as defined above. The derivation of this bound is found in Appendix A.3. Finally, the followingresult establishes the behavior of Simpson’s Index for large n . Proposition 3.3.
Conditioned on N t = n , R ( t ) converges to zero in probability as n → ∞ . This result tells us that as the number of clones increases, the probability ofselecting two cells from the same clone goes to zero.Next, we use Monte Carlo simulations to evaluate (5) and study the temporalevolution of Simpson’s Index (see Appendix A.5 for details on evaluating (5)). InFigure 1, we observe that the index first increases until it reaches a maximum, andthen starts decaying in a monotone fashion. Essentially, this result stems from thefact that in the early phase of the model the first few clones are developing andexpanding, so samples are becoming more likely to be from the same clone. As6
100 200 300 400 500 600 700 800 900 100000.10.20.30.40.50.60.7
Time S i m p s on s I nde x s=0.05s=0.1s=0.15s=0.2s=0.25 Time S i m p s on s I nde x u =7.5e − =7.5e − =7.5e − =7.5e − =7.5e − A B
Figure 1:
Time-dependence of non-spatial Simpson’s Index R ( t ) . The tem-poral evolution of the expected value of the non-spatial Simpson’s Index is shown (A) for varying values of the mutation rate u , and (B) for varying values of thefitness advantage s of preneoplastic cells over normal cells. In both simulations: M = M = 500 , N = 10 , and u = 7 . × − , s = 0 . u increases, this process of establishing mutant families anddiversification occurs faster. In particular, the maximum Simpson’s Index decreaseswith increasing mutation rate because the time periods in which only a single cloneexists become shorter. Finally, Figure 1B shows that as the selective advantage s is increased, the growth of mutant families and diversification occurs faster due tothe faster spread of mutant cells. A possible implication of this result is that moreaggressive tumors will have a higher level of heterogeneity. In the clinical setting, the spatial heterogeneity of premalignant tissues poses con-siderable challenges. It is standard practice to take one biopsy sample from anarbitrary location in a suspected premalignant tissue, and to use molecular infor-mation from this sample to determine whether it is (pre)cancerous as well as itsspecific cancer sub-type, if applicable. This information is used to help guide the di-agnosis, prediction of prognosis and treatment strategies. Due to the heterogeneityof premalignant tissue, such a single-biopsy approach may lead to incorrect subtypelabeling or diagnoses, and subsequently to suboptimal therapeutic measures. Forexample, the spatial extent of this clone is unknown and thus surgical excision orprognosis prediction may be difficult. In view of these issue, the analysis of severalbiopsies across the tumor mass upon excision seems necessary. However, this raisesanother question of the length-scale of heterogeneity : how fine or coarse should thespatial sampling be, i.e. how many sections are required for a representative geneticfingerprint of the heterogeneous tissue?In order to gain insight into these issues, we focus here on two specific clinical7uestions and introduce corresponding measures of spatial heterogeneity. • Question 1: Given a region of premalignant tissue, what is the expected length-scale of heterogeneity? (i.e. how far apart should biopsy samples be taken?) • Question 2: Provided that only a single point biopsy is available, what is theexpected size of the clone present at the biopsy?Before we introduce analytical expressions for these two measures of spatial hetero-geneity I and I , we introduce notation that will be useful below. Suppose twotype-1 mutations occur at space-time points ( x , t ) and ( y , t ), respectively. Thenthe two clones will collide at time t ∗ = t + t || x − y || c d . Define the vector v = ( y − x ) / || y − x || . Then the first interaction between the twoclones occurs at location v ∗ = x + c d v ( t ∗ − t ) = y − c d v ( t ∗ − t ) . Next define the half-spaces H + = { x ∈ R d : (cid:104) x, v ∗ (cid:105) > } and H − = { x ∈ R d : (cid:104) x, v ∗ (cid:105) < } , where (cid:104)· , ·(cid:105) is the inner product in R d . If x ∈ H + and y ∈ H − then the region ofspace influenced by the mutation that occurred at ( x , t ) is B x ,c d ( t − t ) ∩ H + , and similarly the region influenced by the mutation that occurred at ( y , t ) is givenby B y ,c d ( t − t ) ∩ H − . Note that we still have χ t = B x ,c d ( t − t ) ∪ B y ,c d ( t − t ) but we have decomposed χ t into regions influenced by the two distinct mutations. I : length scale of heterogeneity A mutation at point ( t i , x i ) generates a ball B x i ,c d ( s )( t − t i ) growing linearly in t . Thusbarring interference, at time t > t i , the type 1 family is of size Y ,i ( t ) = γ d c dd ( s ) ( t − t i ) d , (7)where γ d is the volume of the d -dimensional unit sphere. To determine the length-scale of spatial heterogeneity, consider a fixed distance r > r uniformly at random. We define I ( r, t ) to be the probability thatthese two cells are genetically identical (from the same mutant clonal expansion)8t time t . The functional dependence of I ( r, t ) on r provides an estimate of thelength scale of heterogeneity and thus may provide guidance on sampling proce-dures. For example, a suggested sampling distance r ≡ { argmin r> I ( r, t ) < . } between biopsies would ensure that sampled clones would be genetically differentfrom neighboring samples 50 percent of the time. The measure I ( r, t ) is a spatialanalog of the Simpson’s Index.The actual analysis of I ( r, t ) is quite technical so we will leave the details to theAppendix for interested readers. However, here we will provide some intuition forour approach, and also provide some graphs demonstrating the dependence of I onparameters and time. We will also provide some comparisons between our analysis(based on the mesoscopic model approximation) and simulations of I in the fullcell-based stochastic evolutionary model. Idea behind calculation.
Let a, b be the positions of two cell samples taken attime t and assume that (cid:107) a − b (cid:107) = r . Define D ab as the event that the cells atpositions a and b are genetically different at time t . Calculations in sections B.1 andB.2 of the Appendix demonstrate that P ( D ab ) only depends on the distance betweenthe samples (cid:107) a − b (cid:107) , as long as c d t + r (cid:28) L . Thus, conditioning on the location of a and b we can conclude that I ( r, t ) = 1 − P ( D ab ).Next we discuss the idea behind calculating P ( D ab ). Recall that if two clonesmeet, then each continues to spread in all directions away from the interacting clone.We will denote a cell in position x at time t by the coordinate ( x, t ). V a ( t ) and V b ( t ) are the space-time regions in which a mutation can influence the genetic stateof the samples at a and b , respectively. The union of these regions ( V a ( t ) ∪ V b ( t ))represents the space-time region in which mutations can influence the genetic stateof the samples examined at locations a or b at time t . Let E ( A ) be the number ofmutations that occur in a region A and E k the event { E ( V a ( t ) ∪ V b ( t )) = k } . Then,the event D ab can be divided into sub-events according to how many mutations haveoccurred in the spacetime region V a ( t ) ∪ V b ( t ) P ( D ab ) = ∞ (cid:88) i =1 P ( D ab ∩ E i )The following simple calculation demonstrates that the probability of more thantwo type-1 mutations occurring in the region V a ( t ) ∪ V b ( t ) is small for the carcino-genesis setting. First, we note that the volume | V a ( t ) ∪ V b ( t ) | is bounded above by | V a ( t ) | + | V b ( t ) | . In 1D this sum of volumes is c ( s ) t and in 2D it is π c ( s ) t o ,where c ( s ) and c ( s ) are the spreading speeds of the single mutant clones in di-mensions 1 and 2 respectively, provided in Section 2. Since type-1 mutations arriveinto V a ( t ) ∪ V b ( t ) as a Poisson process with rate u s , the number of mutations in V a ( t ) ∪ V b ( t ) is stochastically dominated by a Poisson random variable with rate λ = u s t in 1D and λ = u st πs log(1 /s ) 2 π t is a time prior to carcinogenesis when premalignant tissueis sampled, and in our model tumor initiation occurs at time σ when the first9uccessful cell with two mutations arises. In [8] we found that the appropriate timescale of this process is 1 /N u s . Replacing t by this in the Poisson rate λ in eachdimension we obtain λ ≡ /N u in 1D and λ ≡ π N u s log(1 /s ) in 2D.We assume the point mutation rate in healthy tissue is within the range of10 − to 10 − per base pair per cell division [22, 23, 24]. Selection advantages aremore difficult to ascertain experimentally but one study has estimated the averageadvantage s to be approximately 0.001 [25]. Lastly, the cell population sizes ofinterest in tissues at risk of initiating cancer are in the range of 10 and upward.Using these estimates we can easily calculate that the probability of more than twomutations arriving within the region of interest is negligible across all reasonableparameter ranges. Thus we can approximate P ( D ab ) with the first two terms of thesum above: P ( D ab ) ≈ P ( D ab ∩ E ) + P ( D ab ∩ E ) . (8)The exact calculations for these quantities P ( D ab ∩ E ) and P ( D ab ∩ E ) are providedin the Section B.1 and B.2 of the Appendix. Agreement with microscopic model simulations.
To verify our results we sim-ulated the full cell-based stochastic evolutionary model and compared the spatialmeasure I ( r, t ) with our derivations from the previous section. Table 1 shows theresults of these comparisons. Since the cell-based model is very computationallyintensive, only 100 simulations were performed in each set of parameter values;however the Wald confidence intervals are provided for each set of simulations. Inthis table we see a close agreement between our theoretical values for I based onthe mesoscopic model and the simulations of I based on the microscopic model.For all cases, N = 10000 and r = 10.Table 1: Comparison of I ( r, t ) between theory and simulation of the cell-basedstochastic model. u s t Theory I ( r, t ) Simulation I ( r, t ) 95% CI0.001 0.01 30 0.98 0.97 (0.94, 1)0.001 0.01 40 0.96 0.92 (0.87, 0.97)0.001 0.01 50 0.94 0.9 (0.84, 0.96)0.001 0.01 60 0.90 0.9 (0.84, 0.96)0.001 0.01 70 0.88 0.83 (0.76, 0.90)0.0001 0.1 30 0.9 0.95 (0.91, 0.99)0.0001 0.1 40 0.84 0.91 (0.85, 0.97)0.0001 0.1 50 0.8 0.84 (0.77, 0.91)0.0001 0.1 60 0.79 0.8 (0.72, 0.88)In Figures 2 and 3 we demonstrate how the spatial measure of heterogeneityvaries in time for different parameters. In Figure 2 we observe that as time in-creases I ( r, t ) decreases; this reflects the clonal expansion of existing mutant fam-ilies leading to an increase in heterogeneity over time. In addition, as the distance10 increases the probability of the two samples being genetically the same decreases,as expected. Figure 3 shows this result as a function of u . As the mutation rateincreases, I ( r, t ) decreases since the heterogeneity of the tissue increases with moremutant clones emerging. In both figures we see that sensitivity of I to parameterchanges increases as r increases. This is natural since the likelihood that two cellsare identical is more likely to be altered the further apart the two cells are.
100 200 300 400 5000.80.850.90.951 Sampling Time (t ) I r=2r=5r=15r=30 Figure 2: I in 2D as a function of sampling time t . We vary the sampling radius r and set s = 0 .
01 and u = 1 e −
5, so the mutation rate is 1 e −
7. We also set themutant growth rate c d = 0 .
25. 11 −5 u I r=2r=5r=15r=30 Figure 3: I in 2D as a function of u , which contributes to the mutation rate.Mutations arise according to a Poisson process with rate u s , and we set s = 0 . r and set the sampling time t = 300 and the mutantgrowth rate c d = 0 . I : extent of a premalignant lesion Next, suppose we have obtained a premalignant (type-1) biopsy at a single point inthe tissue at time t . We would like to estimate the expected size of the correspond-ing clone. In particular we define I ( r, t ) to be the probability that an arbitrarilysampled cell at distance r is from the same clone as the original sample.In order to study I it is necessary to define the concept of ‘size-biased pick’from a sequence of random variables { X i } i ≥ . Definition 1.
A size-biased pick from the sequence ( X i ) is a random variable X [1] such that P (cid:0) X [1] = X i (cid:12)(cid:12) X , . . . X n (cid:1) = X i X + . . . + X n . Suppose that by time t there have been successful mutations at space time points { ( x i , t i ) } Ni =1 initiating populations C ,i , 1 ≤ i ≤ N . In the model description westated that the assumptions (A0)-(A3) hold throughout the paper. A consequenceof these assumptions, proved in Theorem 4 of [8], is that overlaps between distincttype-1 cell is unlikely. Thus we assume here that 1 ≤ i ≤ N , C ,i = B x i ,c d ( t − t i ) . In order to calculate I , we first choose a clone C [1] via size biased pick from thedifferent clones { C ,i } . The radius of this pick is denoted R [1] . For ease of notationwe will make the following substitution throughout the rest of the section C = C [1] and R = R [1] . We next choose a point, p , at random from the picked clone C . Wechoose a second point p , at random a distance r away from p . In other words, thepoint p is chosen at random from the circle S = { x ∈ R : | x − p | = r } .
12o calculate I we are interested in determining the probability that p is containedin C . More specifically, let us denote the center of C by x o . It is useful to define X = | p − x o | which is a random variable with state space [0 , R ]. The heterogeneitymeasure I ( r, t ) is given by: I ( r, t ) ≡ P ( p ∈ C ) = E [ P ( p ∈ C | R, X )] . (9)The following two properties are useful in determining I :(i) If X + r ≤ R then S ⊂ C . To see this consider z ∈ S then | z − x o | = | z − p + p − x o | ≤ | z − p | + | p − x o | ≤ r + X ≤ R. (ii) If R + X < r then
S ∩ C = ∅ . To see this take z ∈ C which of course implies | z − x o | ≤ R and thus | z − p | = | z − x o + x o − p | ≤ R + X < r.
We then have that P ( p ∈ C | R, X ) = , R + X < r , X + r ≤ Rφ ( X, R ) , otherwise. (10)We can use the cosine rule to see that φ ( X, R ) = 1 π cos − (cid:18) X + r − R Xr (cid:19) . (11)Substituting expressions (10) and (11) into (9) results in a formula for I that canbe easily approximated via Monte Carlo simulation. Details of this procedure areprovided in Appendix C.The heterogeneity measure I is designed to be an estimate of the extent of apremalignant lesion that has already been detected via one point biopsy. Thus it isof interest to determine the value of I at the time of detection of the premalignantcondition (which itself may be random). We hypothesize that detection of thepremalignancy may occur at a random time τ which occurs with a rate proportionalto the total man-hours of premalignant lesions. In other words, detection of thecondition is driven by the size and duration of premalignant lesion presence. Let usdefine τ with the following: P ( τ > t ) = E exp (cid:18) − µ (cid:90) t | χ t | dt (cid:19) , (12)where we recall that | χ t | is the volume of type-1 cells at time t . Display (12) tellus that detection occurs at rate µ . Note we assume that (A1-A3) hold with u s replaced by µ . 13n Section C of the Appendix we also develop a numerical approach for estimat-ing I at the random detection time τ . Interestingly enough, it is computationallyeasier to compute I ( r, τ ) than I ( r, t ). Numerical examples.
In Figure 4, I ( r, t ) is plotted as a function of r for variousvalues of u , s and t . Figure 5 shows analogous plots of I ( r, τ ) at the randomdetection time. Comparing Figures 4 and 5 we observe an interesting phenomenon.In particular when looking at I at a fixed time in Figure 4 we see that for each r and t , I is an increasing function of both u and s . This makes sense if we consider thesystem at a fixed time. Then increasing the mutation rate will increase the expectedgrowing time of any clones, i.e., they are more likely to be born earlier; therefore anyclone we select is likely to be larger, so it is more likely that the second point selecteda distance r away will be in the original clone. Similarly increasing s increases theexpected size of clones present at time t , and thus increases I ( r, t ). However,when we look at Figure 5 we see that I ( r, τ ) is decreasing in u . Interestingly byobserving the process at the random time τ we flip the dependence on the parameters u . This phenomenon results from the fact that increasing the mutation rate allowsfor detection to be caused by multiple clones, which will therefore be smaller atdetection than if the detection were driven by a single clone. A B C
Sampling Distance I (r ,t ) s=0.05s=0.1s=0.15s=0.2s=0.25 Sampling Distance I (r ,t ) u =7.5e − =7.5e − =7.5e − =7.5e − Sampling Distance I (r ,t ) t= 32.6242t= 65.2483t= 130.4967t= 195.745t= 260.9933 Figure 4: Plot of I ( r, t ) in 2D as a function of sampling radius for (A) varyingselection strength, s , (B) varying u , and (C) varying t . In all panels N = 2 e e s = 0 . u = 7 . e − t is the median of the detection time τ with µ = 2 e − B Sampling Distance I (r , (cid:111) ) s=0.05s=0.1s=0.15s=0.2s=0.25 Sampling Distance I (r , (cid:111) ) u =7.5e − =7.5e − =7.5e − =7.5e − Figure 5: Plot of I ( r, τ ) in 2D as a function of sampling radius. In panel (A) wevary the selection strength, s and in panel B we vary u . In all panels N = 2 e e τ we use µ = 2 e −
6. If not mentioned we set s = 0 . u = 7 . e − In this work we have analyzed and examined several measures of heterogeneity in aspatially structured model of carcinogenesis from healthy tissue. In particular, wefirst derived estimates of the traditionally nonspatial measure of diversity, Simpson’sIndex, in the premalignant tissue and studied how the Simpson’s Index changes intime and varies with parameters. We observed that as expected, the Simpson’sIndex decreases over time as more mutants are produced, and that this processoccurs faster in a higher mutation rate setting. The effect of selection overall is alsoto speed up this process. We also formulated and analyzed two spatially-dependentmeasures of population heterogeneity, motivated by clinical questions. In particularwe analyzed a measure ( I ) that can identify the length scale of genetic heterogeneityduring the carcinogenesis process as well as a measure ( I ) that can estimate extentof a surrounding premalignant clone given a premalignant point biopsy.We note that in this work we have confined our analysis to a two-step modelof carcinogenesis. The results can be used in a setting where more genetic hits arerequired for full malignant transformation. However our heterogeneity estimateswould apply to the population of cells with a single mutation and thus can be usedin the setting of early stages of carcinogenesis only. Work on further mutations willbe the subject of future work.These analyses facilitate a better understanding of how to interpret discrete (inboth time and space) samples from a spatially evolving population during carcino-genesis. For example, the quantity I can be calculated to help determine the ex-pected size of a premalignant lesion, given a point biopsy that is premalignant. In ad-dition the quantity I may be used to generate suggestions for optimal sample spac-ing in situations where multiple biopsies or samples are possible. Finally, we notethat although it is possible to calculate these diversity indices using computationalsimulation of similar cell-based or agent-based models, it can be extremely computa-tionally onerous to simulate such models for even small sized lattices (100x100 sites)15or a lengthy period of time such as during the process of carcinogenesis. Thereforethe heterogeneity estimates we derive based on our mesoscopic model in many casesprovide the only feasible way to estimate spatial diversity in models living on a largerlattice e.g., 1000x1000 sites or larger. Given the large number of epithelial cells in asmall area, realistic simulations to determine statistical properties of diversity mea-sures during carcinogenesis may be completely infeasible. Our results here provideanalytical or rapidly computable expressions that enable a detailed assessment ofhow these heterogeneity measures vary depending on time as well as tissue/geneticparameters such as mutation rate and selective advantage conferred by the geneticalteration. These tools can be utilized to study how tissue heterogeneity in prema-lignant conditions varies between sites and tissue types, and thus guide sampling orbiopsy procedures across various cancer types. A Details for non-spatial Simpson’s Index
A.1 Preliminary definitions and results
To characterize the distributions of the Simpson’s Index, we introduce two defini-tions. Let L , L , . . . , L n be independent, identically distributed random variableswith distribution F . Recall the definition of a sized-biased pick from Definition 1 inSection 4.2. Then we define a size-biased permutation as follows. Definition 2.
We call (cid:0) L [1] , . . . L [ n ] (cid:1) a size-biased permutation (s.b.p) of the se-quence ( L i ) if L [1] is a size-biased pick of the sequence, and for ≤ k ≤ n and P (cid:0) L [ k ] = L j (cid:12)(cid:12) L [1] , . . . , L [ k − ; L , . . . L n (cid:1) = L j ( L j (cid:54) = L [ i ] , ∀ ≤ i Proposition A.1. (Proposition 2 in [26]) For ≤ k ≤ n , let ν k be the density of S k , the sum of k i.i.d random variables with distribution F . Then P (cid:0) X [1] ∈ dx , . . . , X [ k ] ∈ dx k (cid:1) = n !( n − k )! k (cid:89) j =1 x j ν ( x j ) dx j . . .. . . (cid:90) ∞ ν n − k ( s ) k (cid:89) j =1 ( x j + . . . + x k + s ) − ds. (13) Corollary A.2. (Corollary 3 in [26]) Let T n − k = X [ k +1] + . . . + X [ n ] denote thesum of the last n − k terms in an i.i.d s.b.p of length n . Then for k = 1 , . . . , n − , P ( T n − k ∈ ds | T n − k +1 = t ) = ( n − k + 1) t − st ν ( t − s ) ν n − k ( s ) ν n − k +1 ( t ) ds. (14)16 .2 Conditional expectation Recalling that type 1 clones have a linear radial growth rate, Simpson’s index (4)can be rewritten explicitly R ( t ) = N t (cid:88) i =1 (cid:34) (1 − t i /t ) d (cid:80) N t j =1 (1 − t j /t ) d (cid:35) , (15)where { t i } N t i =1 are the points of a Poisson process with constant intensity N u st . Inparticular, we note that conditioned on N t , the t i are i.i.d and( t /t | N t ) ∼ U (0 , . We define now X i := (1 − t i /t ) d +1 and let T := N t (cid:88) i =1 X i , (16)which allows us to rewrite (15) as R ( t ) = N t (cid:88) i =1 (cid:18) X i T (cid:19) . (17)Note that conditioned on N t , the X i are i.i.d with( X | N t ) ∼ Beta (cid:18) d , (cid:19) . To see this, note first that by symmetry X i ∼ ( t i /t ) d ; using characteristic functions,it is then easy to verify that for X ∼ Beta ( α, 1) and Y ∼ U (0 , X ∼ Y n if and only if α = 1 /n . Using the above notation and recalling the notion of asize-biased pick in Definition 1, we condition (17) on N t to find( R ( t ) | N t = n ) = n (cid:88) i =1 X i T P (cid:0) X [1] = X i | X , . . . , X n (cid:1) = E (cid:18) X [1] T (cid:12)(cid:12)(cid:12)(cid:12) X , . . . , X n (cid:19) . (18)To compute the conditional expectation of Simpson’s Inded R ( t ), we take theexpectation of (18) to find E ( R ( t ) | N t = n ) = E (cid:18) X [1] T (cid:19) = (cid:90) ∞ (cid:90) ∞ (cid:16) rx (cid:17) P (cid:0) X [1] ∈ dr, T ∈ dx (cid:1) . (19)Setting k = 1, it follows now from Corollary A.2 E ( R ( t ) | N t = n ) = (cid:90) ∞ (cid:90) ∞ (cid:16) rx (cid:17) P ( T n − ∈ d ( x − r ) , T ∈ dx )= (cid:90) ∞ (cid:90) ∞ (cid:16) rx (cid:17) P ( T ∈ dx ) P ( T n − ∈ d ( x − r ) | T = x )= n (cid:90) ∞ (cid:90) ∞ (cid:16) rx (cid:17) ν ( r ) ν n − ( x − r ) dxdr. (20)17ote that the support of ν is over [0 , ν n − is over [0 , n ]. Now,by definition, ν is the pdf of Beta (cid:0) d , (cid:1) , i.e. ν ( x ) = 1 d x d − . On the other hand, ν n − is the density of the sum of n − Beta ( d , 1) randomvariables, i.e. ν n − ( x ) = (cid:16) ν ∗ ( n − (cid:17) ( x ) . For positive integer n let S n = B + . . . + B n where B i are independent Beta ( d , E [ R ( t ) | N t = n ] = n E (cid:34)(cid:18) S S n (cid:19) (cid:35) , (21) S n := B + . . . + B n where B i are independent Beta ( d , 1) random variables. A.3 Upper bound for variance We derieve an upper bound for the variance of the conditional Simpson’s Index asfollows:[ E ( R ( t ) | N t = n )] ≤ E (cid:0) R ( t ) (cid:12)(cid:12) N t = n (cid:1) = E (cid:32)(cid:20) E (cid:18) X [1] T (cid:12)(cid:12)(cid:12)(cid:12) X , . . . , X n (cid:19)(cid:21) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) N t = n (cid:33) ≤ E (cid:32) E (cid:32)(cid:20) X [1] T (cid:21) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) X , . . . , X n (cid:33) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) N t = n (cid:33) = E (cid:32)(cid:20) X [1] T (cid:21) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) N t = n (cid:33) = n (cid:90) ∞ (cid:90) ∞ (cid:16) rx (cid:17) ν ( r ) ν n − ( x − r ) dx dr, (22)where the second to last equality follows from the fact that the sub-sigma algebra σ ( N t = n ) is coarser than σ ( X , . . . , X n ). A.4 Proof of Proposition 3.3 Proof. Let Y n = ( R ( t ) | N ( t ) = n ), and note that by definition Y n ≥ 0. Thus itsuffices to show that E [ Y n ] → n → ∞ . Note that E [ Y n ] = 1 n E (cid:34)(cid:18) S S n /n (cid:19) (cid:35) , S / ( S n /n ) → S / E [ B ] as n → ∞ . Thus if weestablish that sup n< ∞ E (cid:34)(cid:18) S S n /n (cid:19) (cid:35) < ∞ , (23)then by uniform integrability we will have E [( S / ( S n /n )) ] → 1, and thus E [ Y n ] → S ,n = B + . . . + B n and for ε > A n = { S ,n > (1 − ε )( n − E [ B ] } . We then have that E (cid:34)(cid:18) S S n /n (cid:19) (cid:35) = n E (cid:34)(cid:18) S S + S ,n (cid:19) (cid:35) = n E (cid:34)(cid:18) S S + S ,n (cid:19) ; A n (cid:35) + n E (cid:34)(cid:18) S S + S ,n (cid:19) ; A cn (cid:35) ≤ O (1) + n P ( A cn ) . From Azuma-Hoeffding inequality we know that there exists a k independent of n such that P ( A cn ) ≤ e − kn , thus establishing (23). (cid:4) A.5 Monte Carlo Simulations We evaluate the conditional expectation of Simpson’s Index for fixed time t usingMonte Carlo simulations. Based on the representation in (5), we first generate M independent copies of the vector ( S , S n ) denoted by { ( S ( i )1 , S ( i ) n ) } Mi =1 , and form theestimator ˆ µ ( n, M ) = nM M (cid:88) i =1 (cid:32) S ( i )1 S ( i ) n (cid:33) , which satisfies E [ˆ µ ( n, M )] = E [ R ( t ) | N t = n ] and V ar [ˆ µ ( n, M )] = O (1 /M ). If wesimulate M copies of N t , denoted by { N ( i ) t } M i =1 and then for each realization N t = n we form the estimator ˆ µ ( n, M ) then we have an unbiased estimator of E [ R ( t )] viaˆ R ( M , M ) = 1 M M (cid:88) j =1 ∞ (cid:88) n =0 { N ( j ) t = n } ˆ µ ( n, M ) , since N ( j ) t is independent of ˆ µ ( n, M ). Note that simulating the mesoscopic model M times and averaging R ( t ) over those simulations is equivalent to using the estimatorˆ R ( M, I Calculations Recall from Section 4.1 that I ( r, t ) is approximated by (8). It is therefore necessaryto calculate P ( D ab ∩ E ) and P ( D ab ∩ E ). Also recall that V x ( t ) is the space-timecone centered at x and has radius c d t at time 0 and radius 0 at time t . For twopoints a and b in our spatial domain we will be interested in the sets D ( r, t ) = V a ( t )∆ V b ( t ) M ( r, t ) = V a ( t ) ∩ V b ( t ) , where r = (cid:107) a − b (cid:107) . We suppress the dependence on a and b in D and M to emphasizethat the volume of these sets depends only on the distance (cid:107) a − b (cid:107) . Denote theLebesgue measure of a set A ∈ R d × [0 , ∞ ) by | A | . In order to calculate I ( r, t ) itwill be necessary to compute | D ( r, t ) | and M ( r, t ) | . Note that | D ( r, t ) | = 2 ( | V a ( t ) | − | M ( r, t ) | ) . (24)In the next two subsections we compute I in one and two dimensions. For ease ofnotation we define µ = u s . For real number a define a + = max { a, } . B.1 I in 1 dimension We will first calculate the volumes | V a ( t ) | , | M ( r, t ) | , and from (24) | D ( r, t ) | . Inone dimension these calculations are simple: | V x ( t ) | = t c d , and | M ( r, t ) | = [(2 t c d − r ) + ] c d , so we have that: D ( r, t ) = 2 t c d − t c d − r ) + ] c d . Note that if t < r/ (2 c d ) then the only way sites a and b are the same at time t isif there are zero mutations in V a ( t ) ∪ V b ( t ), i.e., I ( r, t ) = exp ( − µ | V a ( t ) ∪ V b ( t ) | ) = exp (cid:0) − µt c d (cid:1) . Thus assume for the remainder of the subsection that t > r/ (2 c d ), in which case | V a ( t ) ∪ V b ( t ) | = 2 t c d − (2 t c d − r ) c d . And since the mutations arise according to a Poisson process with parameter µ , P ( E k ) = ( µ | V a ( t ) ∪ V b ( t ) | ) k e − µ | V a ( t ) ∪ V b ( t ) | k ! . From (8) it thus remains to compute P ( D ab | E ) and P ( D ab | E ). Note that ifevent E occurs then D ab can only occur if the single mutation occurs in the set D ( r, t ) and therefore P ( D ab | E ) = | D ( r, t ) || V a ( t ) ∪ V b ( t ) | = 1 − (2 t c d − r ) t c d − (2 t c d − r ) . R R R R R R R a b t Figure 6: RegionsIn order to calculate P ( D ab | E ), we must split V a ( t ) ∪ V b ( t ) into 7 differentregions because the probabilities will differ, depending on where the first mutationoccurs (as shown in Figure 6). By conditioning on E we assume that two mutationsoccur in the space-time region V a ( t ) ∪ V b ( t ). Denote the space time coordinates ofthe first mutation by ( x , t ).If ( x , t ) occurs outside of M ( r, t ) but between a and b (i.e. in regions R or R ), then the cells will definitely be different, regardless of where the secondmutation occurs. However, if the first mutation occurs in R i , 1 ≤ i ≤ 5, thenthe location of the second mutation will determine whether the sampled cells aredifferent. Thus each region R i , 1 ≤ i ≤ 5, will have an associated region Z i that willbe used to calculate P ( D ab | E ). If the first mutation occurs at the point ( x , t ) ∈ R i ,then the shape and size of Z i ( x , t ) depends on i and ( x , t ).First, we will consider the regions inside M ( r, t ), which are R , R , and R .For i = 1 , , Z i ( x , t ) represents the region in which the occurrence of a secondmutation would make the sampled cells different at time t , i.e. the two clones willmeet between a and b and then will each spread to one of the cells.If ( x , t ) ∈ R , then ( x , t ) is in M ( r, t ) and between a and b . In this case, Z ( x , t ) consists of two triangles, whose upper vertices occur at positions a and b (see Figure 7 a ). The base of the triangle on the left is 2( x − a ), and the base of thetriangle on the right is 2( b − x ), so the total area of Z ( x , t ) is c − d [( x − a ) +( b − x ) ].If ( x , t ) ∈ R , then ( x , t ) is in M ( r, t ) but to the left of a . In this case, Z ( x , t ) is a trapezoidal region. This trapezoidal region can be constructed by21 t a Z x (a) Z ( x , t ) t a b x Z (b) Z ( x , t ) Figure 7: The region in which the occurrence of a second mutation would make thecells located at a and b different. a bx t Z Figure 8: Z ( x , t ): The region in which the occurrence of a second mutation wouldmake the cells located at a and b the same.22aking the triangle whose upper vertex is at position b and subtracting the smallertriangle with upper vertex at position a (see Figure 7 b ). The base of the largertriangle is 2( b − x ), and the base of the smaller triangle is 2( a − x ). Hence, thearea of Z ( x , t ) is c − d [( b − x ) − ( a − x ) ]. Z ( x , t ) is constructed analogouslyto Z ( x , t ). Z ( x , t ) and Z ( x , t ) have a slightly different meaning. Given that the firstmutation occurs in region 4 or 5, respectively, Z ( x , t ) and Z ( x , t ) each representthe region in which the occurrence of a second mutation would make the sampledcells the genetically identical.If ( x , t ) ∈ R , then ( x , t ) is outside of M ( r, t ) and to the left of a . Inorder for a and b to be the same in this case, the second clone must meet the firstclone before it reaches a , and the second clone must spread to b before t . Hence, Z ( x , t ) is a triangle inside M ( r, t ) (see Figure 8). In the next paragraph we willexplain how the area of Z ( x , t ) is calculated.The distance between the right vertex of Z ( x , t ) and a is equal to the distancebetween a and x , so the position of that vertex is a + ( a − x ) = 2 a − x . Let V (cid:48) b be the portion of V b that falls between the t -values t and t . Then we can findthe position of the left vertex of Z ( x , t ) by considering it as the left corner of V (cid:48) b .The height of V (cid:48) b is t − t , so its base is 2 c d ( t − t ). Then the left vertex of V (cid:48) b ,and consequently the left vertex of Z ( x , t ) is b − c d ( t − t ). Hence, the base of Z ( x , t ) has length 2 a − x − b + c d ( t − t ). Therefore, the area of Z ( x , t ) is:(2 a − x − b + c d ( t − t )) c d . Analogously, the area of Z ( x , t ) is:( a − b + x + c d ( t − t )) c d . In summary, we have the following areas:If ( x , t ) ∈ R , then | Z ( x , t ) | = c − d [( x − a ) + ( b − x ) ]If ( x , t ) ∈ R , then | Z ( x , t ) | = c − d [( b − x ) − ( a − x ) ]If ( x , t ) ∈ R , then | Z ( x , t ) | = c − d [( x − a ) − ( x − b ) ]If ( x , t ) ∈ R , then | Z ( x , t ) | = (2 a − x − b + c d ( t − t )) c d If ( x , t ) ∈ R , then | Z ( x , t ) | = ( a + c d ( t − t ) − b + x ) c d X n be the position of the n th mutation. Then: P ( D ab | E ) = (cid:88) i =1 P ( X ∈ Z i | X ∈ R i ) P ( X ∈ R i )+ (cid:88) i =4 P ( X / ∈ Z i | X ∈ R i ) P ( X ∈ R i ) + (cid:88) i =6 P ( X ∈ R i ) . Thus to calculate P ( D ab ) it remains to calculate P ( X ∈ Z i | X ∈ R i ) and P ( X ∈ R i )for i ∈ { , . . . , } . t a b Ax 11 1 (a) A ( x , t ) bx t a A (b) A ( x , t ) Figure 9: The region inside V a ∪ V b that is affected by a mutation at ( x, t ) ∈ R i , andthus is not susceptible to subsequent mutation.Let A i ( x, t ) be the region inside V a ∪ V b that is affected by a mutation at ( x, t ) ∈ R i . Since type-1 mutations must occur in cells that have not yet mutated, thesecond type-1 mutation cannot occur inside A i ( x, t ).The area of A i ( x, t ) depends on whether (x,t) is in M ( r, t ) , V a \ V b , or V b \ V a . Thefollowing are the areas | A i ( x, t ) | , which will be used to calculate P ( X ∈ Z i | X ∈ R i ): | A ( x, t ) | = | A ( x, t ) | = | A ( x, t ) | = c d ( t − t + r c d ) − r + ( b − x + c d ( t − t )) + ( x − a + c d ( t − t )) c d | A ( x, t ) | = | A ( x, t ) | = c d ( t − t ) − ( x − a + c d ( t − t )) + ( a − x + c d ( t − t )) c d | A ( x, t ) | = | A ( x, t ) | = c d ( t − t ) − ( x − b + c d ( t − t )) + ( b − x + c d ( t − t )) c d . We will explain how | A ( x, t ) | is calculated and leave out the calculations for | A ( x, t ) | and | A ( x, t ) | , which can be done similarly.24 A ( x, t ) | is calculated by taking the area of the truncated triangle V (cid:48) a (the por-tion of V a that lies between times t and t ) and then subtracting the area of twosmaller triangles that are not in A (see Figure 9 b ). The bases of these triangles liealong line t , between x and the two lower vertices of V (cid:48) a . The height of V (cid:48) a is t − t , sothe base is 2 c d ( t − t ). Hence the lower left vertex of V (cid:48) a is at position a − c d ( t − t ),and the lower right vertex of V (cid:48) a is at position a + c d ( t − t ). Therefore the base of theleft small triangle is x − a + c d ( t − t ), so its area is ( x − a + c d ( t − t )) c d . The baseof the right small triangle is a + c d ( t − t ) − x , so its area is ( a − x + c d ( t − t )) c d .Since | V (cid:48) a | = c d ( t − t ) , we get the area listed above for A and A .If X = ( x, t ) ∈ R i , then P ( X ∈ Z i ) = | Z i ( x, t ) || V a ∪ V b \ A i ( x, t ) | . We can integratethis quantity over the places where the first mutation could have occurred, which isall of R i , and then divide by | R i | to get: P ( X ∈ Z i | X ∈ R i ) = 1 | R i | (cid:90) R i | Z i ( x, t ) || V a ∪ V b \ A i ( x, t ) | dxdt. Now it remains to calculate P ( X ∈ R i ). Since mutations arrive according to aPoisson process, we have P ( X ∈ R i ) = µ ( | R i | ) e − µ ( | R i | ) , and it suffices to know thefollowing areas: | R | = (2 t c d − r ) c d − ( a − b + c d t ) c d | R | = | R | = ( a − b + c d t ) c d | R | = | R | = c d t − ( a − b + c d t ) c d | R | = | R | = r c d . The expression for | R | listed above is calculated by considering R as a triangleinside V b . The height of V b is t , so the left vertex is at position b − c d t . Then thebase of R is a − b + c d t , which means its area is ( a − b + c d t ) c d .Then we can use | R | to calculate | R | and | R | : | R | = | M ( r, t ) | − | R | , and | R | = | V a | − | R | . And the height of R is t minus the height of R , so | R | simplifies to r c d .25ll of the equations above can be used to calculate P ( D ab ): P ( D ab ) ≈ P ( D ab | E ) P ( E ) + P ( E ) (cid:0) (cid:88) i =1 P ( X ∈ Z i | X ∈ R i ) P ( X ∈ R i )+ (cid:88) i =4 P ( X / ∈ Z i | X ∈ R i ) P ( X ∈ R i ) + (cid:88) i =6 P ( X ∈ R i ) (cid:1) . B.2 I in 2 dimensions Similary to the one dimensional case, we will first calculate | V a ( t ) | and | M ( r, t ) | .In the two dimensional setting this is slightly more difficult. First we know that | V a ( t ) | = πt c d / 3, so it remains to find | M ( r, t ) | . r/ θ x c d ( t − s ) Figure 10: Overlap of space time cones at time s .First observe that if r > c d t then M ( r, t ) = ∅ , so we only need to calculate | M ( r, t ) | in the case r < c d t . If we consider the overlap of space-time cones atthe fixed time s ∈ [0 , t − r/ (2 c d )] then looking at Figure 10 it can be seen thathalf the area of the overlap of their cones at this specific time is given by taking thedifference between the area of the circular section with radius c d ( t − s ) and angle θ and twice the area of the triangle with side lengths x, r/ , c d ( t − s ). The area ofthe circular section is given by c d ( t − s ) cos − (cid:18) r c d ( t − s ) (cid:19) , and twice the area of the triangle is given by r (cid:113) c d ( t − s ) − r / . Thus the area of overlap between the two cones at time s is given by a ( s ) = 2 (cid:18) c d ( t − s ) cos − (cid:18) r c d ( t − s ) (cid:19) − r (cid:113) c d ( t − s ) − r / (cid:19) . M ( r, t ) is therefore given by | M ( r, t ) | = (cid:90) t − r/ c d a ( s ) ds =2 (cid:90) t − r/ c d (cid:18) c d ( t − s ) cos − (cid:18) r c d ( t − s ) (cid:19) − r (cid:113) c d ( t − s ) − r / (cid:19) ds = 2 c d (cid:90) c d t r/ y cos − (cid:18) r y (cid:19) dy − rc d (cid:90) c d t r/ (cid:112) y − r / dy. Applying integration by parts to the first integral we see that2 c d (cid:90) c d tr/ y cos − (cid:18) r y (cid:19) dy = 2 c d t − (cid:18) r c d t (cid:19) − r c d (cid:90) c d t r/ y (cid:112) y − r / dy. We thus have that | M ( r, t ) | = 2 c d t − (cid:18) r c d t (cid:19) − r c d (cid:90) c d t r/ y − r (cid:112) y − r dy = 2 c d t − (cid:18) r c d t (cid:19) − rt (cid:113) c d t − r − r c d log r c d t + (cid:113) c d t − r , which we can combine with (24) to see that for r < c d t | D ( r, t ) | = 2 c d t (cid:18) π − − (cid:18) r c d t (cid:19)(cid:19) + 2 rt (cid:113) c d t − r + r c d log r c d t + (cid:113) c d t − r . With these calculations we see that | V a ( t ) ∪ V b ( t ) | = 2 c d t (cid:18) π − cos − (cid:18) r c d t (cid:19)(cid:19) + rt (cid:113) c d t − r + r c d log r c d t + (cid:113) c d t − r , We can now explicitly calculate P ( D ab | E ) = | D ( r,t ) || V a ( t ) ∪ V b ( t ) | . The remainder of thissection will deal with the calculation of P ( D ab | E ).The approach here will be slightly different from the one-dimensional case be-cause it is easier to look at the two-dimensional cross sections of | V a ( t ) ∪ V b ( t ) | ,rather than the entire three-dimensional space-time cones. Therefore, we will splitthe cross sections into just two regions, and then when calculating the relevant vol-umes involved in I , we will split the regions into multiple cases. In the end, the27 b m m Figure 11: When two mutation circles collide, they will continue to expand alongthe line perpendicular to the line segment joining the two mutation origins.process is similar, but the setup will be simpler, and then the volume calculationswill be more complicated in the two-dimensional setting.If two events occur in V a ( t ) ∪ V b ( t ), then the probabilities will differ, dependingon whether the first event occurs in M ( r, t ) or in D ( r, t ). We will assume that if twomutation circles collide, then we can draw a line through that point, perpendicularto the line segment connecting the two mutations (as show in Figure 11). The circleswill not extend beyond that line but will continue to expand in all other directions. r a r b a b Figure 12: The cross-section of the cones V a ( t ), V b ( t ), C a ( t ), and C b ( t ) at themoment when a mutation occurs in the intersection, M ( r, t ). If a second mutationoccurs in the shaded mutation, then the cells located at a and b will be different.If the first event occurs in M ( r, t ) at position ( x , y ) at time t , then let r a be the distance between ( x , y ) and a , and let r b be the distance between ( x , y )and b . Then let C a ( t ) be the cone centered at a that extends to the edge of theexpanding clone, so C a ( t ) will have radius r a at time t and radius 0 at time28 t + r a c d ). Similarly, C b ( t ) will be the cone centered at b with radius r b at time t and radius 0 at time ( t + r b c d ). Cross-sections of these cones are shown in Figure 12.If the second mutation occurs outside of C a ( t ) ∪ C b ( t ), then the first clone willreach both a and b before interacting with the second clone. If the second mutationoccurs in C a ( t ) \ C b ( t ), then the line dividing the two clones will separate a from b , so the second clone will affect a , and the first will affect b , making the two cellsdifferent. Similarly, if the second mutation occurs in C b ( t ) \ C a ( t ), then the firstclone will affect a , and the second will affect b .However if the second mutation occurs in C b ( t ) ∩ C a ( t ), then both a and b willbe on the same side of the line dividing the mutation circles, so the second clonewill affect both a and b . Therefore, the two cells will only be different if the secondmutation occurs in C b ( t ) (cid:52) C a ( t ). a b r a Figure 13: The cross-section of the cones V a ( t ), V b ( t ), and C a ( t ) at the momentwhen a mutation occurs in D ( r, t ). If a second mutation occurs in the shadedmutation, then the cells located at a and b will be the same.If the first mutation occurs in D ( r, t ), then its position ( x , y ) is closer to either a or b . Without loss of generality, say that ( x , y ) is closer to a . Again let r a bethe distance between ( x , y ) and a , and let C a ( t ) be the cone centered at a withradius r a at time t and radius 0 at time ( t + r a c d ). A cross-section of this cone isshown in Figure 13.If the second mutation occurs outside of C a ( t ), then the first mutation willreach a before interacting with the second mutation. Since the first mutation isoutside V b ( t ), it cannot reach b by time t , so the sampled cells will be different.If the second mutation occurs inside C a ( t ) but outside M ( r, t ), then the twomutations will interact before the first mutation reaches a , meaning that the secondmutation will affect a . However, the second mutation will not spread to b , since itdoes not start in V b ( t ). Hence, the cells located at a and b will only be the same ifthe second mutation occurs in C a ( t ) ∩ M ( r, t ).29n summary, we have: P ( D ab | E ) = P ( X ∈ C b ( t ) (cid:52) C a ( t ) | X ∈ M ( r, t )) P ( X ∈ M ( r, t )) (25)+ P ( X / ∈ C a ( t ) ∩ M ( r, t ) | X ∈ V a ( t ) \ V b ( t )) P ( X ∈ V a ( t ) \ V b ( t ))+ P ( X / ∈ C b ( t ) ∩ M ( r, t ) | X ∈ V b ( t ) \ V a ( t )) P ( X ∈ V b ( t ) \ V a ( t ))Since the mutations arise according to a homogenenous Poisson process, we canuse the volume calculations for M ( r, t ) and V a ( t ) \ V b ( t ) to calculate the followingprobabilities: P ( X ∈ M ( r, t )) = µ ( | M ( r, t ) | ) e − µ ( | M ( r,t ) | ) P ( X ∈ V a ( t ) \ V b ( t )) = P ( X ∈ V b ( t ) \ V a ( t )) = µ ( | V a ( t ) \ V b ( t ) | ) e − µ ( | V a ( t ) \ V b ( t ) | ) Similarly to the 1-D case: P ( X ∈ C b ( t ) (cid:52) C a ( t ) | X ∈ M ( r, t )) (26)= 1 | M ( r, t ) | (cid:90) M ( r,t ) | C b ( t ) (cid:52) C a ( t ) || V a ( t ) ∪ V b ( t ) \ A ( x, y, t ) | dxdydt, where A ( x, y, t ) is the cone-shaped region inside V a ( t ) ∪ V b ( t ) that is affected by amutation at ( x, y, t ).In addition: P ( X ∈ C a ( t ) ∩ M ( r, t ) | X ∈ V a ( t ) \ V b ( t )) (27)= 1 | V a ( t ) \ V b ( t ) | (cid:90) V a ( t ) \ V b ( t ) | C a ( t ) ∩ M ( r, t ) || V a ( t ) ∪ V b ( t ) \ A ( x, y, t ) | dxdydt. We next develop formulas to compute the volumes in the previous two displays.Note that | C b ( t ) (cid:52) C a ( t ) | = | C b ( t ) | + | C a ( t ) | − | C b ( t ) ∩ C a ( t ) | and that | C a ( t ) | = π r a (cid:18) t + r a c d (cid:19) , since C a ( t ) is a cone with radius r a and height t + r a c d . Similarly, | C b ( t ) | = π r b (cid:18) t + r b c d (cid:19) .We next compute | C a ( t ) ∩ C a ( t ) | . A cross-section of C a ( t ) has radius r a − ( s − t ) c d at time s , and a cross-section of C b ( t ) has radius r b − ( s − t ) c d at time s . C a ( t )and C b ( t ) will have a nonempty intersection until r a − ( s − t ) c d + r b − ( s − t ) c d = r ,i.e. when s = r a + r b − r c d + t . 30f we denote the area of intersection of the cross-sections of C a ( t ) and C b ( t )at time s by I ( s ), then | C b ( t ) ∩ C a ( t ) | = (cid:90) ra + rb − r cd + t t I ( s ) ds. (28) I ( s ) is calculated by summing the areas of the two circular segments, each ofwhich can be calculated by subtracting the area of a triangle from the area of awedge of the circle I ( s ) = R a ( s ) cos − (cid:18) d a ( s ) R a ( s ) (cid:19) − d a ( s ) (cid:112) R a ( s ) − d a ( s )+ R b ( s ) cos − (cid:18) d b ( s ) R b ( s ) (cid:19) − d b ( s ) (cid:113) R b ( s ) − d b ( s ) , where R a ( s ) = r a − ( s − t ) c d , R b ( s ) = r b − ( s − t ) c d , d a ( s ) = r − R b ( s ) + R a ( s )2 r ,and d b = r + R b − R a r . In order to compute the quantity | C a ( t ) ∩ M ( r, t ) | used in equation 27, we needto first determine when the time cross sections of C a ( t ) and V b ( t ) have nonemptyintersection. This occurs when r a − ( s − t ) c d + c d ( t − s ) > r , i.e. when s < (cid:18) r a − rc d + t + t (cid:19) . Hence: | C a ( t ) ∩ M ( r, t ) | = (cid:90) ( ra − rcd + t + t ) t ˆ I ( s ) ds, (29)where ˆ I ( s ) = R a ( s ) cos − (cid:32) ˆ d a ( s ) R a ( s ) (cid:33) − ˆ d a ( s ) (cid:113) R a ( s ) − ˆ d a ( s )+ ˆ R b ( s ) cos − (cid:32) ˆ d b ( s )ˆ R b ( s ) (cid:33) − ˆ d b ( s ) (cid:113) ˆ R b ( s ) − ˆ d b ( s ) .R a is defined above, ˆ R b ( s ) = c d ( t − s ), ˆ d a ( s ) = r − ˆ R b ( s ) + R a ( s )2 r , and ˆ d b ( s ) = r + ˆ R b ( s ) − R a ( s )2 r .We finally compute | ( V a ( t ) ∪ V b ( t )) \ A ( x, y, t ) | . In pursuit of this define U asthe region that is affected by the mutation at ( x , y , t ), i.e., U = { ( x, y, s ) : | ( x, y ) − ( x , y ) | ≤ c d ( s − t ) , t ≤ s ≤ t } . u ( s ) be the cross-section of U at time s , i.e., u ( s ) = { ( x, y ) : | ( x, y ) − ( x , y ) | ≤ c d ( s − t ) } . Observe that A ( x , y , t ) is the region inside V a ( t ) ∪ V b ( t ) that is affected by amutation at ( x , y , t ), so A ( x , y , t ) = U ∩ ( V a ( t ) ∪ V b ( t )). This of courseimplies that | ( V a ( t ) ∪ V b ( t )) \ A ( x , y , t ) | = | V a ( t ) ∪ V b ( t ) | − | A ( x , y , t ) | , It thus remains to find | A ( x , y , t ) | . This will be accomplished by looking atthe cross-sections of this set for each fixed time s . Define v a ( s ) and v b ( s ) as thecross section of V a and V b , respectively, at time s , i.e., v a ( s ) = { ( x, y ) : | ( x, y ) − a | ≤ c d ( t − s ) } v b ( s ) = { ( x, y ) : | ( x, y ) − b | ≤ c d ( t − s ) } . If ( x , y , t ) ∈ V a ( t ) \ V b ( t ), then U will not intersect V b ( t ), so in this case A ( x , y , t ) = U ∩ V a ( t ). In order to compute the volume of this set we look at thearea of the cross section for each fixed time point. There are three distinct regionsfor the behavior of the area of this cross section. In the first section of time thecross section of U is contained in the cross section of V a : u ( s ) ∩ v a ( s ) = u ( s ) ⇐⇒ u ( s ) ⊂ v a ( s ) ⇐⇒ r a + c d ( s − t ) < c d ( t − s ) ⇐⇒ s < t + t − r a c d . In the final section of time the cross section of V a is contained in the cross sectionof U : u ( s ) ∩ v a ( s ) = v a ( s ) ⇐⇒ v a ( s ) ⊂ u ( s ) ⇐⇒ r a + c d ( t − s ) < c d ( s − t ) ⇐⇒ s > t + t r a c d . When t + t − r a c d < s < t + t r a c d : | u ( s ) ∩ v a ( s ) | = R u s ) cos − (cid:18) d u ( s ) R u ( s ) (cid:19) − d u ( s ) (cid:112) R u ( s ) − d u ( s ) + ˆ R a ( s ) cos − (cid:32) ˆ d a ( s )ˆ R a ( s ) (cid:33) − ˆ d a ( s ) (cid:113) ˆ R a ( s ) − ˆ d a ( s ) , R u ( s ) = c d ( s − t ), ˆ R a ( s ) = c d ( t − s ), d u ( s ) = r a − ˆ R a ( s ) + R u ( s )2 r a , andˆ d a ( s ) = r a + ˆ R a ( s ) − R u ( s )2 r a .Thus for ( x , y , t ) ∈ V a ( t ) \ V b ( t ): | A ( x , y , t ) | (30)= (cid:90) t t − ra cd t | u ( s ) | ds + (cid:90) t t + ra cdt t − ra cd | u ( s ) ∩ v a ( s ) | ds + (cid:90) t t t + ra cd | v a ( s ) | ds | A ( x , y , t ) | is computed analogously when ( x , y , t ) ∈ V b ( t ) \ V a ( t ).It remains to compute | A ( x , y , t ) | when ( x , y , t ) ∈ V b ( t ) ∩ V a ( t ). Firstnote that if ( x , y , t ) ∈ V b ( t ) ∩ V a ( t ): u ( s ) ∩ ( v a ( s ) ∪ v b ( s )) = u ( s ) ⇐⇒ u ( s ) ⊂ ( v a ( s ) ∪ v b ( s )) ⇐⇒ s < t + t − min { r a , r b } c d . = s . Once the cross-sections v a ( s ) and v b ( s ) are no longer intersecting i.e., when s > t − r c d , then | u ( s ) ∩ ( v a ( s ) ∪ v b ( s )) | = | u ( s ) ∩ v a ( s ) | + | u ( s ) ∩ v b ( s ) | . The two quantities | u ( s ) ∩ v a ( s ) | and | u ( s ) ∩ v b ( s ) | can be calculated as shownabove.Then for t + t − min { r a , r b } c d < s < t − r c d : | u ( s ) ∩ ( v a ( s ) ∪ v b ( s )) | = | u ( s ) ∩ v a ( s ) | + | u ( s ) ∩ v b ( s ) | − | u ( s ) ∩ v a ( s ) ∩ v b ( s ) | . | u ( s ) ∩ v a ( s ) | and | u ( s ) ∩ v b ( s ) | can be calculated as shown above, and thequantity | u ( s ) ∩ v a ( s ) ∩ v b ( s ) | can be calculated as shown in [27].Thus if ( x , y , t ) ∈ V a ( t ) ∩ V b ( t ), then | A ( x , y , t ) | (31)= (cid:90) s t | u ( s ) | ds + (cid:90) t − r/ (2 c d ) s | u ( s ) ∩ ( v a ( s ) ∪ v b ( s )) | ds. With (30) and (31) we can compute | A ( x , y , t ) | for arbitrary ( x , y , t ). We canthen use | A ( x , y , t ) | with (29) and (28) to compute (26) and (27). Finally we use(26) and (27) compute P ( D ab | E ) based on (25).33 I Calculations In this section we describe how to compute I ( r, t ) and I ( r, τ ). First recall fromSection 4.2 that R is the radius of the clone, Y , chosen according to a size-biasedpick, and X is the distance of p (a point selected at random from Y ) from the centerof Y .We first describe how to estimate I ( r, t ) based on (11). In particular, we cangenerate K i.i.d copies of the vector ( X, R ), denoted by { ( X i , R i ) } Ki =1 . Our methodfor generating ( X , R ) based on the time interval [0 , t ] is as follows. First generatethe arrival times of mutations based on a Poisson process with rate N u s , denotethis set of times by t , . . . , t n .. Then for each mutation calculate the size of itsfamily at time t using the formula (7), this gives us the collection of family sizes Y , , . . . , Y ,n of clones C , , . . . , C ,n . Choose a clone C = C [1] via a size biasedpick from the collection C , , . . . , C ,n , and set R to be the radius of C . Let U be auniform random variable on [0 , 1] independent of R and set X = R √ U . With thesesamples form the estimatorˆ I ( r, t ) = 1 K K (cid:88) i =1 P ( p ∈ C | R i , X i ) . We can also derive an alternative representation for P ( p ∈ C ) that is moresuitable for mathematical analysis. Denote the conditional density of X given R = y by f X ( x | R = y ) and the density of R by f R . It’s easy to see that f X ( x | R = y ) = xy for x ∈ (0 , y ) and 0 otherwise, and therefore P ( p ∈ C ) = (cid:90) ∞ r (cid:90) y − r xy f R ( y ) dxdy + (cid:90) rr/ (cid:90) yr − y φ ( x, y ) 2 xy f R ( y ) dxdy + (cid:90) ∞ r (cid:90) yy − r φ ( x, y ) 2 xy f R ( y ) dxdy = (cid:90) ∞ r ( y − r ) y f R ( y ) dy + (cid:90) rr/ (cid:90) yr − y φ ( x, y ) 2 xy f R ( y ) dxdy + (cid:90) ∞ r (cid:90) yy − r φ ( x, y ) 2 xy f R ( y ) dxdy. (32)Define ψ r ( R ) = 2 R (cid:90) R | R − r | xφ ( x, R ) dx and Φ r ( R ) = ( R − r ) R + ψ r ( R ) , R ≥ rψ r ( R ) , R ∈ ( r/ , r )0 , R ≤ r/ . Therefore we see from (32) that we have P ( p ∈ C ) = E [Φ r ( R )].The formula P ( p ∈ C ) = E [Φ r ( R )] is difficult to work with due to the complexdistribution of R . However, an interesting observation is that the distribution of34 becomes much simpler if we assume that the sampling occurs at the randomdetection time τ . In this case define R ( τ ) to be the radius of the clone that wechoose at time τ , then we can use equation (9) in [9] to see that conditional on τ = t , R ( τ ) has density f ( x | t ) = µγ d x d c d (1 − e − θt d +1 ) exp (cid:20) − µγ d r d +1 c d ( d + 1) (cid:21) for x ≤ c d t and zero otherwise. In the conditional density above θ = µγc dd / ( d + 1).In order to describe the distribution of R ( τ ) we then need the distribution of τ ,which we can get from (4) of [9]. In particular define φ ( t ) = 1 t (cid:90) t exp (cid:16) − θr d +1 (cid:17) dr, and λ = N u s then τ has density f τ ( t ) = λe tλ ( φ ( t ) − (cid:16) − e − θt d +1 (cid:17) . Therefore we can calculate that P ( R ( τ ) > z ) = µγ d λc d (cid:90) ∞ z/c d (cid:90) c d tz r d exp (cid:20) − µγ d r d +1 c d ( d + 1) (cid:21) dre tλ ( φ ( t ) − dt = λ (cid:90) ∞ z/c d (cid:16) exp (cid:16) − θ ( z/c d ) d +1 (cid:17) − exp (cid:16) − θt d +1 (cid:17)(cid:17) e tλ ( φ ( t ) − dt = exp [ λ ( z/c d )( φ ( z/c d ) − − λ (cid:104) − exp (cid:16) − θ ( z/c d ) d +1 (cid:17)(cid:105) (cid:90) ∞ z/c d e λt ( φ ( t ) − dt. Furthermore we can take derivatives to find that R ( τ ) has density given byˆ f R ( z ) = λθ ( d + 1) c d +1 d z d exp (cid:34) − θ (cid:18) zc d (cid:19) d +1 (cid:35) (cid:90) ∞ z/c d e tλ ( φ ( t ) − dt. Note that the density ˆ f is very similar to the Weibull density, and thus we cangenerate samples from ˆ f by using the acceptance rejection algorithm with a proposaldistribution based on the Weibull distribution. With these samples from the densityˆ f then we can use the function Φ r to estimate I ( r, τ ).Note that when approximating I ( r, τ ) it is not necessary to simulate the meso-scopic model. We simply generate random variables according to the density ˆ f R and then evaluate the function Φ r ( R ). However, when approximating I ( r, t ) itis necessary to simulate the mesoscopic model and thus a greater computationalburden. References [1] N. Rahman. Realizing the promise of cancer predisposition genes. Nature ,505(7483):302–308, 2014. 352] C.P. Wild, A. Scalbert, and Z. Herceg. Measuring the exposome: a powerfulbasis for evaluating environmental exposures and cancer risk. Environmentaland molecular mutagenesis , 54(7):480–499, 2013.[3] R.J. Gillies and R.A. Gatenby. Metabolism and its sequelae in cancer evolutionand therapy. The Cancer Journal , 21(2):88–96, 2015.[4] I. Bozic, T. Antal, H. Ohtsuki, and et. al. Accumulation of driver and passengermutations during tumor progression. Proceedings of the National Academy ofSciences , 107(43):18545–18550, 2010.[5] N. McGranahan and C. Swanton. Biological and therapeutic impact of intra-tumor heterogeneity in cancer evolution. Cancer cell , 27(1):15–26, 2015.[6] C.C. Maley, P.C. Galipeau, J.C. Finley, and et. al. Genetic clonal diversitypredicts progression to esophageal adenocarcinoma. Nature genetics , 38(4):468–473, 2006.[7] R. Durrett and S. Moseley. Spatial Moran models I. stochastic tunneling in theneutral case. The Annals of Applied Probability , 25(1):104–115, 2015.[8] R. Durrett, J. Foo, and K. Leder. Spatial Moran models II. Cancer initiationin spatially structured tissue. Journal of Mathematical Biology , page in press,2015.[9] J. Foo, K. Leder, and M.D. Ryser. Multifocality and recurrence risk: a quanti-tative model of field cancerization. Journal of theoretical biology , 355:170–184,2014.[10] R.H Whittaker. Evolution and measurement of species diversity. Taxon , pages213–251, 1972.[11] Y. Iwasa and F. Michor. Evolutionary dynamics of intratumor heterogeneity. PLoS One , 6:e17866, 2011.[12] R. Durrett, J. Foo, K. Leder, J. Mayberry, and F. Michor. Intratumor hetero-geneity in evolutionary models of tumor progression. Genetics , 188:461–477,2011.[13] A. Dhawan, T.A. Graham, and A.G. Fletcher. A computational modellingapproach for deriving biomarkers to predict cancer risk in premalignant disease. bioRxiv , page 020222, 2015.[14] N. Komarova. Spatial stochastic models for cancer initiation and progression. Bull. Math. Biol. , 68:1573–1599, 2006.[15] M. Nowak, Y. Michor, and Y. Iwasa. The linear process of somatic evolution. PNAS , 100:14966–14969, 2003.[16] T. Williams and R. Bjerknes. Stochastic model for abnormal clone spreadthrough epithelial basal layer. Nature , 236:19–21, 1972.3617] C. Thalhauser, J. Lowengrub, D. Stupack, and N. Komarova. Selection inspatial stochastic models of cancer: Migration as a key modulator of fitness. Biology Direct , 5:21, 2010.[18] N. Komarova. Spatial stochastic models of cancer: Fitness, migration, invasion. Mathematical Biosciences and Engineering , 10:761–775, 2013.[19] K. Sprouffske, J.W. Pepper, and C.C. Maley. Accurate reconstruction of thetemporal order of mutations in neoplastic progression. Cancer Prevention Re-search , 4(7):1135–1144, 2011.[20] M. Bramson and D. Griffeath. On the Williams-Bjerknes tumor growth model:II. Mathematical Proceedings of the Cambridge Philosophical Society , 88:339–357, 1980.[21] M. Bramson and D. Griffeath. On the Williams-Bjerknes tumour growth model:I. Annals of Probability , 9:173–185, 1981.[22] Nachmann M and Crowell S. Estimate of the mutation rate per nucleotide inhumans. Genetics , 156(1):287–304, 2000.[23] S Kumar and S Subramanian. Mutation rates in mammalian genomes. PNAS ,99(2):803–808, 2002.[24] Miyamoto MM Baer CF and Denver DR. Mutation rate variation in multi-cellular eukaryotes: causes and consequences. Nature Rev. Genet. , 8:619?631,2007.[25] I Bozic and et al. Accumulation of driver and passenger mutations duringtumor progressio. PNAS , 107(43):18545?18550, 2010.[26] J. Pitman and N.M. Tran. Size biased permutations of a finite sequence withindependent and identically distributed terms.