Combined changes in Wnt signaling response and contact inhibition induce altered proliferation in radiation-treated intestinal crypts.

Curative intervention is possible if colorectal cancer is identified early, underscoring the need to detect the earliest stages of malignant transformation. A candidate biomarker is the expanded proliferative zone observed in crypts before adenoma formation, also found in irradiated crypts. However, the underlying driving mechanism for this is not known. Wnt signaling is a key regulator of proliferation, and elevated Wnt signaling is implicated in cancer. Nonetheless, how cells differentiate Wnt signals of varying strengths is not understood. We use computational modeling to compare alternative hypotheses about how Wnt signaling and contact inhibition affect proliferation. Direct comparison of simulations with published experimental data revealed that the model that best reproduces proliferation patterns in normal crypts stipulates that proliferative fate and cell cycle duration are set by the Wnt stimulus experienced at birth. The model also showed that the broadened proliferation zone induced by tumorigenic radiation can be attributed to cells responding to lower Wnt concentrations and dividing at smaller volumes. Application of the model to data from irradiated crypts after an extended recovery period permitted deductions about the extent of the initial insult. Application of computational modeling to experimental data revealed how mechanisms that control cell dynamics are altered at the earliest stages of carcinogenesis.

INTRODUCTION

The intestinal crypts of Lieberkühn are closely packed, test tube–shaped invaginations that cover the surface of the intestine. Crypts are lined with a monolayer of epithelial cells arranged in a proliferative hierarchy (Figure 1a) and house stem cells that are responsible for the rapid and constant renewal of the intestinal surface. The stem cell compartment is positioned at the crypt base, comprising slow-cycling stem cells interspersed between Paneth cells in the small intestine, and related secretory cells in the colon (Sato ). Stem cells produce multipotent progenitor cells that migrate upward and differentiate to produce both absorptive and secretory cells. In the large intestine, once cells reach the crypt collar, they lose contact with the underlying basement membrane, extrude, and undergo apoptosis (Watson ; Eisenhoffer ); in the small intestine, they continue to migrate toward the tip of the villus, a finger-like projection that is connected to at least six surrounding crypts, before they are extruded.

FIGURE 1:

The structure of intestinal crypts. (a) Cartoon image of a single crypt, illustrating the decreasing concentration gradient of Wnt along the long crypt axis and highlighting the stem cell compartment, which consists of the stem and Paneth cells. Nuclei displaced to the apical surface represent mitotic cells. (b) 3D reconstruction of a single crypt, with the red surface corresponding to the lumen and the blue surface the basal surface of a crypt outlining its shape. This is used to define the dimensions of the computational crypt model (Materials and Methods). (c) 3D computational crypt model. Stem cells are blue, transit cells are yellow, differentiated cells are pink, and Paneth cells are black.

Several signaling pathways contribute to the normal regulation of these processes. A decreasing concentration gradient of Wnt signaling factors, produced by Paneth cells and the mesenchymal cells surrounding the stem cell compartment, regulates cell proliferation along the crypt–villus axis (Gaspar and Fodde, 2004). Stem cells are positioned closest to the source of Wnt ligands, which diffuse upward, so that Wnt exposure decreases for transit-amplifying and differentiated cells, triggering growth arrest and differentiation. Another important signaling pathway involved in the maintenance of proliferating cells is Notch (Fre ). Notch proteins and their receptors are active at the crypt base in the stem cell niche (Crosnier ). The combination of these and other signals coordinates maintenance of the stem and transit-amplifying cell populations and also directs binary cell fate decisions between secretory or absorptive lineages (Riccio ).

The coordinated program of cell division, migration, differentiation, and death/exfoliation ensures that the epithelial monolayer is completely renewed every few days. Furthermore, the regular clearance of cells from the epithelium ensures that cells carrying transforming mutations, which can occur frequently in this highly proliferative environment, do not remain sufficiently long enough to disrupt homeostasis. Nonetheless, colorectal cancer (CRC) is a common human disease. Tumors originate in crypts and are usually initiated by inactivation or mutation of the adenomatous polyposis coli (Apc) gene (Fodde and Brabletz, 2007; Humphries and Wright, 2008). Heterozygous germline mutation in Apc is responsible for the heritable condition familial adenomatous polyposis (FAP), which causes patients to develop numerous benign polyps in their gut lumen. These polyps typically progress to CRC, and FAP patients present with CRC earlier than sporadic cases (Alberts ; Boman ). The mechanisms responsible for the cancerous changes induced by Apc mutations involve its role as a scaffold protein in the β-catenin destruction complex: loss of the wild-type APC protein activates the canonical Wnt pathway by stabilizing β-catenin. In addition, loss of APC also directly causes defects in cell migration and adhesion due to the stabilizing effects of APC on cytoskeletal proteins, including F-actin and microtubules (Näthke, 2006).

Curative intervention is possible if CRC is identified early, which makes identifying biomarkers that permit the detection of early stages of tumor development important. One characteristic of adenomatous polyps is an expansion of the proliferative zone in crypts. In healthy crypts, the mitotic distribution peaks at positions 40–50% along the long crypt axis (Wright and Alison, 1984; Trani ). However, in hyperplastic and adenomatous crypts, this distribution broadens, such that mitotic cells are more evenly distributed along the entire crypt length (Figure 2, left and middle; Wiebecke ; Wong ; Fatehullah, Sharma, Newton, Lay, Nelson, McMahon, McIlvenny, Appleton, Cochran, Nathke, unpublished data). Of interest, the same effect has been observed in crypts after radiation treatment with both high-energy 56Fe ions and γ-rays, but before the formation of adenomas (Trani ). Such radiation treatment initially leads to increased proliferation in surviving crypts—a regenerative response that replaces cells that were killed (Wright and Alison, 1984). Immediately after radiation, tissue is highly disorganized (François ), making it difficult to investigate any change in cell size or dynamic behavior until tissue structure is restored, within 2 wk for nonlethal irradiation (Potten, 1990). However, as shown in Figure 2 (third column), the broadening of the mitotic zone persists and is still detectable 90 d after irradiation, suggesting that mutated cells generated as a result of radiation have persisted. This provides a useful biological model to determine how radiation causes changes in the distribution of mitotic cells along surviving crypts and how this leads to tumor formation. In our investigations, we concentrate on the period up to 3 mo after the initial regenerative response, when tissue organization is restored.

