Towards a multiscale model of colorectal cancer

Abstract

Colorectal cancer (CRC) is one of the best characterised cancers, with extensive data documenting the sequential gene mutations that underlie its development. Complementary datasets are also being generated describing changes in protein and RNA expression, tumour biology and clinical outcome. Both the quantity and the variety of information are inexorably increasing and there is now an accompanying need to integrate these highly disparate datasets. In this article we aim to explain why we believe that mathematical modelling represents a natural tool or language with which to integrate these data and, in so doing, to provide insight into CRC.

Key Words: Intestinal epithelium; Crypt fission; APC mutations; Mathematical modelling; Stem cell niche; Wnt signalling

INTRODUCTION

Colorectal cancers (CRCs) originate from the rapidly renewing epithelium that covers the luminal surface of the bowel and lines the colonic crypts (Figure 1). At the bottom of these crypts, stem cells proliferate continuously producing amplifying transit cells that divide rapidly several times before differentiating into the various cell types that constitute the epithelium (e.g. absorptive columnar cells, mucinous Goblet cells, and neuroendocrine cells). The production of new cells is synchronised with migration of the epithelium upwards along the crypt axis. Upon reaching the orifice, cells undergo apoptosis and are shed into the lumen. Under normal conditions, these cellular processes are tightly regulated by biochemical and biomechanical signals (e.g. Wnt factors, cell-cell contacts and growth factors)[1,2]. CRC arises when genetic alterations, for example in the Wnt pathway, disrupt normal crypt dynamics. Proliferating cells are then no longer confined to the crypt’s base and the associated proliferative excess generates biomechanical stress on the crypt structure, which may deform in order to accommodate the newborn cells. The dysplastic cell population will expand further by invading neighbouring crypts and/or inducing crypt fission, leading to the formation of an adenoma[3,4].

Figure 1 Normal colonic mucosa. The normal colonic mucosa consists of crypts (lined with epithelium) embedded within a connective tissue stroma. The mucosa rests on a layer of smooth muscle, called the muscularis mucosae, below, which is a loose connective tissue layer, called the submucosa. The epithelium consists of a number of different cell types including goblet cells (containing mucin), enterocytes and neuroendocrine cells (which secrete hormones). The crypts are enveloped by a sheath of myofibroblasts that are separated from the epithelium by basement membrane. The stroma contains a number of different cell types including macrophages, fibroblasts and lymphoid cells. The crypt is a dynamic structure that turns over every 3-6 d. The base of the crypts contains the stem cells that are thought to give rise to transit amplifying cells. The cells migrate up the crypt and undergo differentiation along the way. The goblet cells secrete mucin into the crypt and intestinal lumen and the surface epithelium usually contains few goblet cells. Upon reaching the surface, the cells undergo apoptosis and are shed into the lumen. The migration along the crypt by the epithelial cells is accompanied by a similar migration of myofibroblasts. The microenvironment is thought to change along the crypt axis and this is thought to give environmental cues for induction of maturation and inhibition of proliferation.

Cancers of the colon and rectum are the third most frequent malignancies in the United Kingdom, being responsible for 16 148 deaths in 2004[5]. Most CRCs occur sporadically, with only 5% attributable to hereditary syndromes[6], such as familial adenomatous polyposis (FAP). The first detailed genetic model for the development of CRC was published in 1990[7], with adenomatous polyposis coli (APC), K-ras, DCC and p53 as candidate cancer genes that undergo mutation in an ordered sequence. Due to its high socio-economic impact, CRC is studied extensively world-wide. While some of the mechanisms of CRC development have now been revealed, as more biochemical, genetic, histological and clinical information is made available, new layers of complexity are becoming apparent. We believe that CRC research has now reached a stage at which theoretical thinking and mathematical modeling have become essential to integrate and make sense of the biological information being generated and, more importantly, to generate new biological hypotheses that can then be tested in the laboratory. Among the benefits of modeling are: (1) mathematical models act as in silico tools with which to carry out and iterate virtual experiments that might be dismissed by experimentalists for being unethical, expensive, time consuming, or simply impossible; (2) modeling necessitates the statement of explicit hypotheses, a process which often enhances comprehension of the biological system; and (3) simulations can reveal hidden patterns and/or counter-intuitive mechanisms in complex systems. The use of mathematical modeling to tackle complex biological problems is becoming increasingly popular within the life sciences[8,9]. Hence, cancer modeling is a rapidly growing field within applied mathematics (critical reviews[10,11]) and, in particular, a number of models have been developed specifically to study issues associated with the normal gut and CRC (reviewed in Van Leeuwen et al[12]). Epidemiological models, for instance, have sought to describe CRC incidence in human populations[13,14]. Notably, such models have provided insight into the multistage nature of carcinogenesis and the importance of clonal expansion. Furthermore, discrete models characterising the position and behaviour of individual cells within an intestinal crypt have been formulated to advance our understanding of the mechanisms of cell migration and differentiation[15-17], while stochastic models have been used to investigate the accumulation of genetic (or epigenetic) alterations in epithelial cells[18-21].

