The role of cell location and spatial gradients in the evolutionary dynamics of colon and intestinal crypts

Background Colon and intestinal crypts serve as an important model system for adult stem cell proliferation and differentiation. We develop a spatial stochastic model to study the rate of somatic evolution in a normal crypt, focusing on the production of two-hit mutants that inactivate a tumor suppressor gene. We investigate the effect of cell division pattern along the crypt on mutant production, assuming that the division rate of each cell depends on its location. Results We find that higher probability of division at the bottom of the crypt, where the stem cells are located, leads to a higher rate of double-hit mutant production. The optimal case for delaying mutations occurs when most of the cell divisions happen at the top of the crypt. We further consider an optimization problem where the “evolutionary” penalty for double-hit mutant generation is complemented with a “functional” penalty that assures that fully differentiated cells at the top of the crypt cannot divide. Conclusion The trade-off between the two types of objectives leads to the selection of an intermediate division pattern, where the cells in the middle of the crypt divide with the highest rate. This matches the pattern of cell divisions obtained experimentally in murine crypts. Reviewers This article was reviewed by David Axelrod (nominated by an Editorial Board member, Marek Kimmel), Yang Kuang and Anna Marciniak-Czochra. For the full reviews, please go to the Reviewers’ comments section. Electronic supplementary material The online version of this article (doi:10.1186/s13062-016-0141-6) contains supplementary material, which is available to authorized users.


Supplementary material
Considering two-hit mutants at location x = 0 In the main text, our algorithm assumed that cells at the location x = 0 die at each updating time step.
Hence, the calculated probability of two-hit mutant generation excluded the two-hit mutants that may have been generated at location x = 0. Here, we compute the amended probability of two-hit mutant generation including two-hit mutants at location x = 0. Figure S1b shows the result of this assumption comparing the probability of 2-hit mutant production with P gen 2hit (T ) that excluded two-hit mutants at x = 0 ( Figure S1a). In Figure S1b, for the top high division pattern become higher than P gen 2hit (T ) in the previous model, but remains the optimal probability distribution function.

Analytical tools
In order to calculate the probability of double-hit mutant generation at each location, we first obtain the probability of 2-hit mutant production at the bottom of the crypt, i.e. the location x = n − 1.

Probability of 2-hit mutant production at the bottom of the crypt
Let p be the probability of division at the bottom of the crypt, i.e. p = p div (n − 1). We calculate the probability of two-hit mutant production for symmetric and asymmetric divisions patterns separately.

Asymmetric division
Since in the model for the asymmetric division, there is no difference between two cells in the location x = n − 1, so we assume there is only one cell at this location. At each updating step, with probability p the cell at the location n − 1 will divide. With probability u1 2 , the offspring which is placed in the same location as its wild type mother cell is a 1-hit mutant . Moreover, the offspring replacing the divided 1-hit mutant mother cell is a two-hit mutant with probability u2 2 . Thus the probability of production of 2-hit mutant in the case of asymmetric division is obtained by the following system of ordinary differential equations (homogeneous state approximation): P 0 (n − 1, 0) = 1, P 1 (n − 1, 0) = 0, P 2 (n − 1, 0) = 0, where P 0 (n − 1, t), P 1 (n − 1, t), P 2 (n − 1, t) are respectively the probability of having wild type, one-hit mutants, and double-hit mutants at time t at the location x = n − 1. This system of equations implies the probabilities of 2-hit mutant and 1-hit mutant production at the bottom of the crypt are

Symmetric division
Here for the symmetric pattern, we obtain the transition matrix for the cells at the position x = n − 1.
At each location, there are two cells, so for the location x = n − 1, there are 9 different states; S 1 := , where 0 means no mutation, 1: 1-hit mutant, and 2: two-hit mutant. We calculate the transition matrix, the value at the row i and the column j corresponds to the probability of transforming from the state i to the state j.
This transition matrix implies that the cell at the location x = n − 1 is a two-hit mutant or 1-hit mutant at time t, with the following probabilities: v = (T s ) t < 1, 0, 0, 0, 0, 0, 0, 0, 0 > (12) where v i is the i-th element of v. We numerically calculate v after t time steps.

