Multiscale Model of Colorectal Cancer Using the Cellular Potts Framework

Abstract

Colorectal cancer (CRC) is one of the major causes of death in the developed world and forms a canonical example of tumorigenesis. CRC arises from a string of mutations of individual cells in the colorectal crypt, making it particularly suited for multiscale multicellular modeling, where mutations of individual cells can be clearly represented and their effects readily tracked. In this paper, we present a multicellular model of the onset of colorectal cancer, utilizing the cellular Potts model (CPM). We use the model to investigate how, through the modification of their mechanical properties, mutant cells colonize the crypt. Moreover, we study the influence of mutations on the shape of cells in the crypt, suggesting possible cell- and tissue-level indicators for identifying early-stage cancerous crypts. Crucially, we discuss the effect that the motility parameters of the model (key factors in the behavior of the CPM) have on the distribution of cells within a homeostatic crypt, resulting in an optimal parameter regime that accurately reflects biological assumptions. In summary, the key results of this paper are 1) how to couple the CPM with processes occurring on other spatial scales, using the example of the crypt to motivate suitable motility parameters; 2) modeling mutant cells with the CPM; 3) and investigating how mutations influence the shape of cells in the crypt.

Introduction

Colorectal cancer (CRC) arises from the buildup of mutations in the cells of the colonic epithelium that lines the walls of the gut. In particular, these mutations occur in structures known as colonic crypts, a series of test-tube-shaped invaginations in the wall of the large intestine. These mutations cause cells to proliferate without the usual cell cycle regulation and to alter their mechanical properties, eventually leading to polyp-type formations where all structure of the epithelium is lost. ¹

Much work has been done to try to understand the processes underlying the onset of CRC, both experimental and theoretical. ² Early theoretical work mainly focused on the development and use of continuum models of the crypt, and CRC, and consisted primarily of compartmental models.^{3, 4} However, some of the first modeling work in the crypt was actually undertaken using a two-dimensional cellular automata (CA) model to investigate the distribution of cell types within the healthy crypt. ⁵ More recently, multiscale multicellular models, the spiritual descendants of the CA models, have been used to study the crypt and the onset of CRC. These approaches model the interactions of individual cells and study the resulting ensemble behaviors. They are the ideal modeling approach when you are looking to consider the effect of mutant cells, as you are able to adjust the properties of individual cells and see the effect on the tissue as a whole.

A multiscale model of the crypt based on the off lattice cell center was presented in Van Leeuwen et al. ⁶ This work coupled a previous model of the crypt ⁷ (which extended the earlier CA model but used force balance, instead of probabilities, to define cell-cell interactions) with a model of the cell cycle ⁸ developed for crypt epithelial cells. This model of a healthy crypt has since been extended to consider the effect of geometry on cell distribution, ⁹ to represent deformable substrates, ¹⁰ and to test hypotheses about active migration within the crypt. ¹¹ The model has also been used to study the onset of CRC, by including cells with mutant phenotypes, and has given insight into how mutant cells take over the crypt.^{12, 13} Recently, additional three-dimensional multicellular models have been developed that consider cells of multiple types (absorptive and secretory), where the cell type is determined by the internal state of the cell, which is updated at each timestep according to the state of neighboring cells, allowing the composition of the epithelial layer to be investigated. ¹⁴

Rather than considering force balance, a popular method for representing interacting cells is the cellular Potts model (CPM), which is derived from the Ising model of ferromagnetics and the large-q Potts model. ¹⁵ The first use of the CPM was by Graner and Glazier ¹⁶ and by Glazier and Graner, ¹⁷ where they used it to represent the aggregation and sorting of cells in vitro. The sorting was based on interaction energies, and this sorting property is fundamental to the CPM. Since its introduction, the CPM has been used to model many other biological systems, for example, morphogenesis,^{18, 19} blood clot formation, ²⁰ and vascular tumor growth and invasion.^{21, 22} because of its popularity there are many software packages that can run simulations using the CPM, two of the most established of these being CompuCell3D ²³ and Morpheus. ²⁴ Additionally, due to the lattice-based structure of the model, it is possible to derive continuum representations of the system, allowing larger parameter regimes to be investigated. ²⁵