Most theoretical studies of CRC, including those exemplified above, were designed to investigate specific biological questions and, consequently, they typically focus on a single time and length scale. In reality, the dynamics of normal tissue, colorectal carcinogenesis and cancer progression involves the integration of a large number of processes occurring over different levels of organisation; ranging from the fast enzymatic reactions taking place inside individual cells to the relatively slow dynamics of tumour growth and metastasis (Figure 2). In order to understand the ways in which subcellular events influence macroscopic tumour progression, it is necessary to develop models that incorporate multiple time and length scales. Such multiscale models are also needed to account for feedback within the system; for example, hypoxic cues, mechanical stress, immunological response and chemotherapeutic interventions all act on the macro- and cell-scales but have a profound influence on subcellular dynamics[22,23]. A natural consequence of this coupling across scales is that models cannot be built from the bottom-up by coarse-graining proteomic data, even if this were computationally feasible[24,25]. We propose that a multilevel approach be adopted, with (sub)models assembled at different scales before being coupled together so that more complex emergent behaviour can be explained.

Figure 2 Prototype multiscale model of the gut. Subcellular, cellular and tissue elements are represented in black, grey and white boxes, respectively. The schematic shows possible interactions and feed-back loops coupling the different levels of organisation and illustrates the complexity of the biological system. An additional challenge in building such multiscale models is coupling phenomena that act on very short time scales (e.g. the diffusion of metabolites and signalling molecules in the extracellular medium) with processes that occur on longer time scales (e.g. tumour growth).

The potential use of multiscale modeling to tackle complex biological problems is best exemplified by the sophisticated computational models of the heart developed by Noble and co-workers[26,27] where the simulation of normal and abnormal physiology, from the molecular level to the whole-organ, is well established. Only recently have these efforts been mirrored by similar projects focused on solid tumour growth. For example, Swanson and co-workers[28] have developed three-dimensional models that describe the spread of gliomas through the brain. While the mathematical model that they use is highly idealised (Fisher’s equation describes how the tumour cells migrate and proliferate), they have gone to considerable efforts to replicate the detailed anatomical structure of the brain and to calibrate their models against MRI scans. In particular, the random motility and proliferation rates of tumour cells are assumed to differ between grey and white matter. Their models show excellent agreement with clinical data and are being used to predict the response of patients to different treatment protocols. The limitations of the models described in Swanson et al[28] are that they do not explicitly account for nutrient supply to the expanding tumour mass and they neglect the multiscale nature of solid tumour growth. Several multiscale models of tumour growth that attempt to resolve these issues have been developed recently[22,29-31]. For example, Alarcón and co-workers have formulated a hybrid cellular automaton of vascular tumour growth that links submodels that describe processes occurring at the subcellular, cellular and tissue scales[22]. In more detail, progress through the cell cycle and the production of angiogenic factors are incorporated at the subcellular level and competition for resources are considered at the cellular level. The transport of nutrients and angiogenic factors, blood flow and vascular adaptation are included at the tissue scale. Coupling between the different submodels occurs in several ways. For example, oxygen levels, which are determined at the macroscale, influence subcellular processes such as progress through the cell cycle and the production of angiogenic factors. Equally, subcellular production of angiogenic factors modulates vascular adaptation at the tissue scale. Model simulations reveal that the tumour’s growth dynamics and, in particular, its morphology, are heavily influenced by the way in which the haematocrit distribution in the vasculature is modeled. For example, if it is heterogeneously distributed, then the tumour’s growth is highly irregular, being confined to well-oxygenated regions. By contrast, if the haematocrit is uniformly distributed then the tumour expands radially outwards as a compact mass.