FIGURE 2:

Experimental data illustrating the change in the distribution of mitotic cells that occurs in response to tumor-inducing radiation. Mitotic distributions from Trani ) in crypts from the middle of the small intestine (jejunum) of male control mice (left), male mice irradiated with 4 Gy of γ-radiation after a 48-h recovery period (middle), and mice irradiated and allowed to recover for 3 mo (right). Top, raw data are plotted as a bar histogram together with a smoothed data distribution (blue curve), which is the data fitted to a nonparametric kernel-smoothing distribution (with normal distribution and a bandwidth of 10). Bottom, a sample from the smoothed distribution to illustrate the “smoothed data” that are subsequently used for parameter fitting.

It is well established that Wnt signaling is a key regulator of proliferation, but exactly how cells “interpret” Wnt signals of different strengths to decide different fates is not well understood. In addition, although increased levels of Wnt signaling resulting from mutations in Apc or β-catenin are clearly implicated in CRC, details about how much Wnt signal strength is required to affect specific cellular processes are not clear. Similarly, it is not known how subtle changes in Wnt signaling can contribute to early tumorigenesis. The broadening of the mitotic distribution, together with the known contribution of Wnt signaling to proliferation, is the focus of the investigations presented here. Specifically, using computational modeling, we determine the cellular behaviors that account for the measurable changes in the distribution of mitotic cells in response to γ-radiation. We explore different hypotheses for the effect of Wnt signaling and altered contact inhibition on proliferation.

Using experimental data that show the altered distribution of dividing cells in intestinal crypts after irradiation and before tumor development, we compare different hypotheses about factors that govern proliferation in this epithelial tissue. Each hypothesis is reflected in a unique model of division and is applied to a three-dimensional (3D) computational crypt model with a geometry constructed from 3D image data of small intestinal crypts (Appleton ; Figure 1, b and c, and Supplemental Figure S7). First, we use parameter-fitting techniques to identify the model that most accurately produces the pattern of division in untreated and irradiated crypts postrecovery. Subsequently, we compare these two cases to identify the cause of the shift in the distribution of mitotic cells after radiation treatment. Finally, we examine crypts after an extended recovery period following irradiation to determine how many “radiation-damaged” cells and/or crypts, defined as “mutant,” remain.

We find that in the optimal model for cell division, cell cycle duration is Wnt dependent, and proliferative status—whether a cell stays in cycle or exits—is determined by the Wnt concentration it experiences at division. Moreover, we find that the broadening of the proliferative zone along the crypt axis after irradiation is due to a lowering of the threshold of Wnt required for cells to proliferate and/or a loss of contact inhibition, such that cells can withstand greater compression and divide at smaller sizes. Combining these effects gives the optimal fit. Thus the effects of radiation on proliferation can be explained by changes in both cell cycle control and contact inhibition. After a 3-mo recovery period, our simulations predict that a heterogeneous population of crypts will exist consisting of individual monoclonal crypts that contain either healthy or mutated cells. Of note, in our models, we exclude changes in mechanical properties, suggesting that changes induced by radiation can be explained solely by differences in the response of cells to Wnt and to compression.

RESULTS

Experimental data

We consider the experimental data generated and published by Trani ), which we reproduce here for convenience. These data were generated in mice heterozygous for a low-penetrance Apc mutation (Apc). Unlike in animals with the more penetrant Apc genotype, the distribution of mitotic cells in crypts in control (untreated) Apc mice is indistinguishable from that in wild-type animals. At least four male mice were exposed at 6–8 wk of age with 4 Gy of whole-body γ-radiation. Their intestinal tissue was used to prepare 3D images that were then analyzed. Twenty jejunal crypts per mouse were selected randomly in 3D images and the position of mitotic cells measured relative to total crypt length after a 48-h and a 3-mo recovery period. The raw data are shown as a bar histogram in Figure 2.

For comparison with computational simulations, we applied a nonparametric kernel-smoothing algorithm to these data. This generates a distribution (Figure 2, blue curve) that represents the pattern of mitotic events in a crypt, smoothing out the noise present in the original data, which is due to the relatively small number of observations. We used the Matlab routine fitdist with a Gaussian kernel, using the default bandwidth (theoretically optimal for estimating densities for the normal distribution). The smoothed data are used in the following for parameter fitting.

Computational models

In silico experiments are conducted on a 3D lattice-free, agent-based crypt model (Materials and Methods). The geometry of the in silico crypt is defined by the dimensions of healthy small intestinal crypts from the jejunum of a 6-wk-old male wild-type (WT) mouse (Supplemental Table S1 and Supplemental Figure S7). Cells are represented by deformable spheres constrained to lie on this surface. They exert forces on one another, calculated using a log-exponential law. Cell death occurs at the crypt collar, above a threshold height. Model parameters are summarized in Supplemental Table S2.