Wong et al. ²⁶ presented the first model of the colonic epithelium using the CPM. This work focused on the pedigree model for cell proliferation. The pedigree model has an explicit population of immortal stem cells at the base of the crypt; these cells divide asymmetrically to produce an additional transit amplifying cell, which will subsequently divide a specified number of generations before differentiating. This pedigree model was compared to a positional-based model (where proliferative type is determined by external signaling factors) in Van Leeuwen et al. ⁶ which was shown to be unable to reproduce the colonizing of crypts seen experimentally. More recently, the random Boolean network models for cell-type specification have been coupled to a CPM to represent the spatial structure of the crypt. ²⁷ This model includes cell mutations through the evolution of hyperproliferative cell phenotypes in the random networks, but does not change the mechanical properties of cells. For more details of the history of experiments and modeling of CRC, see the reviews by Kershaw et al. ² and De Matteis et al. ²⁸

The key biological aim of this work is to investigate whether mutant cells can colonize the crypt by altering their mechanical properties. Moreover, we will study how different these mutant cells need to be (compared to healthy cells) in order to become persistent in the crypt. In addition, we aim to demonstrate how to include mutant cell phenotypes (with altered mechanical properties) in to the cellular Potts framework.

The remainder of this paper is structured as follows. First we introduce the mathematics underlying the CPM and discuss how to couple this model to temporal models of subcellular processes such as those for the cell cycle. The result of this section is a multicellular model of a healthy crypt based on the CPM framework which is subsequently investigated. We then discuss a method for implementing mutant cell phenotypes in the framework and present a model for the onset of CRC. We conclude by discussing future avenues for this work.

A CPM Wnt-Dependent Model of a Healthy Crypt

We now detail the elements of our multicellular model of the crypt. The model is similar to that used in Osborne et al. ¹² ; however, here physical interactions between cells are modeled using the CPM. As in Osborne et al. ¹² , the cell cycle is represented stochastically with proliferation and cell cycle duration controlled by the level of an external stimulus, Wnt. This stimulus can affect different cell types in different ways, allowing mutant cells to be represented. Moreover, due to the multiscale construction of the model, a more complicated cell cycle model (like the one in Van Leeuwen et al. ⁶ ) or other subcellular processes (as reviewed in Kershaw et al. ² ) could be included.

A CPM for cellular interaction

Here we use the CPM to represent the physical interactions between cells in the crypt. The traditional CPM has been used for many biological processes, including cell sorting, ¹⁶ morphogenesis, ¹⁸ and vascular tumor growth. ²¹ The CPM is described in detail in Izaguirre et al. ¹⁸ , among others, so here we just give an outline of the CPM and detail the changes to the general model, used previously, in order to represent the crypt and the onset of CRC. In the CPM, the spatial domain is discretized in to a finite collection of lattice sites, as shown in Figure 1A. Here, these sites form a regular lattice but they could be a set of arbitrarily connected sites of different shapes. Each cell is represented by a collection of such sites, and each site is contained in at most one cell. The shape of the cells is evolved by minimizing the energy of the system. This is done in a stepwise manner through evaluation of a Hamiltonian, H, (which represents the energy of the various cell-cell interactions) and metropolis sampling. If an appropriate Hamiltonian is given, then the lattice sites associated with each cell are evolved at each Monte Carlo step (MCS, a component of the metropolis sampling algorithm representing a sampling of the domain) using the following algorithm.

Figure 1

Potts model and the multiscale coupling. (A) Representation of cells in the cellular Potts model. Cells comprise a collection of lattice sites (here a regular square lattice) and evolve over time by selecting neighboring (or contained) lattice sites to be included in (or removed from) the cell. Cell boundaries are given by bold lines, and the central cell is shaded in gray. (B) Schematic drawing of the components of the model: cellular Potts model (CPM) for cell movement; cell cycle model (CCM) for control of cell proliferation; and Wnt stimulus for controlling the cell cycle. Arrows indicate coupling between components.

where the neighborhood of a lattice site is defined by the Moore neighborhood (ie, the eight sites that surround it, with appropriate adjustments for periodicity on the left and right domain boundaries). The parameter T (known as the “temperature”) is a residual parameter from the origins of the model in ferromagnetism and, while it has no specific biological significance here, it can be thought of as a motility or fluctuation parameter. ¹⁵ As T is increased, cells will move more freely, but the likelihood of the cell (defined as the set of lattice sites comprising it) becoming non-connected also increases. Therefore, a key element of any CPM simulation is choosing the value of the parameter T in order to balance it with the other model parameters in order to get realistic results. In a wider sense, T (specifically in relation to the magnitude of the terms in the Hamiltonian) is analogous to the viscosity or drag parameters η and μ found in cell center and cell vertex models, respectively, ¹² and we will use this property later to represent mutant cell phenotypes. The number of samples per MCS, N _S , how this MCS relates to an actual timestep dt, and the “temperature” T will dictate the dynamical properties of the system. For clarity, here we let N _S = N, where N is the total number of lattice sites, so we are aiming to sample each site at every MCS in an averaged sense, and we couple the CPM to the cell cycle by choosing