Most existing multiscale models of solid tumour growth are generic in nature, in the sense that they are not specialised to a particular tissue or a specific type of cancer. We believe that it is now timely to develop mathematical models of this type to describe CRC. We aim in this article to explain the potentially valuable role of multiscale modeling in advancing our understanding CRC. We start by describing briefly some areas where modeling has generated insight into CRC (for a more exhaustive review, see Van Leeuwen et al[12]). These include signalling pathways, CRC initiation and crypt buckling. We conclude by explaining how we intend to integrate such models within a multiscale framework and, in so doing, enable experimentalists to explore, in a controlled in silico environment, the ways in which phenomena occurring on different space and time scales interact (Figure 2).

BIOCHEMICAL NETWORKS IN CRC

APC inactivation is believed to be an early, initiating event in CRC. Genetic alterations in the APC gene have been found in about 85% of all CRC cases and germline mutations are responsible for FAP[6,7]. Although APC is a multifunctional protein, it is generally accepted that APC’s tumour suppressor potential relies upon its ability to inhibit Wnt signalling[32]. This is best illustrated by the fact that activation of Wnt signalling by other means obviates the need for APC inactivation. Wnt signalling is essential for maintaining the stem-cell niche at the base of the crypt[33] and it has been proposed that, under normal physiological conditions, spatial gradients in extracellular factors acting on Wnt signalling, regulate cell proliferation, migration, differentiation and apoptosis along the vertical crypt axis[34]. This maintains the integrity of the normal crypt and allows renewal of the intestinal epithelium to occur in an ordered and controlled fashion. Inappropriate activation of Wnt signalling, however, results in loss of this normal control, leading to tumour development. The importance of Wnt signalling in carcinogenesis means that a model of this subcellular biochemical network should play an integral role in the development of any multiscale model of CRC.

Since the biochemical details and biological implications of the Wnt pathway have been reviewed elsewhere[32,35], here we summarise the key points. In the absence of an external Wnt stimulus, a destruction complex (including APC and axin) keeps cytoplasmic and nuclear levels of β-catenin low, thus allowing the repressor Groucho to bind preferentially to TCF/LEF transcription factors, thereby preventing the expression of Wnt target genes. In contrast, when extracellular Wnt factors bind to specific cell-surface receptors, a cascade of protein activation/inactivation events is triggered, leading eventually to the inhibition and destabilisation of the β-catenin destruction complex. Consequently, β-catenin protein accumulates in the cytoplasm and travels to the nucleus, whereby displacing Groucho from TCF/LEF proteins, it can induce the expression of target genes. Lee et al[36] have formulated a mathematical model of the Wnt pathway that describes normal protein turnover, step-wise assembly/dissociation of the destruction complex, binding and phosphorylation of β-catenin by the destruction complex, and assembly/dissociation of the transcription complex formed by β-catenin and TCF. The model provides quantitative expressions for dynamic changes in the level of 15 molecular components, based on parameter values either measured or estimated from experiments with Xenopus extracts. The combination of theoretical and experimental results revealed the importance of a novel regulatory loop involving APC-dependent degradation of axin. Hence, when this feedback mechanism is active, β-catenin levels are relatively resistant to reductions in active APC.