Within this 3D framework, we seek to identify a model of cell proliferation that most accurately predicts the distribution of mitotic cells observed in control (unirradiated) crypts. This model should also inform us of the changes responsible for the alterations observed in irradiated tissue postrecovery, preceding tumor formation (Figure 2, left and middle). We compare six models that vary in the rules governing proliferative status and cell cycle duration. For comparison we include models previously applied to the study of crypt dynamics (Meineke ; van Leeuwen ; Osborne ; Buske ; Dunn ).

We first consider a simple pedigree model that assumes a generation-based approach to cell division: each proliferative cell undergoes a fixed integer number of symmetric divisions before terminally differentiating. In all cells, the length of the cell cycle is uniformly distributed, U[10, 14] hours. From here on, we refer to this as model 1. The second model (model 2) is an extension to the first and includes a longer cell cycle duration for stem cells, U[22, 26] hours, which are included as a distinct population from the transit-amplifying cells. This model was implemented in the first cell-center model of a two-dimensional (2D) cylindrical crypt, and we include it here to permit comparison (Meineke ).

The third model is also based on previous work. It assumes a linear, decreasing gradient of Wnt along the crypt axis that is normalized to 1 at the crypt base and 0 at the crypt collar (Figure 1a; Gregorieff and Clevers, 2005). A Wnt concentration threshold is defined, such that above this threshold, Wnt stimulates a full Wnt response. In this model, a cell remains proliferative only if it resides in the region with concentrations of Wnt above the threshold. Similar to models 1 and 2, we compare two scenarios (corresponding to models 3 and 4, respectively): one with uniform cell cycle durations (U[10, 14] hours) for all proliferative cells, and one with Wnt-dependent cell cycle duration (so that the length of the cell cycle is proportional to the Wnt stimulus a cell receives). In the latter case, cells located toward the base of the crypt where stem cells reside have a longer cell cycle than those positioned further upward. The cell cycle duration is U[22, 26] hours for cells at the crypt base and U[10, 14] hours at the position where the Wnt threshold is reached. Cell cycle length decreases linearly between these two spatial limits.

A limitation of models 3 and 4 is that loss of the Wnt stimulus (exposure to concentrations lower than the threshold) causes a cell to differentiate immediately. Such an abrupt exit from the cell cycle as a cell moves to a region of subthreshold Wnt is biologically unrealistic, as it permits a cell to abandon the cell cycle at any point. Despite this limitation, these models are evaluated here to allow a direct comparison to previously described computational crypt models that implement this assumption (e.g., van Leeuwen ; Osborne ; Buske ; Dunn ; Mirams ).

To overcome the limitations of models 3 and 4, we define a new model, which stipulates that proliferative status be assigned at birth, depending on whether a newly born cell receives a sufficiently high Wnt signal to remain in cycle. Unlike models 3 and 4, in this scenario, proliferative cells will always complete the cell cycle, even if they move into a region of Wnt that is below the threshold. By coupling either uniform cell cycle duration (regardless of spatial location) or Wnt-dependent cell cycle duration (as in models 3 and 4) to this Wnt response, we define two new models, models 5 and 6.

A summary of these six models is provided in Table 1. To evaluate and compare them, we sweep across the range of possible proliferation parameters and identify the optimal values for each model by comparing them to data acquired in tissue.

TABLE 1:

The six alternative models that we consider.

	Proliferative capacity			Cell cycle duration
Model	Pedigree	Spatially dependent on Wnt concentration	Spatially dependent on Wnt concentration at birth	Uniform	Wnt dependent
1	√	–	–	√	–
2	√	–	–	–	√
3	–	√	–	√	–
4	–	√	–	–	√
5	–	–	√	√	–
6	–	–	√	–	√

Each model is constructed from five separate components, or rules. The first three rules (columns 2–4) determine the proliferative status of each cell, and the final two rules (columns 5 and 6) determine cell cycle duration. The results shown in Figure 3 and Supplemental Figures S6 and S7 correspond to model 6.

FIGURE 3:

Example of a 2D parameter sweep for model 6. The effect of increasing the volume threshold for contact inhibition (x-axis) and decreasing the Wnt concentration threshold (y-axis) on the distribution of mitotic cells. The optimal parameter set to fit to the control data are highlighted in blue (0.9, 0.6), the 48-h γ-irradiated crypts in red (0.6, 0.5), and the 3-mo recovered crypts in green (0.8, 0.5). The shaded regions have an error that is within 25% of the optimal parameters. These results illustrate that the effect of irradiation within model 6 is both to decrease the Wnt threshold concentration and lower the volume threshold for contact inhibition: cells can divide at much lower volumes and under a lower Wnt stimulus to cause widening of the mitotic distribution.

We also examine the effect of density-dependent inhibition of proliferation: in each model, compression of a cell by its neighbors beyond a set amount will prevent it from attaining a threshold volume and halt progression through the cell cycle (Nurse, 1985; Gao, 1997; Dietrich ). The cell will undergo mitosis once it has attained the threshold volume (Küppers ; Dunn ; Leontieva ). We sweep across possible threshold volumes normalized to the maximum cell volume.

Optimal model of Wnt response

We seek to identify the model and parameters that most accurately reproduce the distribution of mitotic cells observed in crypts from control and irradiated mice (Figure 2, left and middle). We conduct a 2D parameter sweep for each model, varying both the Wnt concentration threshold and the volume threshold for contact inhibition. For models 1 and 2, we represent the response to a varying Wnt concentration threshold by sweeping over the average number of generations that mitotic cells are allocated before differentiating (Materials and Methods). For each parameter set, we grow an in silico crypt to a homeostatic steady state and then simulate for a further 1000 h, recording the position of each mitotic event in this 1000-h window.