$$1\text{MCS}=dt.$$

(1)

This choice is arbitrary, as we are just specifying that one sweep of the domain corresponds to the timestep, dt, and we can subsequently vary dt to change the coupling scales. Therefore we use the terms MCS and timestep interchangeably, using whichever is more appropriate.

The dependence of the CPM, and therefore the simulations, on these key model parameters is an important property of the CPM and allows multiple behaviors to be considered and also allows simulations to be tailored to the system being investigated and the questions being asked. A detailed investigation in to this dependence will be undertaken for our system of the colorectal crypt in the next section.

The Hamiltonian is central to all CPM simulations and specifies how the cells evolve as the simulation progresses. It is comprised of multiple components, which can be included in the model to represent different biological properties or phenomena depending on the problem you are investigating. This is analogous to the idea of including forces in off-lattice models.^10,12,14,29 The most commonly used Hamiltonian components are cell-cell adhesion and area (volume in three dimensions) constraint. ¹⁷

In this model, we use

$$H=H_{\text{Adhesion}}+H_{\text{Area}},$$

(2) where H_Adhesion and H_Area are the adhesion energy and area constraint contributions, respectively. The adhesion energy contribution is given by

$$H_{\text{Adhesion}}=\sum _{(i,j)\text{neighbors}}\left(1-{\delta }_{\sigma (i),\sigma (j)}\right)J\left(\tau (i),\tau (j)\right),$$

(3) where σ(k) is the cell containing lattice site k, τ(K) is the type of cell σ(k), and J represents the interaction energies between cells and captures the effect of adhesion molecules on the surface of cells. This contribution represents a cell's propensity to remain as one solid unit rather than a set of disconnected regions and is common to all CPM simulations. Different adhesion energies J can be chosen to represent cell aggregation and sorting (this was one of the first applications of Potts type models in biology). ¹⁷ Here we choose to have cells of the same type (either healthy or mutant) have a stronger adhesion, J _nn = 0.1, than those of different types, J _nm = 0.2 (note that a smaller energy contribution implies a stronger adhesion); these choices are inspired by the values used in Wong et al. ²⁶

We include an area constraint in the Hamiltonian, and this leads to cells tending to a set target area A_Target(i, t), which can depend on internal properties of the cell allowing detailed models of cell growth. The contribution to the Hamiltonian is given by

$$H_{\text{Area}}=\underset{Cells}{\sum \alpha }{\left(A\left(i,t\right)-A_{\text{Target}}\left(i,t\right)\right)}^2,$$

(4) where A(i, t) is the area (or volume in three dimensions) of cell i at time t, and A_Target (i, t) is the corresponding target area (volume), α is a parameter that influences how fast cells react to the constraint. ¹⁷ The balance between α, the interaction energies J, and the other dynamic parameter Tinfluences the speed at which the cells move and rearrange, (dt also affects the speed but this contribution can be isolated by looking at the ratio between T and dt.) The larger the value of α, the more quickly the cells area will tend to the target area. Additional components can be added to the Hamiltonian to represent many other cellular processes, for example, a surface area constraint, Chemotaxis, or haptotaxis. ¹⁸

Crypt domain and multicellular coupling

Following Osborne et al. ¹² , we consider the crypt to be an unwrapped cylinder. Here we restrict cells to lie on the periodic regular lattice of width 50 and height 110 lattice sites. Once the center of a cell (defined by its center of mass) is above the height of 100 lattice sites, then it is assumed to be sloughed from the crypt and is removed from the simulation.

Cell proliferation is controlled by an imposed gradient of a signaling molecule Wnt (shown in Fig. 3). Cells of different types respond to the gradient in different ways. Normal cells can proliferate only in the lower third of the crypt where there is a plentiful supply of Wnt. While proliferating, each cell has a model of the cell cycle, and here we draw the duration of the whole cell cycle from a normal distribution with mean T_CCD and unit standard deviation. ¹² This duration includes all phases of the cell cycle. When a cell gets to the prespecified age (if there is still enough Wnt present, ie, it is in the lower third of the crypt), it divides by evenly distributing its lattice sites into two daughter cells (if there are an odd number of lattice sites, then the extra site is distributed at random) and each daughter cell draws its own cell cycle duration from the same normal distribution as the parent cell. A schematic of the components of the model and their coupling is given in Figure IB. This is a simplified approach to modeling the cell cycle and is used here as we wish to keep the model as simple as possible in order to isolate the influence of the mechanical properties of cells. In future work we could use a more complicated cell cycle model that is fitted to experimental data.