While Lee and co-workers[36] focused on the role of β-catenin in signalling, Van Leeuwen et al[37] have recently developed an alternative model of the Wnt pathway that accounts for its dual role in signalling and cadherin-mediated cell-cell adhesion. Instead of restricting their analysis to a limited set of parameter values, they used mathematical methods to gain insight into the general properties of the model system. In Figure 3, we exemplify this approach with a simple network involving only three components, i.e. A, B and AB. The reader can reduce the level of abstraction by considering the often presumed hypothesis that E-cadherin (A) can sequester β-catenin (B) from signalling. Both proteins are subject to normal protein turnover and to reversible binding to form a complex AB. The time-dependent change in the concentration of each molecular component is determined by its production and elimination rates. However, once enough time has elapsed, the compounds will reach a steady-state and no further changes in concentration will occur. For the given system, the final levels can easily be deduced from the kinetic equations. The final concentration of B, for example, is given by the ratio of its synthesis and degradation rates and does not, therefore, depend on any parameters associated with either A or the initial concentrations. This implies that the long-term, equilibrium level of cytoplasmic β-catenin is independent of E-cadherin and, consequently, changes in E-cadherin expression can only transiently affect the cytoplasmic β-catenin level. This example illustrates the potential value of adopting a theoretical approach; even in the absence of detailed parameter values, it is possible to extract meaningful conclusions and to generate new (experimentally-testable) hypotheses.

Figure 3 Example of a simple kinetic model. The proteins A and B undergo normal protein turnover (s = synthesis rates, d = degradation rates) and associate to form a protein complex, AB (k₊ = assembly rate, k_- = dissociation rate). The behaviour of the system is captured by the following system of ordinary differential equations: d [A]/dt = s_A - k₊ [A] [B] + k_- [AB] - d_A [A]; d [B]/dt = s_B - k₊ [A] [B] + k_- [AB] - d_B [B]; d [AB]/dt = k₊ [A] [B] - k_- [AB]; with initial conditions [A] (0) = A₀, [B] (0) = B₀ and [AB] (0) = 0. It can be shown that, eventually, the components stabilise at levels given by: [A] = s_A/d_A; [B] = s_B/d_B, and [AB] = s_A s_B/(d_A d_B K), with K = k_-/k₊ the dissociation constant.

The approach of Van Leeuwen et al[37] made it possible to compare two potential mechanisms of interaction between signalling, transcription and cell-cell adhesion. According to the first hypothesis (H.1), E-cadherin and TCF simply compete for β-catenin. Alternatively, under H.2, the Wnt signal induces the production of a ‘closed’ molecular form of β-catenin that binds preferentially to TCF. This might occur, for instance, via the activation of a kinase that phosphorylates β-catenin at Tyr-142. As might be expected, in both cases the model predicted an increase in the level of transcription complexes in response to Wnt. Surprisingly, however, the final equilibrium level of β-catenin-TCF complexes was independent of whether H.1 or H.2 was active and, consequently, it is not possible to discriminate between the two alternatives based on target gene expression alone. Under H.1, the model always predicted an increase in cell-cell adhesion in response to Wnt, whereas under H.2 this was not always the case: under H.2, the degree of cell-cell adhesion at equilibrium decreased as the rate of β-catenin Tyr-phosphorylation increased.

Van Leeuwen et al[37] also used their model to investigate the impact of two genetic alterations commonly observed in CRC; i.e. aberrant inactivation of APC and mutations that render β-catenin resistant to APC-mediated degradation. Simulations showed that in wild-type cells and single mutants, the level of transcription complexes depends on the strength of the Wnt gradient; by contrast, double mutants achieve the same (high) level of target gene expression and this value does not depend on the strength of the Wnt stimulus. As a consequence, the proportion of functional APC needed to maintain a normal phenotype increases with increasing extracellular Wnt levels. This result illustrates that the environment can have a substantial impact on malignant transformation and tumour phenotype.

STOCHASTIC MODELS FOR CRYPT DYNAMICS AND CRC-INITIATION