We compare the histograms that describe the mitotic distribution for each of the models with the smoothed mitotic distribution derived from experimental measurements (Figure 2, bottom). This permits comparison of the error in each of the “bins” (using equal width for the simulated and experimental data), which corresponds to the difference between simulated and experimental results. Here we use the sum of squares of the differences in each bin to define the objective function for comparison.

An example of such a sweep for model 6 is shown in Figure 3. The threshold volume for contact inhibition varies along the x-axis, and the Wnt concentration threshold varies along the y-axis. Note that we omit the case for a volume threshold or Wnt threshold of 1, as in either of these scenarios, proliferation will not occur. The set of parameters that minimizes the objective function when compared with the experimental data for the control mice (Figure 2, first column) is shown in blue, corresponding to a volume threshold of 0.9 and a Wnt concentration threshold of 0.6. Thus, in this model, cells with <90% of the maximum volume will not undergo division, and only cells born in the lower 40% of the crypt will be proliferative.

Similarly, we identify the parameters for each division model that produce the best fit for the smoothed mitotic distribution observed in irradiated mice (Figure 2, middle), using the assumption that all proliferating cells will have been affected by irradiation. For model 6, a volume threshold of 0.6 and a Wnt concentration threshold of 0.5 produce the closest match to the experimental data (Figure 3, red). Compared to the parameters identified for the control case, this means that a cell can divide at a smaller volume and when experiencing a lower Wnt concentration. This produces a broadened proliferative zone, with cells dividing higher up the crypt.

Parameter sweeps were conducted for each of the six models (Figure 3 and Supplemental Figures S1–S5). For both scenarios— control and irradiated tissue—model 6 gave the smallest error between experimental and simulated data using the optimal parameters described earlier. The variation in error between the simulated and experimental data across the parameter domain is shown in Figure 4 for all six models evaluated (see also Table 2). To facilitate direct comparison of the best fitting parameters in each case, in each subplot, the blue circle identifies the parameter set that most closely matches the experimental data for untreated crypts, and the red circle is for the irradiated crypts.

FIGURE 4:

The log of the error between the simulated and experimental data varies with the parameters implemented in each model. Model numbers are indicated on each plot. The red and blue circles mark the minimum error for the control and irradiated cases, respectively. The first contour is within 25% of the minimum, the second contour within 200%, the third within 400%, and so on. In all models except model 2, the parameter sets that produce the minimum error in each case reveal that the proliferation (Wnt concentration) threshold and the volume threshold decrease from the control to the γ-irradiated case.

TABLE 2:

The parameter sets and minimum error identified for each model when compared with experimental data.

	Control			Irradiated (48 h)
Model	Volume threshold	Proliferation threshold	Error	Volume threshold	Proliferation threshold	Error
1	0.5	3	2.146 × 10−4	0.5	3.5	4 × 10−5
2	0.5	4	2.264 × 10−4	0.4	3	6.251 × 10−5
3	0.7	0.2	1.786 × 10−4	0.7	0.1	4.775 × 10−5
4	0.5	0.2	1.582 × 10−4	0	0.1	5.526 × 10−5
5	0.9	0.6	1.12 × 10−4	0.9	0.5	2.639 × 10−5
6	0.9	0.6	6.968 × 10−5	0.6	0.5	1.793 × 10−5

On the basis of these results, we can deduce the changes in contact inhibition and proliferation that are required to produce the altered mitotic distribution after irradiation. We observe that in all models except model 2, the threshold volume for contact inhibition decreases, the Wnt concentration threshold decreases, or both. These results are also reflected by the increase in the mean and SD of the simulated positions of mitotic cells in irradiated compared with control crypts (Supplemental Table S3), consistent with the changes observed experimentally. The simulated distribution of mitotic cells in model 2 does not follow this trend (i.e., the mean height decreases), suggesting that it is not an appropriate model to describe changes in proliferation in crypts.

Our result that cells in irradiated crypts have a lower threshold for contact inhibition predicts that cells will divide at smaller volumes. Consistently, for each model, the average cell volume in simulated crypts is smaller for irradiated crypts than for controls by 1–10% (Supplemental Table S4).

Overall, results from our simulations show that broadening of the distribution of mitotic cells in precancerous conditions—in this case, modeled by observations from crypts exposed to γ-irradiation—can be caused by a lowered threshold for proliferating cells to respond to Wnt signals, the ability of such cells to divide at smaller volumes, or a combination of both. We refer to these cells collectively as “mutant” cells. Thus far, we have assumed a scenario in which all proliferating cells are altered by irradiation. In reality, this is unlikely, and therefore to reproduce the shift in mitotic distribution, the characteristics we have identified for mutant cells would be more dominant—that is, more exaggerated—in the more realistically occurring, smaller number of mutant cells. Nonetheless, we can extrapolate from these results to predict that cells colonizing the epithelium after irradiation carry mutations in genes contributing to Wnt responses and cell size.

Mutant colonization of the crypt

Crypts usually are clonal (Ponder ; Bjerknes and Cheng, 1999). This means that we can ask whether the altered cell properties of mutant cells we have identified are sufficient to ensure dominance of a mutant population over healthy epithelial cells. Using model 6, which most accurately reproduces experimental data, we determined the probability that a crypt with an initially heterogeneous population of control and mutant epithelial cells becomes colonized solely by mutant cells. We compare the scenario in which mutant cells adopt the parameters identified for the 48-h γ-irradiated data (Figure 3, red) with the control case in which mutant cells are identical to untreated, healthy cells (Figure 3, blue). For each initial proportion of mutant cells, we performed 500 simulations and tracked the probability that mutant cells colonize the entire crypt in each case (Figure 5a). Note that we assume that the mechanical properties of mutant and healthy cells remain identical.