Figure 2

Cellular Potts model of a healthy crypt. Illustrative snapshots of the crypt simulation at t = 100 hours with: (A) sensibly chosen “motility parameters” (dt = 0.01, T = 0.1) and the crypt is in dynamic homestasis, and (B) poorly chosen parameters (dt or T too high, here dt = 0.1, T = 0.1, but similar behavior is seen for increased T). Cells proliferate (yellow cells) when there is a plentiful supply of Wnt stimulus (which linearly decreases up the crypt as shown on the left). Cells no longer proliferate (pink cells) when the Wnt value is below a given threshold (here this is a third up the crypt). Videos of the simulations for (A) and (B) can be seen in Supplementary Movies 1 and 2, respectively.

Figure 3

Comparing the CPM to existing models for the crypt. (A) Variation of the averaged compression, (A_Target – A)/A_Target with height up the crypt y (B) Variation of the average vertical velocity v _y with height up the crypt y Results are compared for the model presented here with dt ranging from 0.00125 to 0.04 in factors of 2 (solid lines and the arrow shows the direction of increasing dt) to the cell center model presented in Osborne et al. ¹² (dashed line) and to experimental results from Kaur and Potten ³¹ (dot dashed line). (C) and (D) Dependence of the cell number at dynamic equilibrium as dt varies: all cells are considered in (C), and proliferating and differentiated cells are shown separately in (D) (in red and blue, respectively). For cell numbers, one standard deviation to the measurements is shown by the shaded regions. Results are averaged from making observations every hour (after a 200-hour burn-in) in a 2,200-hour simulation for each model and parameter set.

Note that cell division can occur only between MCS and not during them, as each MCS represents the rearrangement of cells during one timestep. Also note that, if a cell is only occupying one lattice site when the subcellular machinery dictates the cell should divide (as could be the case for poorly chosen parameters or an unrealistic proliferation model), then the cell cannot divide and the simulation is halted.

In order to represent cell growth when a cell divides, its daughter cells have a target area of half the mature target area, and this increases linearly until it reaches the mature target area. In line with previous models of the crypt, we increase this over the first hour of the cell cycle.^{6, 12, 13}

Traditionally, within the CPM framework, cell division is usually modeled by constantly increasing the target area of individual cells, and once the actual cell area reaches a predefined threshold, it divides. ¹⁸ The approach taken here is more in line with established multiscale multicellular frameworks and allows cells to divide based on internal cues, permitting a full coupling between cellular and subcellular levels.

Parameters, initial conditions, and simulation

As specified above, all simulations presented here are on a domain of 110 lattice sites tall and 50 lattice sites wide, with periodic boundaries applied on the left and right, and cells are removed from the simulation when their center is more than 100 lattice sites from the base of the crypt. Cells have a mature target area of 25 lattice sites, and so the crypt is approximately 10 × 20 cell diameters (CDs) in size (in agreement with Osborne et al. ¹² ). All model parameters are given in Table 1.

Table 1

Parameters in the model.

N x , N y	Width and height of lattice	50,110 lattice sites
H Slough	Height at which cells are removed	100 lattice sites
T CCD	Mean cell cycle duration	16 hours
α	Deformation energy parameter	0.1
A Target	Mature cell target area	25 lattice sites
J hh , J mm	Adhesion energy - same cell types	0.1
J hm , J mv , J hv	Adhesion energy – different cell types	0.2
T	“Temperature”	0.1*
dt	Timestep	0.01* hours
v	Mutant cell “drag” multiplier	1–10

Note: “Default values – these are varied in Figure 3.

Simulations are initialized by creating a regular array of 10 × 20 square cells, each containing 25 lattice sites. Cell ages and cell cycle durations are drawn from a normal distribution with mean T_CCD and unit standard deviation. All cells are initially healthy and are specified as proliferative or differentiated based upon the Wnt level (ie, position up the crypt).

All simulations are undertaken using the open source Chaste (Cancer, Heart, and Soft Tissue Environment) framework.^{29, 30} Further details on Chaste are available at http://www.cs.ox.ac.uk/chaste. All the codes used to run the simulations presented here are released under an open source BSD license. The code and tutorials on how to use it can be found at https://chaste.cs.ox.ac.uk/trac/wiki/PaperTutorials/PottsCrypt2015. The time taken to run simulations depends on the size of the domain, the number of cells, and the timestep dt. As an illustration, the simulation presented in Figure 2A took less than 2 minutes to run on a single-core laptop computer (note this is for 100 hours of crypt time). The same amount of time is taken for each simulation run in the mutation examples in later figures.