According to Knudson’s classic model for inactivation of tumour suppressors, two mutations (one in each copy of the gene) are required to completely abrogate APC’s anti-cancer activity. Given the extensive period of time required for a cancer to develop, stem cells have been widely regarded as the most likely targets for these mutations due to their apparent immortality. Komarova and Wang[18] used mathematical modeling to investigate whether neoplastic change could, instead, take place in the semi-differentiated, transit cell population. For this purpose, they considered an idealised crypt containing a single stem cell and investigated three complementary scenarios: (SS) both hits happen in stem cell; (ST) the first hit occurs in the stem cell and the second in a transit cell; and (TT) the first hit is in a transit cell and the second in one of its progeny. Scenarios SS and ST share an intermediate state in which the crypt is monoclonal, being filled with single APC mutants (APC^+/- cells). Challenging a priori intuition, the model revealed that the probability that the first double mutant appeared in the stem cell pool was negligible. The authors thus concluded that at least one of the hits occurs in the migrating, transit population and remarked that “a possible role of differentiated cells in cancer initiation cannot be discarded simply based on the fact that they are short-lived[18].”

In Komarova and Wang[18], the ‘immortal’ stem cell always divided asymmetrically, renewing itself and producing a transit cell. Their model can therefore easily be extended for crypts with multiple immortal stem cells, as such crypts can be seen as multiple single-stem cell crypts. Alternatively, Yatabe et al[21] considered ‘mortal’ stem cells that underwent asymmetric division, symmetric division (giving two stem cells), or cell death, with certain probabilities. To gain insight into the number and patterns of division of human stem cells, they monitored in silico the accumulation of random methylation changes in individual crypts. By comparing the predicted intra- and inter-crypt variability with experimental data on the methylation content of specific CpG islands, they concluded that “the complex epigenetic patterns were more consistent with a crypt niche model wherein multiple stem cells were present and replaced through periodic symmetric divisions[21]”. An important consequence is the concept of niche succession, which is, after a random waiting time, the crypt will be filled with the progeny from a single stem cell. More recently, it has been suggested that these ‘lineage bottlenecks’ enhance malignant transformation, because the progeny of an APC^+/- stem cell could take over the stem-cell pool, thereby increasing the number of candidates for the second hit. Van Leeuwen et al[12] performed a series of stochastic simulations to test this hypothesis. Their model was built upon the following basic assumptions: (1) stem cells divide synchronously; (2) for any given generation, two random sets of m stem cells divide symmetrically and die/differentiate, respectively; (3) newborn cells acquire an APC mutation with a certain probability. For a crypt with 64 stem cells, simulations were iterated for values of m between 0 and 16. As expected, the average niche succession time decreased as the number of symmetric divisions per generation increased. The average time to a first double mutant, however, appeared to increase with m. In particular, crypts with immortal stem cells only (i.e. m = 0), on average, developed cancer earlier than crypts with mortal stem cells (i.e. m > 0). As expressed by the authors, “the stem cell lineage bearing a single APC mutation has more chance of becoming extinct during the next niche succession cycle than of benefiting from the advantages of fixation in the crypt”. It is therefore proposed that symmetric stem cell division and niche succession protect stem cells against malignant transformation[12].

Lobachevsky and Radford[38] have recently proposed a different view of the crypt’s architecture in that their in silico stem cells cannot divide indefinitely. This assumption was motivated by experimental evidence supporting the occurrence of significant telomere-shortening in stem cells. The authors classified stem cells into two pools; i.e. deep, slowly proliferating stem cells that produce rapidly dividing proximate stem cells. The former always divide asymmetrically whereas the latter either divide symmetrically or die, as necessary, to maintain the homeostasis and integrity of the crypt. Furthermore, the probability of symmetric division or death depends on the cell’s age. Comparison of model predictions with experimental data suggested a replicative lifespan for stem cells of about 50-150 divisions, which differs markedly from the 5000 divisions estimated by other modellers. The new model was also able to reproduce the observed time-dependent appearance of partially and monoclonally mutant crypts in mice following treatment with ethyl nitrosourea.

SPATIALLY-STRUCTURED MODELS OF CRYPT-DYNAMICS

None of the models discussed thus far incorporate spatial effects. For example, the models for Wnt signalling view cells as ‘well-stirred buckets’ and, thus, consider only time-dependent changes in protein levels[36]. Equally, the carcinogenesis models record the accumulation of genetic or epigenetic alterations in a population of cells, without considering the relative positions of individual cells within an individual crypt[12,18,21]. Given that external cues to such subcellular phenomena are often spatially distributed (c.f. variations in extracellular Wnt levels along the crypt axis), it is natural to include spatial effects in order to describe accurately the dynamics of normal and pathological crypts.