FIGURE 5:

Colonization of a crypt by mutant cells. (a) The probability that a population of mutant cells will colonize an entire crypt for different starting sizes of mutant populations. Results for mutants with parameters that mimic control cells (blue) and for mutant cells that adopt the parameters identified for cells in the 48-h γ-irradiated case (red). The shaded region represents 1 SD. (b) Simulation snapshots of a crypt with an initial heterogeneous population of 10% mutant and 90% healthy epithelial cells (blue and red, respectively; black shows Paneth cells). After 400 h, the mutant cells have colonized the crypt. (c) The average time taken (hours) for either mutant or control cells to colonize the crypt with increasing initial proportion of mutants. Shaded red region represents 1 SD.

The blue curve in Figure 5a shows that, if assigned the same proliferative properties as control cells, the probability that mutant cells colonize a crypt is equal to the initial proportion of mutant cells, as expected. Comparing the blue and red curves shows that when assigned the properties identified from the 48-h γ-irradiated data (Figure 3, red), mutant cells are more effective at colonizing the crypt than control cells, despite identical mechanical properties. Moreover, if the radiation insult affected ≥40% of proliferating cells, this is sufficient to guarantee complete conversion of a crypt to mutant cells. Under these conditions, 499 of the 500 simulations were colonized by mutant cells, and therefore the 95% confidence interval for the probability of a crypt being taken over by nonmutant cells (when 40% of the initial cells were mutant) is (−0.0019, 0.0059). It follows that there is a <1% chance (at the 95% significance level) of a crypt initially comprised of 40% mutant cells to become colonized by nonmutant cells. If a lower percentage of cells were affected by radiation, either the mutant or healthy cells could take over. Figure 6b shows an example of a simulated crypt that initially contained 10% mutant cells, which became colonized entirely by mutants after 400 h.

FIGURE 6:

Identifying the initial proportion of mutant cells that explains the mitotic distribution in recovered crypts. (a) Experimental data for the distribution of mitotic cells in irradiated crypts after a recovery period of 3 mo (smoothed data have been fit to a nonparametric kernel-smoothing distribution). (b) The log error between simulated and experimental data for the indicated percentage of mutant crypts (x-axis). The lowest error between experimental and simulated data occurs when 64% of crypts are homogeneously mutated (black circle). (c) Mitotic distribution derived for a heterogeneous population of crypts, with 64% mutant crypts and 36% control crypts (error between simulated and experimental data indicated). (d) Best-fit histogram of mitotic events in the simulated crypt to the recovered data, assuming a homogeneous population of cells (error between simulated and experimental data indicated).

Figure 5c shows the average time taken for either mutant or healthy cells to colonize a crypt, which decreases as the proportion of the corresponding cell type increases. The longest time required to reach clonality was ∼800 h (33 d), and the mean time was 586 h (24 d). This means that in the tissue used by Trani ), to investigate the mitotic distribution in crypts 3 mo postirradiation, all crypts in the samples consisted of either healthy or mutant cells. On the basis of our simulation results, we predict that all crypts in the 3-mo “recovered” tissue were dominated solely by mutant cells if at least 40% of the proliferating cells were affected by irradiation. In contrast, we expect a mixed population of mutant and healthy crypts if the radiation affected no more than 40% of proliferating cells, as either cell type could win out. That the tumor burden in the irradiated animals is 10 times higher than in untreated controls (Trani ) suggests that the number of mutant crypts was indeed high. However, since only macroscopically detectable tumors were scored, it is impossible to know how many crypts were transformed.

Recovery postirradiation

Immediately after irradiation, many cells in crypts die, and there is a rapid proliferative response to repair this injury that restores normal architecture by 2 wk for the most severe nonlethal dose (Wright and Alison, 1984; Maj ). The data showing an altered mitotic distribution even 3 mo later are consistent with the idea that mutations persisted in some of the cells during the recovery period. Our models show that by this time, individual crypts are populated entirely by either such mutant or healthy cells (Figures 5, a and c). Next we used model 6 to predict the composition of recovered tissue and sought to determine the proportion of crypts that are dominated entirely by mutant cells that most closely fits the experimental data.

We calculated the mitotic distribution for a mixed population of crypts using the individual distributions identified for single control and irradiated crypts (Figure 3, blue and red histograms). To calculate the distribution for a heterogeneous crypt population (e.g., 10% mutant crypts and 90% healthy crypts), we combined the distributions in the correct ratio (e.g., 0.1 times the mutant distribution with 0.9 times the healthy distribution). In this way, we constructed a mixed population of crypts for all ratios of healthy to mutant and compared the resulting mitotic distributions with the recovered data 3 mo postirradiation. Figure 6b shows that the error between the simulated data and the recovered experimental data is minimum for a population consisting of 64% mutant crypts and 36% control crypts. The histogram corresponding to this ratio is shown in Figure 6c.

Crucially, the error is one order of magnitude smaller than if we assumed a homogeneous population of crypts containing cells with unique proliferation parameters (i.e., if we derived new parameters following the approach used in Figure 3). This suggests that a mixed population of homogeneous crypts, individually containing either all control or all mutant cells, explains the experimental data much better than a homogeneous population of crypts. Further, we can infer that the initial burst of radiation produced mutations in ≤40% of proliferative cells in each crypt.

DISCUSSION

Cells within tissues receive instructive cues to regulate proliferation and differentiation. Changes in how this information is processed lie at the core of many diseases—particularly cancer—and can act as biomarkers for detection and therapy response. Concurrently, preceding full transformation, tissue aberrations are already detectable, such as the altered distribution of proliferating cells in preadenomatous intestinal crypts exposed to γ-radiation (Trani ). In addition to providing potential biomarkers, such changes can reveal how mechanisms that control cell behavior in healthy crypts are altered to initiate tumors.