Probability of two-hit mutation at each location
Recall that x = 0 corresponds to top of the crypt, and x = n − 1 to the bottom. We can write a recursive equation to obtain the probability of 2-hit mutant production at each location x. Note that for both symmetric and asymmetric divisions' patterns, if a division occurs at each location i < x then position of cells at the location x does not change. Furthermore, due to the spatial architecture of this model, both asymmetric and symmetric divisions result in two cells at position x and two cells position x − 1. Also if the divisions happen at the location i > x + 1, then cells at the location x + 1 would move to the location x. The probability of having no mutations at position x at time t is: Using Taking a Taylor expansion, and assuming u i (p div (x + 1) results in the following advection-reaction equation The probability of having one mutation at time x at time t is given by This can be rewritten as Again, assuming p div (x) ≈ p div (x + 1), and keeping only the highest order terms in the Taylor expansion, result in the following advection-reaction equation Similarly, we obtain an equation for P 2 (x, t) Equations (17) (8) and initial condition is given by (5).
Note that we can numerically obtain the probability of 2-hit mutant production by solving equations (17), (18), (19) in steady state: Equation (22) implies that for large t if u 2 p div (x) 1, ∂P2 ∂x ≈ 0. Hence, when the number of cells n is large and mutation rates are small P 2 (i + 1, t) ≈ P 2 (i, t) for all 1 < i < n .

Time for mutant removal is minimized when most cell divisions occur at the bottom
We also computed the removal velocity n−1 i=x+1 p div (i) of cells at position x per cell division in Figure S3.
In accord with [20] we also find that the velocity at which cells are removed from the crypt increases for cells higher in the crypt. Note, that with exception of the delta pattern, the fastest velocity occurs for the experimentally observed division pattern, and the cells at position x are removed faster when divisions are more likely in the lower part of the crypt. From the velocity, we can calculate the estimated time τ i = i/p div (i) needed to remove a mutant at position i ( Figure S3c), which again is minimized when all divisions occur at the bottom of the crypt. Note that we can also calculate the expected time to remove a mutant that arises anywhere in the crypt i p div (i)τ i . The most optimal design for eliminating already occurring mutations is once again delta, where all divisions take place at the lowest position, with experimentally observed division pattern also performing very well. However, as we showed, delta is the worst division pattern for accumulating mutations, with a two-hit mutant guaranteed to arise when the probability of getting a two-hit mutant for the experimentally observed divisions is less than half (see Figure 9).