Several computational models of normal crypt dynamics that account for spatial effects have been proposed. In these models the cylindrical crypt is usually cut open and rolled out to give a two-dimensional surface, which is computationally easier to simulate. For example, Loeffler and co-workers[15,17] have built a series of models in which cells are constrained to move within a rectangular grid. During their simulations, they monitored the progress of individual cells through the cell cycle. This enabled them to compare their results with tritiated thymidine labelling data showing how the proportion of proliferating cells changes with spatial position and over time. The model was used to compare eight alternative hypotheses concerning the positioning of newly-formed daughter cells. Of these, only two of the hypotheses were consistent with the experimental data. The first placed the daughter cell at the site of the oldest neighbour, whereas the second combined random insertion with inevitable final differentiation at a certain position along the crypt axis[15].

A major weakness of the 2D-grid approaches is that insertion of each newborn cell displaces a whole column of cells upwards, breaking many cell-cell contacts. This is biologically unrealistic, as epithelial cells are known to form tight cell-cell junctions. Meineke et al[16] resolved this problem by developing a lattice-free model in which cell movement was governed by the net effect of the repulsive and attractive forces acting on each cell, these forces being modeled by a network of springs connecting the centre of neighbouring cells (Figure 4A). An advantage of this approach is that the insertion of newborn cells causes local rearrangements only. To display the geometry, the nodes from the network of springs were used to generate a Voronoi tessellation (Figure 4B). The resulting model was then used to simulate the dynamics of a healthy crypt and shown to reproduce the experimentally observed polygonal shape of epithelial cells and the spatial distribution of proliferating cells. Although the model by Meineke et al[16] resolved the major weakness of the 2D-grid models, both types of model have several limitations in common. For instance, no biological mechanism regulates the removal of cells from the upper edge of the crypt and the model does not account for cell death occurring in cells other than fully differentiated cells. Furthermore, cell size is independent of growth processes, the time since the cell’s last division and extracellular conditions. Finally, the proliferative capacity is an intrinsic property of each cell and the progression of individual cells through the cell cycle is fixed at the time of their birth, with the duration of each phase being independent of any biochemical or biomechanical control (e.g. Wnt signals and cell-cell contacts). Integrating biologically realistic models of such subcellular processes and biochemical network dynamics with extracellular cues will be challenging and computationally intensive, with the approach depending on the correct coupling of the constituent submodels.

Figure 4 Schematics illustrating the lattice-free approach developed by Meineke et al[16]. (A) In silico cells attached by springs. According to Hooke’s law, the restoring force due to each spring is proportional to its net extension and acts in the opposite direction. The total force acting on a cell i is equal to the sum of all forces applied by neighbouring cells j: F_i (t) = μΣ u_i,j (t) ( d_i,j (t) - L_i,j (t) ), with μ the spring constant, L_i,j the equilibrium length of the spring between i and j, d_i,j the distance between i and j, and u_i,j the unit vector from i to j. The new position of cell i after a small time interval Δt is given by: r_i (t + Δt) = r_i (t) + F_i (t) Δt/η, with η a damping constant. (B) In silico 2D tissue, showing cell centres (black nodes) and the associated Voronoi tesselation (solid lines). The related Delauney triangulation (dashed lines) is obtained by connecting all neighbouring nodes and defines the network of springs.

Smallwood and colleagues[39] have used an agent-based framework to develop an alternative, lattice-free model that describes the expansion of a population of epithelial cells cultured in vitro on a monolayer. Individual cells are represented as three-dimensional, rigid spheres that can proliferate, migrate, differentiate or die. The model also accounts for cell-cell and cell-substrate interactions. At every time step, a set of simple predefined rules determines the behaviour of each cell, based on its position, internal parameters and the composition of the local, extracellular environment. Under the assumption that calcium enhances cell-cell adhesion, the resulting computational model showed substantial differences in monolayer growth in the presence of low and physiological calcium concentrations, which nicely mimicked the results from equivalent in vitro studies with normal human urothelial cells. Agent-based models could, in principle, be modified to describe the intestinal epithelium and, in particular, be used to study the impact of changes in cadherin-mediated (Ca⁺⁺-dependent) cell-cell adhesion on the dynamics of the crypt.