Computational modeling is a powerful tool for testing hypotheses derived from experimental data. We describe a 3D computational model geometrically constrained according to the size, shape, and composition of small intestinal crypts. Within this framework, we compare six alternative hypotheses about the mechanisms that control cell division, which differ in how cells interpret Wnt signals to set their proliferative status and the duration of the cell cycle. We directly compare model simulations with experimental measurements to identify parameters that most accurately reproduce the situation in tissue. We focus on Wnt signaling as the major signaling pathway that regulates proliferation and differentiation in many tissues, particularly in intestinal crypts, where it is absolutely required. Changes in key proteins that regulate Wnt signaling are known to be key drivers of cancer in this tissue (Anastas and Moon, 2012; Polakis, 2012). Concurrently, we examine the role of density-dependent inhibition of mitosis to account for cell size. By sweeping over two parameters—a Wnt, or pedigree, threshold that defines when a cell is no longer proliferative, and a volume threshold that determines when the cell cycle pauses—we identify both the model and parameter set that most closely match experimental data.

Our simulations reveal that the concentration of Wnt a cell experiences when it is generated by a division (i.e., when it is born) dictates its proliferative status. Furthermore, simulations suggest that cell cycle duration is proportional to Wnt stimulus, with cell cycle times decreasing linearly along the crypt axis. We predict the Wnt concentration threshold required to maintain cells in cycle, such that cells residing in the lower 40% of the crypt receive sufficient Wnt to proliferate. We also find that, normally, the cell cycle will pause due to contact inhibition if cell volumes are <90% of equilibrium. Crucially, the same optimal model most closely recapitulates the mitotic distribution of γ-irradiated crypts 48 h postirradiation. In this precancerous situation, both the Wnt concentration threshold and the threshold volume for division are lower than in the control. This means that cells proliferate at lower Wnt concentrations and can divide despite not having reached normal size, which could happen under increased compression or because a checkpoint that links cell growth to mitotic entry is defective.

Mutations that stimulate Wnt signaling are common to almost all human tumors in intestinal tissue (Schneikert and Behrens, 2007; Polakis, 2012). Therefore the finding that precancerous cells are more sensitive to Wnt may appear inconsistent with elevated Wnt signaling in tumors (Anastas and Moon, 2012). However, an alternative way to interpret our results is that the cells act as if they perceive higher Wnt concentrations than they actually receive. Our models assume that Wnt ligands in the environment are unchanged, and parameter sweeps compare cellular response to external Wnt concentrations that vary spatially. Therefore in our simulations, a mutation that produces increased Wnt signaling in a cell is equivalent to a lowered Wnt stimulus threshold, making it appear as if a cell maintains its proliferative state at a lower Wnt concentration. Our modeling work supports the idea that mutations that cause or mimic increased Wnt signaling in addition to decreasing sensitivity to compression are sufficient to produce observable changes in proliferation patterns in precancerous crypts, suggesting that even before overt tissue changes are in place, Wnt signaling is up-regulated so that lower external Wnt concentrations can stimulate proliferation.

Similar considerations apply to the data and models for irradiated tissue. Radiation damage requires tissue repair, which involves up-regulation of Wnt signaling. Intestinal tissue damage causes local up-regulation of Wnt5a to support tissue repair (Miyoshi ). Further support for the idea that Wnt activation is involved in recovering from radiation-induced damage in crypts relates to the finding that lack of Mtg16 causes improved recovery of intestinal crypts and organoids from radiation damage (Poindexter ). Mtg16 competes with beta-catenin for binding to (and thus activating) Tcf4 (Moore ). In the absence of Mtg16 β-catenin can activate Tcf4 more effectively (akin to increased Wnt signaling) and recovery from radiation is improved. In addition, radiation damage produces free oxygen radicals, which causes stabilization of hypoxia-inducible factor α (HIF1a). Hif1a represses APC (Newton ), which in turn activates Wnt target genes. Our finding that mutant cells respond to lower Wnt concentrations creates the same situation, and cells proliferate when they normally would not. However, in our model, it is the concentration of Wnt required to stimulate proliferation that is lowered; in response to injury, an increase in the locally available Wnt could have the same effect.

A novel distinction that our models can make is whether the proliferative state of a cell is decided when it is first produced by cell division or depends critically on (and thus varies with) its spatial position. Cells in the crypt move rapidly and can cover up to 50 μm in 12 h (Nelson ). This means that the Wnt concentration that cells experience could vary significantly between divisions. We find that the Wnt concentration a cell experiences is set when it emerges from a division, which means that the position where a cell is born is crucial for its fate. This is consistent with the emerging idea that differentiation signals are most effective in the G1 stage of the cell cycle, that is, after mitosis (Dalton, 2015).

This behavior is particularly relevant for stem cells, which reside at the crypt base. In the optimal model we identify, the high Wnt concentration at the crypt base induces a long cell cycle time (22–24 h), such that stem cells can theoretically move a significant distance between birth and committing to the next cell cycle. However, stem cells may not move significantly unless they are positioned close to the stem cell niche boundary. This was recently suggested by lineage-tracing experiments, which showed that the probability of a stem cell, or its progeny, of exiting the stem cell zone is highest when it is positioned near the boundary between the stem cell niche and the transit-amplifying compartment (Ritsma ). Consistent with these data, our model shows that where a cell is born predicts whether it or its daughters will reenter the cell cycle or differentiate. Thus our models provide key insights into how signaling events in one compartment affect cellular behavior in another. This in turn explains the hierarchical organization of the intestinal crypt and also the relationship between cell position and fate.