Model System with Two Possible Division Location
Suppose that the crypt has length n = 4, where x = 0 is the top of the crypt and x = 3 is the bottom, as before. Further assume that the division probability distribution only has two nonzero values, p(1) = p 1 and p(2) = p 2 , such that p 1 + p 2 = 1 (see Figure 10). Here p div (i) = p i . We assume strictly asymmetric divisions and therefore only considering one row of cells is sufficient. What is the optimal pair (p 1 , p 2 ) to minimize double-hit mutant production?
Since cells can only divide at locations x = 1 and x = 2, the system can be described as a Markov walk over the following 5 states: (0, 0), (1, 0), (0, 1), (1, 1), and E, where the first four states are scribed as the pairs of zeros and ones, where zero corresponds to the wild-type cells and ones to one-hit mutant cells, at locations (1, 2). The last state corresponds to the existence of a double-hit mutant. The 5 × 5 transition matrix is given by The initial state is given by a row-vector X 0 = (1, 0, 0, 0, 0), where each entry is the probability to find the system in each of the 5 states. The probabilities after t divisions are given by X 0 M t .
We can also consider a 4-state system by setting u 2 = 0 (the same states as before except E). We This matrix is the key to understanding the difference between the top-high and bottom high scenarios.
Large values of p 2 create relatively large entries above the main diagonal, creating and retaining one-hit mutants.
Starting from X 0 = (1, 0, 0, 0), the probability to produce a double-hit mutant at step t (given that it hasn't been produced so far) is given by where vector S is defined as, S = (0, p 1 , p 2 , 1) T It is possible to show that the first term in the Taylor expansion of H t in terms of small u 1 is given by where we used p 2 = 1 − p 1 . Note that functions H n satisfy H t = t/2 both for p 1 = 0 and p 1 = 0, and for intermediate values of p 1 they have a "dip"-shape with one minimum. The minimum corresponds to p 1 = 1/2 for n = 1, and it corresponds to increasing values of p 1 > 1/2 as n increases. For large values of n and 0 < p 1 < 1 we have and the minimum of this function is located at We conclude the following. In the presence of at most two dividing locations next to each other, • The maximum of the hazard function corresponds to p 1 = 0 or p 1 = 1. These are both deltadistributions of division probabilities. They are the worst design. Interestingly ,it does not matter where the delta is located. This result generalizes to any number of cells in the system (n > 3).
• If 0 < p 1 < 1, the optimal solution will have p 1 > 1/2, that is, a "top-high" distribution leads to a lower rate of mutant production. For low rates u 2 , it takes a long time to produce a mutant, so hazard functions for large t become important. In this limit, p 1 = 1/ √ t, that is, is it close to one (but is not quite one).
We can derive an analogous analytical expression for the 1st term in the Taylor expansion of H t in terms of small u 1 : where we used p 3 = 1 − p 1 − p 2 . An example of this function is presented in figure S4, where the quantity H t with t = 10 is plotted as a function of p 1 and p 2 . We can see that it has a minimum at the interior of the simplex, relatively close to the p 1 = 1 corner.
It is interesting to see how this landscape changes with t. As expected, this function is equal to t/2 for p 1 = 1, p 2 = 1, or p 3 = 1. These are points of the maximum of this function. For t = 1, the minimum corresponds to p 1 = p 3 = 1/2. Starting from t = 2, the minimum of his function occurs in the interim of the simplex p 1 + p 2 + p 3 and depends on t. It is shown on figure S5, where it traces a curve on the simplex as t increases. For large values of t, we have p 1 > p 2 > p 3 , and lim t→∞ p 2 (t) = lim t→∞ p 3 (t) = 0.
Another way to investigate the generation of double-hit mutants is to consider the transition matrix with a nonzero u 2 , This is a 9×9 matrix among the 8 states with 0,1, and 2 one-hit mutants and the state E with a double-hit mutant. We can calculate the probability to be in state E after t updates, which is given by {X 0 M t 0 } 9 , the 9th component of the vector X 0 M t 0 . As the hazard function, this quantity reaches its maximum in the corners of the simplex, that is, for p 1 = 1, or p 2 = 1, or p 3 = 1. A simple expression can be derived for this quantity's leading term in the u-Taylor expansion, in the case where p 3 = 0 (i.e. p 2 = 1 − p 1 ): For large values of t, this quantity reaches its minimum at The maximum is given by and corresponds to p 1 = 0 or p 1 = 1. This shows the limit of applicability of the Taylor approximation: we must have For large values of t (but within t his limit), we have a very simple, linear in p 1 , formula: Figure S6 plots the comparison of the exact expression {X 0 M t 0 } 9 (dots) with approximation (25) (lines). Figure S7 plots the comparison of the approximation (25), exact expression {X 0 M t 0 } 9 and numerical simulations of the spatial Moran model. One can see that the approximation becomes better for smaller values of the mutation rate.

Conclusions
Our findings can be summarized as follows.
• As before, the maximum production of double-hit mutants corresponds to p 1 = 1, or p 2 = 1, or p 3 = 1 (the corners of the simplex, the delta-distributions, regardless of their location).
• The optimal distribution must have a nonzero probability at all three locations (x = 1, 2, 3). Restricting the support of the distribution function increases the associated rate of double-hit mutant production.  Figure S2. Histogram of first one hit mutant production at each location in the crypt. The number of one-hit mutants at each location in the crypt after a fixed time is proportional to the probability division function at that location p div (x). Here the total number of runs=1000 and solid lines shown are 1000 · p div (x). Parameters used are u 1 = 0.001, u 2 = 0.005, t = 20, 000. Squares are symmetric divisions only (σ = 1), and circles are asymmetric (σ = 0). p 1 p 2 Figure S4. A heat-plot of the function H t with three possible division locations equation with t = 10, as a function of p 1 and p 2 (with p 3 = 1 − p 1 − p 2 ). Lighter colors correspond to higher values.