Simulations of a CPM Wnt-Dependent Model of a Healthy Crypt

Before looking at the effect of mutant cell phenotypes on the development of the crypt, we investigate the dynamic equilibrium of the crypt model with only healthy cells. In particular, we investigate the model's dependence on the dynamic parameters resulting from the use of the CPM: T and dt.

Figure 2 presents results from two realizations of the model introduced in the previous section. In Figure 2A we present a snapshot at dynamic steady state (t = 100 hours) of a simulation using the default parameter set (dt = 0.01 and T = 0.1, as given in Table 1), and in Figure 2B we present a snapshot of the same simulation with a larger timestep dt = 0.1 at t = 100 hours. Videos of these simulations are given in Supplementary Movies 1 and 2, respectively. For the default parameter set, we see that the cells are all approximately the same size and shape, cells toward the base of the crypt are proliferating (where the Wnt level is sufficient, shown by yellow cells), and subsequently migrate up the crypt. With an increased timestep, Figure 2B, we see that cells are no longer of the same size and they have become fragmented, especially toward the base of the crypt. This is due to the restricted motility of cells (as can be seen in Fig. 3).

A comparison of how the cell compression and vertical velocity vary as we move up a healthy crypt is given in Figure 3. In Figures 3A and 3B, the dependence of compression defined as $\left(A_{\text{Target}}-A\right)/A_{\text{Target}}$ and the vertical velocity, v _y , with height up the crypt, are compared for increasing values of dt (from 0.00125 to 0.04 in factors of 2, solid lines). These results are compared with results for analogous simulations using a cell-center-based off-lattice model from Osborne et al. ¹² (dashed line) and to the experimental data from Kaur and Potten ³¹ (dot dashed line). Figures 3C and 3D show how increasing the timestep affects the number of cells in the crypt.

It is clear from Figures 3A and 3B that with larger values of dt the cells are more compressed, resulting in a larger number of cells in the crypt (this is because having a smaller timestep allows cells more MCSs to relax to their target area between cells dividing). Additionally, from Figure 3D, it is clear that the number of proliferating cells is affected more than the number of differentiated cells, which is due to the compression being greatest in the base of the crypt where the proliferating cells are found. However, interestingly, the overall cell velocity is relatively unaffected by dt, and is in excellent agreement with that for published off-lattice simulations ¹² and in good agreement with experimental data. ³¹ In future work, a better fit to experimental data could be achieved by modifying the other parameters of the model.

In order to fully investigate the way the parameters dt and T affect the simulation, in Figure 4 we present how the number of cells in crypt deviates from the expected value of 210 (taken from Fig. 3D, for dt = 0.01) at time t = 50 hours and also how the number of cells changes between t = 25 hours and t = 50 hours, for varying T and dt. From the bottom contour plot, we see that for large values of dt the crypt is not in homeostasis and the number of cells is increasing; this is backed up by the compressed cells in the red and blue snapshots. From the top contour plot we can see that, for small values of T, there is a larger number of cells in the crypt than expected and this gets worse as we increase dt (also shown in the yellow and red snapshots). In addition, a large value of T leads to a fragmentation of cells (shown in the blue and green snapshots); moreover, if we also have a small value of dt, cells are removed as a result of excessive sampling and movement, and this leads to a decrease in cell number (shown in the green snapshot). Taking all of these behaviors into account, we choose a default parameter set dt = 0.01 hours and T = 0.1. This is shown by the black snapshot in Figure 4, and these parameters are used in the remainder of this work.

Figure 4

Selecting dynamic parameters for the crypt. Center top: deviation of cell number at t = 50 from the expected value of 210 (from Fig. 3D, for dt = 0.01) as dt and T are varied. Center bottom: ratio of the numbers of cells at t = 50 and t = 25. Numbers are averaged from 100 simulations for each parameter set. Left and right: illustrative examples of crypt behavior at t = 25 are shown for the four extremes (red, yellow, blue, and green) and for the optimal parameter set dt = 0.01 and T = 0.1 (black).

Mutations in the Crypt

Following Osborne et al. ¹² , we now perturb the homeostatic healthy crypt by including mutant cells. Mutant cells are considered to have gained independence from proliferative signaling, and have altered their mechanical properties to promote invasion. ³² In fact, there is experimental evidence that Wnt-independent mutant cells (which have an inactivation of Ape) have stronger cell-stroma adhesion than normal Wnt-dependent cells. ³³

Modeling mutant cells