Our simulation results show that if radiation induces a conversion to a mutant state in at least 40% of proliferating cells, these mutants will have a sufficient growth advantage and can colonize an entire crypt, even without a change in mechanical properties. For lower proportions of mutant cells, nonmutant cells can also win out. Previous modeling work found that mutant cells needed to exhibit stronger adhesion to the substrate than healthy cells to colonize the crypt; however, these investigations were constructed in a simpler, 2D geometry (Mirams ). We predict the time for conversion to be short in either case, on average requiring <800 h, consistent with idea that crypts are clonal (Ritsma ). This leads to a situation with some crypts fully mutant and some fully nonmutant after a recovery period, with the ratio between these two types of crypts dependent on how many mutant cells were created by the initial insult.

We examined this situation in the context of such “recovered” tissue that had been exposed to radiation and was then allowed to recover for 3 mo before mitotic events were recorded. In this scenario, the mitotic distribution in the crypts remained broadened, suggesting that the tissue had permanently changed. We predict that the mitotic distribution in recovered tissue reflects division events in monoclonal crypts. Our analysis confirmed the existence of two types of crypts: one comprising solely mutant, the other solely healthy cells, with a relative abundance of 64% and 36%, respectively. On the basis of these numbers, we predict that the initial insult affected <40% of proliferating cells in each crypt.

Our result is consistent with the long time it takes for the development of fully transformed tissue. The fact that the altered mitotic distribution preceded tumor development is also consistent with the idea that mutant crypts can expand and produce adenomas and with the fact that additional mutations are required to develop the necessary growth advantage to fully transform and generate tumors (Fearon and Vogelstein, 1990). Overall, the ability for radiation to induce initial tumorigenic changes is consistent with the delayed onset of colon cancer in patients receiving radiation therapy and also professionals exposed to elevated radiation, such as astronauts (Chancellor ).

Another important finding is that the response to contact inhibition is reduced after tumorigenic insult and that cells are smaller when they divide (Table 2). The complex relationship between cell cycle duration and cell size is affected by many different signaling pathways (Ginzberg ), including Wnt, and what governs the mechanical properties of cells is just beginning to be understood. The limited data available suggest that cancer cells in situ are indeed softer and more readily compressed than healthy cells, consistent with our finding that mutant cells are smaller. However, high-resolution direct measurements of mechanical properties of cancer and healthy cells in situ are currently available only in the context of breast tissue (Plodinec ). Further elucidating the relationships between different signaling pathways, mechanical properties, and cell size requires the ability to measure cell size accurately, which is nontrivial in whole tissue. This is further complicated by the fact that the maximum volume change we predict is only 10% (Supplemental Table S4). Using tissue sections to measure cell size only allows 2D measurement. If we assume that size changes result from a decrease in diameter rather than height, the maximum measurable change would be a 5% reduction in cell diameter. The curvature of crypts together with the natural variability of cell packing makes it unlikely that such a small difference can be measured reliably. An alternative, more sensitive approach is flow-activated cell sorting. However, we predict size changes in situ where compression by neighboring cells affects cell size. Thus interpreting results generated from isolated cells requires first testing the assumption that mechanical properties of mutant and nonmutant cells are identical. Tissue organoids may be a useful experimental system in which to explore this, as they are amenable to experimental manipulation and contain only epithelial cells. Indeed, when epithelial cells with homozygous mutations in Apc form organoids, they appear smaller when examined in cross section (Fatehullah ). However, how these organoids relate to the situation in irradiated crypts needs to be established before valid conclusions can be drawn.

Our results help to identify mechanisms that are disrupted at the earliest stages of tumor development, thus providing potential biomarkers for CRC. Here we consider only one signaling pathway and one of its outputs: proliferation. However, the close fit between experimental and in silico results suggests that how cells interpret Wnt is sufficient to explain their proliferative behavior, confirming that it is at the core of the regulation of intestinal tissue. On the other hand, it is well established that other signaling pathways, including Notch/Delta, bone morphogenic protein (BMP), and epidermal growth factor (EGF), also affect proliferation in the intestinal crypt, and cross-talk between them is likely to fine-tune behavior of cells. Our results suggest that it may be the ability of these pathways to modulate Wnt signaling that is key. Including additional parameters in our model could help to identify their contribution. For example, it is possible that when the Wnt concentration threshold is low, Notch signaling, which normally inhibits differentiation into the secretory lineage, becomes irrelevant, but that if the Wnt threshold is high and cells require a significant Wnt stimulus to remain proliferative, Notch signals may have a more potent effect to maintain cells dedifferentiated and closer to cycling.

MATERIALS AND METHODS

The 3D crypt model

We construct a 3D computational crypt model using a discrete agent–based approach (Osborne ; Dunn , 2012b, 2013). Cells are modeled as deformable spheres, free to move on a 3D surface defined according to experimentally measured crypt dimensions described later (Supplemental Table S1). Cell–cell connectivity is governed by an interaction radius (here taken to be 1.5 cell diameters), such that all cells whose centers are closer than this radius are considered to be connected. This connectivity is used to implement attractive and repulsive forces between neighboring cells.

We define four distinct cell types/proliferative states as shown in Figure 1c: stem (light blue), Paneth (black), transit-amplifying (yellow), and differentiated (pink). Simulations are initiated with 15 Paneth cells, which are not proliferative and are restricted to sit positioned between stem cells at the base of the crypt by applying a retainer force, and 60 proliferating cells (these are stem cells in models 1 and 2 and transit-amplifying cells in models 3–6). Paneth cells are included to ensure that we do not overestimate the number of proliferative cells, and therefore mitotic events, at the crypt base. A linearly decreasing gradient of Wnt stimulus is defined along the long crypt axis (Figure 1), from 1 at the crypt base to 0 at the crypt collar. This is used to identify the Wnt stimulus received by any given cell according to its location along the z-axis and thus to determine its proliferative status.