MODELS FOR CRYPT BUDDING AND FISSION

The genetic changes that underlie colorectal carcinogenesis usually manifest themselves at a tissue level through the formation of an adenoma, which may subsequently develop into a malignant tumour (Figure 5). This progression is known as the adenoma-carcinoma sequence[7]. Two competing theories have been postulated to explain how such adenomas form. According to the top-down hypothesis, dysplastic cells in the superficial portions of the crypt spread laterally, invading healthy crypts from the ‘top-down’. According to Shih et al[4], dysplastic cells are found to reside on the luminal surface either because the initial APC mutant has migrated to the surface, or because genetic precursors of the dysplastic cells reside in the intercryptal zones; the possibility of a second genetic hit in the migrating population, as suggested by Komarova and Wang[18] and discussed above, could also explain the appearance of a dysplastic luminal epithelial layer. The alternative mechanism for adenoma formation is known as the bottom-up hypothesis[3]. Here adenoma formation always starts with a monocryptal lesion, which then expands by crypt fission and/or by colonisation of neighbouring crypts.

Figure 5 The stages of tumour development. Colorectal cancers develop through a variety of stages. Disruption of Wnt signalling (usually through APC mutation) will cause adenoma formation. This is manifest as dysplasia of the epithelial cells (characterised by enlargement of the epithelial nuclei, loss of mucin and deeply staining nuclei). Adenomas can grow by crypt branching and colonisation of other crypts but, importantly, remain bounded by the basement membrane. The branching of the crypts causes adenomas to develop a very complex architecture and, as more mutations accrue, the architecture becomes more complex and the cells become more severely dysplastic. Eventually the tumour cells will break through the basement membrane and invade into the stroma and this represents the transition from and adenoma into a carcinoma. Once within the stroma, the tumour cells can grow, infiltrate more deeply into the bowel wall and also can invade into lymphatic and vascular structures to allow metastasis to distant sites, such as lymph nodes. It is noteworthy that metastatic deposits also contain stroma showing the importance of the relationship between the tumour epithelium and the tumour stroma.

Adenomatous crypts are typically elongated and deformed, exhibiting multiple branching events, although still lined by a single layer of epithelial cells[40-42] (Figure 5). In order to explain how crypt fission may occur, Drasdo and Loeffler[43] developed an agent-based model in which a vertical cross-section through a crypt was treated as a U-shaped chain of growing, deformable cells that interact elastically with each other. The authors show that decreasing the cell cycle time and reducing the Young modulus of the cells can cause the crypt to buckle and undergo crypt fission. These results are consistent with experimental observations, which suggest that an increase in cell proliferation accompanied by a reduction in cell-cell adhesion (both perhaps associated with a mutation in APC) may lead to crypt fission.

More recently Edwards and Chapman[44] have developed a biomechanical model of crypt budding and fission in which the epithelial cells are treated as a continuous tissue layer that is attached to a rigid basal lamina (Figure 6). The model input parameters quantify cell proliferation, adhesion, movement and apoptosis, which are dependent on lower-scale effects. Using a combination of numerical and analytical techniques, Edwards and Chapman were able to identify parameter regimes in which a flat (rolled-out) crypt would buckle. This instability was indicative of crypt fission. For example, they determined the extent to which the proportion of proliferating cells can be increased before the layer buckles. Also, in qualitative agreement with Drasdo and Loeffler[43], they found that decreasing the cell cycle time could destabilise the layer. As Edwards and Chapman treat the epithelial cells as a continuous sheet rather than as individual entries, it may prove easier to extend their model to simulate more realistic three-dimensional geometries and larger tissue regions, possibly containing multiple crypts.

Figure 6 Schematic of the biomechanical model developed in Edwards and Chapman[44]. Proliferation in the epithelium generates stress within the layer through cellular attachment to the underlying basal lamina. For certain choices of the model parameters, the stress that is generated causes the layer to buckle, inducing a displacement y and initiating crypt budding and fission.