When using off-lattice models for cell interactions, the increased adhesion often observed in mutant cells (in order to promote invasion) is represented by increasing the resistance to motion (“drag”) experienced by mutant cells. This is possible because in off-lattice models the movement of cells is determined by a force balance and therefore increasing the “drag” a cell experiences will represent increased adhesion to a substrate (see Osborne et al. ¹² for details). The CPM framework utilized here does not use force balance to define cell movement, so there is no such thing as a “drag” parameter in the model as it stands. However, the reduced motility of specific cells in the CPM can be achieved by modifying the success probability of site swaps involving mutant cells. This can be achieved by rejecting a specified proportion of all moves involving mutant cells, and is implemented by adding the following term to the Hamiltonian, Equation (2):

$$H_{\text{Mutant}}=\{\begin{array}{cc}\infty & \text{If}\tau (i)\hspace{0.17em}\text{is}\hspace{0.17em}\text{mutant}\hspace{0.17em}\text{and}\hspace{0.17em}p<\frac{1}{\mu },\\ 0& otherwise,\end{array}$$

(5) where p is drawn from a uniform U(0,1) distribution, μ is the mutant cell “drag” ratio (which defines how many moves to reject), and τ(i) is the type of the cell containing lattice site i. In addition to having an increased “adhesion” (modeled through the μ parameter), mutant cells proliferate independently of the external Wnt stimulus level (ie, they proliferate throughout the crypt).

The effect of mutant cells

In order to show how our implementation of mutant cells within the CPM framework behaves, we again follow Osborne et al. ¹² , and place a patch of mutant cells with varying “drag” ratios in the crypt and track the size of the patch over time. Simulations of a healthy crypt as in the previous section are run until they are in a dynamic steady state (ie, for 50 hours), a patch of cells is then selected, all cells contained in the patch are converted to mutant cells, and we track the size and shape of the patch over time. We define the patch by specifying its center and radius, and all cells who's centers are within this radius from the center of the patch are defined to be mutant cells. In this work, we consider a radius of 2 CDs and investigate a “high patch”, placed at (5 CDs, 4 CDs) and a “low patch” placed at (5 CDs, 12 CDs); these can be seen in the t = 0 simulation snapshots in Figure 5.

Figure 5

Invasion of crypts by mutant cells. (Center) Plots showing how the vertical height of the center of the lowest cell (solid line) and the highest cell (dashed line) in the mutant patch changes over time as the mutant cell “drag” ratio varies. The initial patch is a circle of radius 2 CDs, centered at a height of 4 CDs (top) and 12 CDs (bottom). Results are plotted at times t = 0, 10, 20, 30, 40, and 50, and the arrow shows the direction of increasing t. Results are averaged from 100 simulations, and illustrative simulations are shown for the low (left) and high (right) initial mutant patches for “drag” ratios μ = 1,5, and 10, (top, middle, and bottom, respectively). Videos of these illustrative simulations for low and high initial mutant patch are given in Supplementary Movies 3 and 4, respectively.

In Figure 5, we track how the position of the mutant patch is affected as we increase the mutant “drag” parameter, μ, from 1 to 10. We do this by tracking the position of the base (and top) of the mutant patch at times t = 0, 10, 20, 30, 40, and 50 hours after the addition of mutant cells, for each μ. Illustrative simulation snapshots are also included for low (μ = 1), medium (μ = 5), and high (μ = 10) mutant “drag” ratios for both the low and high mutant patch. Videos of these illustrative simulations are given in Supplementary Movies 3 and 4, respectively.

We see from the central plots of Figure 5 that, as in the off-lattice case, ¹² as the “drag” μ of the mutant cells increases, the base (and top) of the patch are lower down the crypt. This can also be seen in the illustrative simulation snapshots in Figure 5. Additionally, we see that in the case of a low patch, for a mutant “drag” ratio of greater than about μ = 6, the mutant cells are migrating downward. This downward invasion is more pronounced in these CPM simulations than the off-lattice simulations in Osborne et al. ¹² (as also seen in Fig. 6), and may provide more evidence for the top-down mechanism for morphogenesis. ¹

Figure 6.

Comparison of CPM with center-based model. How the vertical height of the center of the lowest cell in the mutant patch changes overtime as the mutant cell “drag” ratio varies, with initial mutations centered at a height of 4 CDs, in (A) our CPM-based simulations, from Figure 5, presented on different axes for comparison, and (B) in a cell-center-based representation. Reproduced with permission from Osborne et al. ¹² Results are plotted at times t = 0, 10, 20, 30, 40, and 50, and the arrow shows the direction of increasing t. The critical values of the “drag” parameter, where cells begin to invade down the crypt, are marked by ${\mu }_P^{*}\approx 5.75$ and ${\mu }_C^{*}\approx 15$ in (A) and (B), respectively.

In order to see how adjusting the Hamiltonian, as in Equation (5), compares to increasing the drag coefficient in off-lattice models, we present a comparison with some of the results from Osborne et al. ¹² in Figure 6. As discussed above, the qualitative behavior is the same as we increase the mutant “drag” parameter, and we see that increasing the parameter μ has the same effect as increasing the drag coefficient in center-based simulations. By matching the value of μ for which cells invade down the crypt, below the initial height, we can see that (for the other parameters as defined) a value of the “drag” parameter of ${\mu }_P^{*}\approx 5.75$ in the CPM is equivalent to a drag coefficient of ${\mu }_C^{*}\approx 15$ in cell-center-based simulations. Therefore we can relate the two parameters by the relationship ${\mu }_P\approx {0.4}_{{\mu }_C}$ , where μ _P is the “drag” parameter for the CPM in this paper and μ _C is the drag coefficient from cell-center models. ¹²

One key advantage of the CPM over the traditional cell-center off-lattice model used previously is that the cell shape (ie, area and perimeter) is a driver of cell movement, and we are able to easily track how cell shape within the crypt is affected by the presence of mutant cells. In Figure 7 we present how, due to the presence of mutant cells, the average cell area, A, perimeter, P, and circularity, defined as $C=\frac{4\pi A}{P^2}$ vary as you move up the crypt. (Note that, due to the fact that we are using a regular lattice, the maximum circularity is π/4 ≍ 0.7854 rather than the usual maximum of 1). These averages are calculated between 40 and 50 hours after the mutant patch is initialized (initially centered on an height of 4 CDs, ie, the low patch) for both low (μ = 1, in black) and high (μ = 10, in blue) mutant “drag” ratios. In addition, on each plot, the proportion of mutant cells in the crypt is also presented for low and high mutant “drag” ratios; these are given as histograms in black and blue, respectively. Each plot is created by averaging the relevant size and shape properties of cells whose centers lie in each height window (of width 1 cell diameter) between 40 and 50 hours after the mutant patch is initialized. Furthermore, these data are also averaged over a collection of 100 realizations (utilizing different random number seeds) in order to show the typical behavior of the system.

Figure 7.

Changes in cell shapes due the presence of mutant cells in the crypt, showing how the average area (A), perimeter (B), and circularity (C) of cells is affected by presence of mutations. Solid lines give the average for all cells, and dashed lines give the average for healthy cells only. Mutant cells with 10 times the “drag” of healthy cells (μ = 10) are represented in blue, and mutant cells having the same “drag” as healthy cells (μ = 1) are represented in black. In each plot, the average distribution is calculated from the same 100 simulations from Figure 5 between t = 40 and t = 50 hours from the addition of mutant cells. For comparison, the corresponding average shape distributions are also give for a crypt with only healthy cells (dotted line). The distribution of mutant cells up the crypt can be seen at the bottom of each figure, with black and blue representing μ = 1 and μ = 10, respectively.

From Figure 7A, we see again that, in the case of a crypt of only healthy cells (dotted line) there is an increase in the level of compression as you go down the crypt. This was shown to depend on the dynamic parameters of the CPM in Figure 3A, and here the typical values (of T = 0.1 and dt = 0.01) are used. Although the cells are more compressed toward the base of the crypt, the cell perimeter remains relatively constant as shown by the dotted line in Figure 7B. This leads to slight decrease in circularity as you go down the crypt (dotted line in Fig. 7C). The addition of mutant cells with low “drag” ratio does not change the shape of the area, perimeter, or circularity distribution; it just results in slightly smaller cells (the solid and dashed black lines in Figs. 7A-C). This is because mutant cells are predominantly being washed out of the crypt by time t = 40 hours (as seen by the black histogram in Figs. 7A-C). Before being washed out mutant cells are still proliferating in the upper areas of the crypt, increasing the number of cells and lowering the average area slightly.

When including mutant cells with a high “drag” ratio, the mutant cells remain lower in the crypt (seen by the blue histograms in Figs. 7A-C and the simulation snapshots in Fig. 5) and the effect on cell shape is much more marked. We see that the level of cell compression toward the base of the crypt is much larger for all cells, but it is actually slightly less for healthy cells (solid and dashed blue lines, respectively, in Fig. 7A). (Note that the average compression decreases again at the base of the crypt compared to that of healthy cells. This is because there are fewer mutant cells there. In fact, there are none at the very base, as seen in the blue histogram, and therefore the contribution from mutant cells is less.) This is due to the fact that a by-product of cells having an increased “drag” in the CPM is that it will take longer for cells to grow and therefore they may divide before reaching full mature size, resulting in slightly more compressed cells. However, here the perimeter for healthy cells is smaller than the average for all cells (dashed and solid blue lines, respectively, in Fig. 7B). These two observations combine together when looking at the circularity, and we see that healthy cells have a relatively constant circularity in the crypt, whereas we see a decrease in the average circularity from all cells where mutant cells are present (dashed and solid blue lines, respectively, in Fig. 7C). The result of this is that, when a population of mutant cells is fixated in the crypt (as is the case for μ = 10), the circularity of the mutant cells is lower, suggesting that the cells are no longer approximately spherical in shape and are more irregular. This can be seen in Figure 8, where we present magnifications of the cell boundaries in a region of healthy cells and a region of mutant cells in an illustrative simulation with a large mutant “drag” ratio (μ = 10). We see that the cell boundaries are much more irregular in the mutant region, explaining the lower circularity when mutant cells are present.

Figure 8.

Loss of mutant cell shape regularity. Snapshot of illustrative simulation for a low initial mutant patch with μ = 10 at t = 50 hours, from Figure 5. Cell shapes magnified for a healthy region (blue) and mutant region (red). The color of cells has been removed to clarify the cell boundaries.

Discussion

In this paper we have presented a model of mutations in the colorectal crypt. This model builds upon previous work,^{6, 12, 13} by utilizing a different cell-based modeling paradigm, the CPM. The presented crypt model is the first, to the authors knowledge, to represent mutations, with different mechanical properties, in the colorectal crypt using a CPM. The implementation of mutations here represents a mechanism by which mutations could be represented in the CPM for other biological systems. The use of a CPM has allowed us to not only verify the results of earlier publications using off-lattice models ¹² but has also enabled us to undertake a detailed investigation of the effect that mutant cells have on the shape of cells within the crypt. While parts of this would have been possible using the vertex-based model presented previously, ¹² the level of cell deformations would of be much less observable.

This work could be extended in two key ways. First, further investigation into the workings of the colorectal crypt could be undertaken; in particular, other types of mutations could be investigated. In addition, the crypt model could be extended to include more realistic models of the cell cycle, as discussed in Kershaw et al. ² or using a different representation for the cell cycle and cell fate specification, for example, a random Boolean network ²⁷ or a state-based model. ¹⁴ Further models of cell interactions, for example, the subcellular element method (SEM), ³⁴ could be used to represent the mechanical components of the system. The second main avenue of future work would include using the implementation of the CPM, developed here, for other biological systems, for example, problems in development or avascular tumor growth. This would involve introducing other components into to the Hamiltonian to represent other biological constraints, for example, surface area constraints and Chemotaxis,^{18, 21} but this would be straightforward given the modular framework developed and used here.

In summary, the work presented here gives new insight into the effect of mutations on the structure of the colorectal crypt and, moreover, presents a framework for including cell phenotypes with varied mechanical properties within the CPM.

Author Contributions

Conceived and designed the experiments: JMO. Analyzed the data: JMO. Wrote the manuscript: JMO. The author reviewed and approved of the final manuscript.

Supplementary Materials

Supplementary Movie 1

Video of simulation of healthy crypt from t = 0 to t = 100 hours with sensible motility parameters. The end point of this movie is shown in Figure 2A. SupplementaryMoviel.mp4 or http://youtu.be/pmj_KOX1LIg.

Supplementary Movie 2

Video of simulation of healthy crypt from t = 0 to t = 100 hours with inappropriate motility parameters. The end point of this movie is shown in Figure 2B. SupplementaryMovie2.mp4 or http://youtu.be/gszAMMQAKUA.

Supplementary Movie 3

Video of a comparison between three simulations of a crypt with mutant cells from t = 0 to t = 50 hours. The video shows ‘drag’ ratios of μ = 1, 5, and 10 (from left to right) with the initial mutant blob placed low in the crypt. Snapshots from this movie are shown on the left of Figure 5. SupplementaryMovie3.mp4 or http://youtu.be/RKiNzip2ZPM.

Supplementary Movie 4

Video of a comparison between three simulations of a crypt with mutant cells from t = 0 to t = 50 hours. The video shows ‘drag’ ratios of = 1; 5; and 10 (from left to right) with the initial mutant blob placed high in the crypt. Snapshots from this movie are shown on the right of Figure 5. SupplementaryMovie4.mp4 or http://youtu.be/8qlLFpV7ugI.