Stem cells are explicitly considered only in models 1 and 2, which are based on the pedigree assumption. However, in the remaining models, cells located toward the base of the crypt may have more “stem-like” properties and a longer cell cycle duration that is dependent on the strength of the Wnt stimulus. Sloughing occurs at the crypt collar, where cells are immediately removed from the simulation.

Crypt geometry

The measurements that define the surface of the crypt model are based on the dimensions of a healthy small intestinal crypt taken from the jejunum in a 6-wk-old male WT mouse (Supplemental Table S1 and Supplemental Figure S7a). The tissue was prepared for imaging as previously described (Appleton ). The sample was imaged on a Zeiss 710 multiphoton microscope and morphological measurements of the 3D data set made using Volocity (PerkinElmer) software. The 3D models of single crypts were made from 3D data sets using Imaris (Bitplane) software, and the dimensions presented here are based on averages of >500 crypts.

We use the cross-sectional area of experimentally measured crypts at 10-μm intervals up the long axis to infer the crypt radius at these points, starting 10 μm from the crypt base (Supplemental Table S1). To specify the dimensions of the crypt for all heights, we fit a fifth-degree polynomial function to these radii (another fit could be used; however, the chosen fit approximates the data points well and is not overspecified): where

This function is used to calculate the radius of the crypt, r, at a given distance along the long axis, z (Supplemental Figure S7b). We define the radius of the crypt base according to the surface of an ellipse, with a minor axis of 15.3973 and a major axis of 16.6968 μm. These axes are determined according to the radius of the crypt 10 μm from the base.

Typically, the crypts of the large intestine are longer. Whereas the dimensions used in our 3D computational model allow for accurate comparison with the radiation data, the overall behavior we observe in simulations does not change if we instead study larger intestinal crypts.

Cell–cell mechanics

Interactive cell forces, which mimic cell–cell adhesion and limited compressibility between neighboring cells, are modeled using a force law (Supplemental Figure S8) that acts along the lines between centers of neighboring cells (Pathmanathan ; Dunn ). Mechanically, all cells behave the same due to adhesion and compressibility. Let r be the position of the center of cell i. The force acting on this cell due to a neighboring cell j is where is the vector from cell center i to cell center j at time t, is the corresponding unit vector, s(t) is the natural length of the spring connecting cell centers i and j, which increases from 0.5 to 1 over the first hour of the cell cycle (van Leeuwen ), μ is the spring constant, and α is a constant that represents the level of volume exclusion, here set to 4. The constant k is a multiplication factor set to 0 if i and j are both Paneth cells and 1 otherwise. This implements differential adhesion between Paneth cells and other cells, similar to Sulsky . In addition to nearest-neighbor forces, a retainer force that acts only on Paneth and stem cells is included: where = (0, 0, 1) and γ = 75 if the cell is a Paneth or stem cell or 0 otherwise. These parameters are chosen to keep stem and Paneth cells restricted to the crypt base without inducing excessive cell compression.

Neglecting inertial terms relative to dissipative terms, the velocity of cell center i is given by where η is the drag coefficient for the motion of cell center i. As that we consider cells to be uniform, η = η for all i. By iterating in small time intervals, cell positions are updated at each time step using the forward Euler method, with the timestep chosen so that reducing it further results in indistinguishable simulations. Supplemental Table S2 gives the parameters used.

Parameter sweeps

For models 1 and 2, the equivalent of sweeping over the Wnt threshold height is to sweep over the number of possible generations each dividing cell can have (Supplemental Figures S1 and S2). The number of generations that can be assigned to dividing cells is bounded: any more than six generations will effectively prevent differentiated cells from being present in the crypt. We sweep over the average number of generations assigned to the dividing population so that when a noninteger number of generations is assigned, the fractional part of the number is used to assign the percentage of the dividing population to have the pedigree found by rounding up, and the remaining percentage will have the pedigree found by rounding down. For example, if the average number of generations is 2.7, then we assign 70% of the population to have pedigree 3 and 30% to have pedigree 2.

Density-dependent inhibition of mitosis

Density-dependent inhibition of mitosis prevents cells from progressing through the cell cycle if they do not have a sufficiently large volume. Instead, such cells “pause” the cell cycle and remain in G1 phase until their volume increases beyond the threshold. The volume of a cell is determined by the balance of all cellular forces and thus serves as a proxy for how crowded the cells are in a specific region, which is believed to affect cell division (Alberts ; Dietrich ; Küppers ).

As implemented in our earlier work (Dunn ), cell volumes are approximated by averaging the interaction distance between each spherical epithelial cell and its overlapping neighbors, with adjustments for when the cell has fewer neighbors than expected, and using this distance as the effective radius of the epithelial cell of interest, r. The volume is then calculated as that of a sphere of radius r. The equilibrium cell volume is calculated as the volume of an undeformed sphere (one cell diameter) at ∼524 μm3.

Simulations

Model simulations are conducted within the Chaste framework (Mirams ; www.cs.ox.ac.uk/chaste). Chaste is an open-source software library written in object-oriented C++ and constructed using agile programming techniques. The code used to run the simulations presented here is available from the Chaste website, where we provide a tutorial specific to this article: https://chaste.cs.ox.ac.uk/trac/wiki/PaperTutorials/CryptProliferationDistribution.