DISCUSSION

In this article we have reviewed a range of mathematical models of CRC that provide insight into different aspects of its development. These include stochastic models designed to investigate where in the crypt the second APC hit associated with CRC initiation is most likely to occur[18], detailed models of Wnt signalling to enhance our understanding of the mechanisms of action of this pathway[36,37] and continuous models to explore the causes of crypt buckling and fission[43,44]. While useful, models that focus on a single time and length scale provide limited insight into crypt dynamics and CRC.

Advances in biotechnology mean that it is now possible to generate vast amounts of different types of experimental data. For example, high-throughput gene and protein expression data can be obtained from microarrays[45-47] and more detailed information from microdissection and tagging specific genes with green fluorescent protein (GFP)[19,48]. Multiscale mathematical models that account for subcellular, cellular and tissue scale phenomena represent a natural framework for integrating and exploiting these different datasets. As part of a UK government funded project in Integrative Biology[49], we are using grid technology and high performance computing to perform multiscale simulations of cardiovascular disease[25,27] and cancer[12,37,44,50]. The protype multiscale model of a colonic crypt that we are developing uses a lattice-free model to describe epithelial cell movement and will incorporate subcellular, biochemical phenomena, such as the Wnt signalling pathway and detailed models of the cell cycle[51]. For example, using our Wnt signalling model[37], it is possible to determine β-catenin levels in the nucleus (and hence involved in gene expression) and the number of molecules forming adhesion complexes, and how these values change with the extracellular Wnt stimulus. By incorporating this information into our multiscale model, we will investigate the way in which Wnt signalling co-ordinates cell movement, proliferation and differentiation along the crypt axis. Additionally, with knowledge of how single and double mutations in the genes coding for APC and β-catenin affect the Wnt signalling pathway, we may use our model to simulate the spread of mutant cells and establish whether their ability to colonise a crypt depends on the spatial location at which the mutation originates. This will provide a comparison of the top-down and bottom-up theories of adenoma formation. By further extending our model to account for the attachment of epithelial cells to the underlying tissue stroma, it will be possible to see how increases in cell proliferation associated with CRC and changes in cell-stroma adhesion contribute to crypt fission. Finally, by upscaling or homogenising our cell-level models, we aim to derive tissue scale models similar in form to that of Edwards and Chapman[44]. In addition to permitting the simulation of large tissue regions (and hence the colonisation of neighbouring crypts by malignant cells), this will enable us to relate tissue level parameters, associated with, for example, cell proliferation and adhesion to experimentally measurable quantities, such as nuclear β-catenin levels and extracellular concentrations of Wnt factors.

Mathematical modeling and computational methods are not only powerful tools for enhancing our understanding of how tumours originate and grow; they may also be used to improve cancer diagnosis and treatment. Hence, modern imaging techniques make it possible to detect and stage the primary tumour and to identify the presence of early metastases, thereby facilitating patient management decisions[52]. Furthermore, a number of mathematical models have been developed to investigate tumour invasion[53,54] and metastasis[55,56]. Such models could, in principle, be adapted to evaluate the aggressiveness of colorectal neoplasms. Finally, we anticipate that in silico studies will be exploited to predict a patient’s prognosis following cancer therapy and/or surgery[23,52,57].

In summary, we believe that mathematical modeling has an important role to play in advancing our under-standing of CRC. In addition to being used to test biological hypotheses (concerning, for example, adenoma formation), well-validated models have the potential for generating new theories (regarding, for example, the ways in which phenomena at different spatial and temporal scales interact) that will themselves stimulate further experimental work. In the longer term, and perhaps most importantly, realistic simulation tools of crypt turnover should assist in the development of new drugs and the optimisation of existing therapies for treating CRC.

ACKNOWLEDGMENTS

The authors gratefully acknowledge financial support from the EPSRC through the Integrative Biology project (IvL, CME) and through an Advanced Research Fellowship (HMB). The authors also thank Professors Philip Maini and John King for helpful discussions.