Homeostasis, injury, and recovery dynamics at multiple scales in a self-organizing mouse intestinal crypt

  1. Louis Gall  Is a corresponding author
  2. Carrie Duckworth
  3. Ferran Jardi
  4. Lieve Lammens
  5. Aimee Parker
  6. Ambra Bianco
  7. Holly Kimko
  8. David Mark Pritchard
  9. Carmen Pin  Is a corresponding author
  1. Clinical Pharmacology and Quantitative Pharmacology, Clinical Pharmacology and Safety Sciences, R&D, AstraZeneca, United Kingdom
  2. Institute of Systems, Molecular and Integrative Biology, University of Liverpool, United Kingdom
  3. Preclinical Sciences and Translational Safety, Janssen, Belgium
  4. Gut Microbes and Health Programme, Quadram Institute, United Kingdom
  5. Clinical Pharmacology and Safety Sciences, AstraZeneca, United Kingdom

Abstract

The maintenance of the functional integrity of the intestinal epithelium requires a tight coordination between cell production, migration, and shedding along the crypt–villus axis. Dysregulation of these processes may result in loss of the intestinal barrier and disease. With the aim of generating a more complete and integrated understanding of how the epithelium maintains homeostasis and recovers after injury, we have built a multi-scale agent-based model (ABM) of the mouse intestinal epithelium. We demonstrate that stable, self-organizing behaviour in the crypt emerges from the dynamic interaction of multiple signalling pathways, such as Wnt, Notch, BMP, ZNRF3/RNF43, and YAP-Hippo pathways, which regulate proliferation and differentiation, respond to environmental mechanical cues, form feedback mechanisms, and modulate the dynamics of the cell cycle protein network. The model recapitulates the crypt phenotype reported after persistent stem cell ablation and after the inhibition of the CDK1 cycle protein. Moreover, we simulated 5-fluorouracil (5-FU)-induced toxicity at multiple scales starting from DNA and RNA damage, which disrupts the cell cycle, cell signalling, proliferation, differentiation, and migration and leads to loss of barrier integrity. During recovery, our in silico crypt regenerates its structure in a self-organizing, dynamic fashion driven by dedifferentiation and enhanced by negative feedback loops. Thus, the model enables the simulation of xenobiotic-, in particular chemotherapy-, induced mechanisms of intestinal toxicity and epithelial recovery. Overall, we present a systems model able to simulate the disruption of molecular events and its impact across multiple levels of epithelial organization and demonstrate its application to epithelial research and drug development.

Editor's evaluation

The proposed model makes an important contribution to the field, allowing a better understanding of the formation and response dynamics of the intestinal crypt through the effective evaluation of healthy, disease and treatment conditions. The authors provided convincing evidence of the validity of their model and their conclusions.

https://doi.org/10.7554/eLife.85478.sa0

Introduction

The intestinal tract is lined by a cellular monolayer which is folded to form invaginations, called crypts, and protrusions, called villi, in the small intestine. The stem cell niche is formed by intermingling Paneth and stem cells located at the base of the crypt (Barker et al., 2007). Stem cells divide symmetrically, forming a pool of equipotent cells that replace each other following neutral drift dynamics (Lopez-Garcia et al., 2010). Continuously dividing stem cells at the base of the crypt give rise to secretory and proliferative absorptive progenitors that migrate towards the villus, driven by proliferation-derived forces (Parker et al., 2017). The transit-amplifying region above the stem cell niche fuels the rapid renewal of the epithelium. The equilibrium of this dynamic system is maintained by cell shedding from the villus tip into the gut lumen (Wright and Alison, 1984).

Epithelial cell dynamics is orchestrated by tightly regulated signalling pathways. Two counteracting gradients run along the crypt–villus axis: the Wnt gradient, secreted by mesenchymal and Paneth cells at the bottom of the crypt, and the bone morphogenetic protein (BMP) gradient generated in the villus mesenchyme, with BMP inhibitors secreted by myofibroblasts and smooth muscle cells located around the stem cell niche (Gehart and Clevers, 2019). These two signalling pathways are also the target of stabilizing negative feedback loops comprising the turnover of Wnt receptors (Hao et al., 2012; Koo et al., 2012; Clevers, 2013b; Clevers and Bevins, 2013c) and the modulation of BMP secretion (Büller et al., 2012; van den Brink et al., 2004). Paneth cells and mesenchymal cells surrounding the niche also secrete other proliferation-enhancing molecules such as epidermal growth factor (EGF) and transforming growth factor-α (Gehart and Clevers, 2019). In addition, Notch signalling-mediated lateral inhibition mechanisms are essential for stem cell maintenance and differentiation into absorptive and secretory progenitors (Gehart and Clevers, 2019). There is also an increasing awareness of the importance of the mechanical regulation of cell proliferation through the Hippo signalling pathway interplaying with several of the key signals, such as EGF, WNT, and Notch, although the exact mechanisms are not currently fully understood (Gehart and Clevers, 2019).

The imbalance of this tightly orchestrated system contributes to pathological conditions, including microbial infections, intestinal inflammatory disorders, extra-intestinal autoimmune diseases, and metabolic disorders (Chelakkot et al., 2018). In addition, critically ill patients and patients receiving chemotherapy/radiotherapy often show severely compromised intestinal barrier integrity (Chelakkot et al., 2018). For instance, oncotherapeutics-induced gastrointestinal toxicity is frequently a life-threatening condition that leads to dose reduction, delay, and cessation of treatment and presents a constant challenge for the development of efficient and tolerable cancer treatments (Stein et al., 2010; Saltz et al., 2000; Saltz et al., 2001; McQuade et al., 2016). This intestinal toxicity often results from the interaction of the drug with its intended molecular target such as cell cycle proteins (Zhang et al., 2021) or the disruption of the cycle through DNA damage (Helleday et al., 2008). Multiscale models integrating our knowledge on how the epithelium maintains homeostasis and responds to injury can contribute to understand epithelial biology and quantify the risk of intestinal toxicity during drug development.

Several agent-based models (ABMs) have been proposed to describe the complexity and dynamic nature of the intestinal crypt. Early models were used as in silico platforms to study the dynamics and cellular organization of the crypt. For instance, one of the pioneering ABMs was used to study the distribution and organization of labelling and mitotic indices (Meineke et al., 2001). This model comprises a fixed ring of Paneth cells beneath a row of stem cells, which divide asymmetrically to produce a stem cell and a transit-amplifying cell that terminally differentiates after a fixed number of divisions. Some subsequent models are lattice-free, recapitulate neutral drift of equipotent stem cells, and describe proliferation and cell fate regulated by a fixed Wnt signalling spatial gradient, which is defined by the distance from the crypt base, with proliferating cells progressing through discrete phases of the cell cycle and showing variable duration of the G1 phase (Pitt-Francis et al., 2009). Further model refinements can be seen in the model of Buske et al., 2011, with stochastic cell growth and division time (Buske et al., 2011), Wnt levels defined by the fixed local curvature of the crypt and lateral inhibition driven by Notch signalling. Here, we present a lattice-free ABM that describes the spatiotemporal dynamics of single cells in the small intestinal crypt driven by the interaction of surface-tethered Wnt signals, cell–cell Notch signalling, BMP-diffusive signals, RNF43/ZNRF3-mediated feedback mechanisms, and the cycle protein network responding to the crypt mechanical environment. We show that our computational model enables the simulation of the ablation and recovery of the stem cell niche as well as of how drug-induced molecular perturbations trigger a cascade of disruptive events spanning from the cell cycle to single-cell arrest and/or apoptosis, altered cell migration and turnover, and ultimately loss of epithelial integrity.

Results

Modelling a self-organizing crypt using an ABM

We have modelled the mouse intestinal crypt as a self-organizing system where cell dynamics and cell composition arise from local interactions between single cells and the mesenchyme through signalling pathways with behaviours (proliferation, differentiation, fate decision, migration, etc.) determined largely by endogenous intracellular and intercellular interactions.

The model describes the spatiotemporal dynamics of stem cells and progenitors undergoing division cycles and responding to intercellular signalling to differentiate into Paneth, goblet, and enteroendocrine cells and enterocytes (Figure 1A). All cells interact physically and biochemically in the geometry of the crypt. Stem cells intermingle with Paneth cells at the bottom of the crypt and randomly replace each other. Progenitors and mature cells migrate towards the villus driven by proliferation forces (Figure 1A). To achieve a stable crypt cell composition under constant cell renewal dynamics, we have implemented several signalling mechanisms which include the Wnt, Notch, and BMP pathways essential for morphogenesis and homeostasis of the intestinal crypt (Gehart and Clevers, 2019; Fevr et al., 2007; VanDussen et al., 2012; Pellegrinet et al., 2011; He et al., 2004), the YAP-Hippo signalling pathway responding to mechanical forces and modulating contact inhibition of proliferation (Gjorevski et al., 2016), and a ZNRF3/RNF43-like-mediated feedback mechanism between Paneth and stem cells to regulate the size of the stem cell niche according to experimental reports (Hao et al., 2012; Koo et al., 2012; Farin et al., 2016; Figure 1B).

Schematics of the small intestinal crypt composition and cell fate signalling pathways included in the agent-based model (ABM).

(A) Depiction of the crypt highlighting key signalling features and cell types in each crypt region. (B) Details of signalling pathways including the formation of the Wnt signalling gradient with high levels of Wnt in the stem cell niche generated by Paneth and mesenchymal cells. Intercellular pressure regulates the duration of the division cycle (YAP-Hippo pathway-mediated contact inhibition of proliferation) which impacts on the accumulation of cell surface-tethered Wnt signals. Notch signalling maintains the balance between Paneth and stem cells through lateral inhibition. A ZNRF3/RNF43-mediated feedback mechanism modulates Wnt signalling in the niche restricting the number of stem and Paneth cells. BMP signals generated by mature villus cells form a feedback loop that regulates maturation and proliferation of absorptive progenitors. (C) Cell fate determination. High Wnt signalling and activation of Notch are required to maintain stemness. Low Notch signalling determines differentiation into secretory fates, including Paneth cells in high Wnt signalling regions, or goblet/enteroendocrine progenitors in low Wnt regions. Absorptive progenitors develop from stem cells in low Wnt conditions and divide 3–5 times, before becoming terminally differentiated when Wnt signal levels are decreased and cells find sufficient BMP signals. (D) Average composition of a simulated healthy/homeostatic crypt (over 100 simulated days), showing the relative proportion of cells at each position.

The Wnt pathway is the primary pathway associated with stem cell maintenance and cell proliferation in the crypt (Fevr et al., 2007; van der Flier and Clevers, 2009). Our model implements two sources of Wnt signals described in the crypt: Paneth cells (Sato et al., 2011) and mesenchymal cells surrounding the stem cell niche at the crypt base (Stzepourginski et al., 2017). Wnt signalling is modelled as a short-range field around Wnt-emitting Paneth and mesenchymal cells with Wnt signals tethered to receptive cells as previously reported (Farin et al., 2016; Clevers and Nusse, 2012). Surface-tethered signals are split between daughter cells upon cell division (Gehart and Clevers, 2019; Farin et al., 2016), which results in a gradual depletion of tethered Wnt signals as cells divide and migrate towards the villus away from Wnt sources (Figure 1A and B). Notch signalling is also implemented in the model with Notch ligands expressed by secretory cells binding to Notch receptors on neighbouring cells and preventing them from differentiating into secretory fates, in a process known as lateral inhibition, that leads to a checkerboard/on-off pattern of Paneth and stem cells in the niche (VanDussen et al., 2012). Specifically, in our model, high Wnt and Notch signalling environments are required to maintain stemness, as reported in the literature (Tian et al., 2015) while under low Notch and high Wnt signalling, stem cells differentiate into secretory cells, including Paneth cells. On the other hand, Notch signalling also mediates the process of Paneth cell de-differentiation into stem cells to regenerate the niche as previously reported (Mei et al., 2020; Yu et al., 2018). Stem cells with decreased levels of Wnt signalling, usually located outside the niche, differentiate into absorptive proliferating progenitors or alternatively into secretory progenitors in the absence of Notch signals (Figure 1C).

In our model, mechanical stimuli, captured through the YAP-Hippo signalling pathway (Gjorevski et al., 2016; Halder et al., 2012; Aragona et al., 2013; Low et al., 2014), indirectly interact with the Notch and Wnt signalling pathways. We recapitulate YAP-mediated contact inhibition of proliferation by using cell compression to modulate the duration of the division cycle which increases when cells are densely squeezed, such as in the stem cell niche, and decreases if cell density falls, for instance, in the transit-amplifying compartment or in cases of crypt damage (Figure 1A and B). In agreement with experimental reports (Pin et al., 2015), in our model, Paneth cells are assumed to be stiffer and larger than other epithelial cells, requiring higher forces to be displaced and generating high intercellular pressure in the niche. Due to the increased mechanical pressure, cells in the niche have longer division cycles and can accumulate more Wnt and Notch signals. These premises imply that Paneth cells enhance their own production by generating Wnt signals and inducing prolonged division times, which increases stem and Paneth cell production and could lead to unlimited expansion of the niche recapitulating the phenotype seen in ZNRF3/RNF43 knockout mice (Koo et al., 2012; see Appendix 1, Section 1.11). To generate a niche of stable size, we implemented a negative Wnt-mediated feedback loop that resembles the reported stem cell production of RNF43/ZNRF3 ligands to increase the turnover of Wnt receptors in nearby cells (Hao et al., 2012; Koo et al., 2012; Clevers, 2013b; Clevers and Bevins, 2013c). Similarly, in our model, a number of stem cells in excess of the homeostatic value reduces cell tethering of Wnt ligands and hence inhibits Paneth and stem cell generation (Figure 1A and B).

The Wnt gradient in the crypt is opposed by a gradient of bone morphogenic protein (BMP) that inhibits cell proliferation and promotes differentiation (Qi et al., 2017). We assume that enterocytes secrete diffusing signals, resembling Indian Hedgehog signals (Büller et al., 2012), that induce mesenchymal cells to generate a BMP signalling gradient effective to prevent proliferative cells from reaching the villus (Figure 1A and B). Based on experimental evidence, we also assume that BMP activity is counteracted by BMP antagonist-secreting mesenchymal cells surrounding the stem cell niche (McCarthy et al., 2020). Proliferative absorptive progenitors migrating towards the villus lose Wnt during every division and eventually meet values of BMP that overcome the proliferation-inducing effect of Wnt signalling (He et al., 2004). We found that a homeostatic crypt cell composition is achieved when BMP and Wnt differentiation thresholds result in progenitors dividing approximately four times before differentiating into enterocytes (Figure 1C). In our model, the BMP signalling gradient responds dynamically to the number of enterocytes, giving rise to a negative feedback loop between enterocytes on the villus and their proliferative progenitors in the crypt that recapitulates the enhanced crypt proliferation observed after epithelial damage (Büller et al., 2012; Pont and Yan, 2018; Sprangers et al., 2021). For instance, a decreased number of enterocytes results in reduced production of BMP, which enables progenitor cells to divide and migrate further up the crypt before meeting BMP levels higher than the differentiation threshold.

Altogether our model describes single cells that generate and respond to signals and mechanical pressures in the crypt–villus geometry to give rise to a self-organizing crypt which has stable spatial cell composition over time (Figure 1D) and reproduces reported experimental data (Buske et al., 2011). An extended description of these modelling features is provided in Appendix 1.

The cell cycle protein network governs proliferation in each single cell of the ABM and responds to mechanical cues

We have used the model of Csikász-Nagy et al., 2006, which is based on the seminal work of Novak and Tyson, 1993; Novak et al., 2001; Novák and Tyson, 2004 and available in BioModels (Le Novère and Csikasz-Nagy, 2006), to recreate the dynamics of the main proteins governing the mammalian cell cycle in each single proliferative cell of the ABM. In this model, a dividing cell begins in G1, with low levels of cyclins A, B, and E and a high level of Wee1, and progresses to S-phase when cyclin E increases. S-phase ends and G2 begins when Wee1 falls. The decrease in cyclin A expression defines the start of M-phase, while falling cyclin B implies the end of M-phase, when the cell divides into two daughter cells with half the final mass value and re-enters the cell cycle (Figure 2A–D).

Multiscale modelling of cell division in single cells of the agent-based model (ABM).

(A, B) Modelled dynamics of the main cell cycle proteins across the phases of division in each single cell over a 24 hr period, according to the cell cycle regulatory protein network model of Csikasz-Nagy (Csikász-Nagy et al., 2006). The protein interaction diagram can be found in the original report of Csikasz-Nagy (Csikász-Nagy et al., 2006). Stem cells in the crowded niche (A) exhibit longer cycles, up to 21.5 hr on average, with elevated levels of p27 regulating the duration of G1 and the starting of S-phase. Cells in the transit-amplifying compartment (B) have shorter cycles, up to 10 hr on average, due to low levels/lack of p27 expression which leads to G1 shortening and early start of the S-phase. A.U., arbitrary units. (C) Observed (dashed line) and simulated (solid line) proportions of BrdU-positive cells at each crypt position at 2 hr (blue), 24 hr (purple), and 80 hr (red) after a single pulse of BrdU. (D) Observed and simulated Ki-67-positive cells at each crypt position assuming that Ki-67 is detected in cycling cells at all phases except G1 and in any recently differentiated and arrested cells. Shadows depict the 95% confidence interval of our simulated staining results assuming that the proportion of staining cells has a beta distribution and estimating its error from experimental data.

To implement YAP-Hippo-mediated contact inhibition of proliferation, we have modified the dynamics of the proteins of the Csikasz-Nagy model to respond to mechanical cues encountered by cells migrating along the crypt. Crowded, constrained environments result in longer cycles, such as in stem cells in the niche, while decreased intercellular forces lead to shortened cycles as cells migrate towards the villus in agreement with experimental reports (Wright and Alison, 1984; Marshman et al., 2002; Potten et al., 1997). The shorter cycle duration in absorptive progenitors has been mainly associated with shortening/omission of G1, while the duration of S-phase is less variable (Wright and Alison, 1984). Using the model of Csikász-Nagy et al., 2006, we modulated the duration of G1 through the production rate of the p27 protein. The p27 protein has been reported to regulate the duration of G1 by preventing the activation of cyclin E-Cdk2 which induces DNA replication and the beginning of S-phase (Morgan and Morgan, 2007). We, hence, hypothesized that rapid cycling absorptive progenitors located in regions of low mechanical pressure outside the stem cell niche have low levels of p27, which bring forward the start of S-phase to shorten G1 (Figure 2D). In support of this hypothesis, it has been demonstrated that p27 inhibition has no effect on the proliferation of absorptive progenitors (Zheng et al., 2008; see Appendix 1 for a full description). These new features of the cell cycle model are updated dynamically and continuously to respond to changes in mechanical pressure experienced by each cell as it migrates along the crypt.

To demonstrate the performance of the model to reproduce the spatiotemporal cell dynamics and composition of a homeostatic crypt, we simulated previous published mouse experiments (Parker et al., 2017; Parker et al., 2019) comprising 5-bromo-29-deoxyuridine (BrdU) tracking (Figure 2E) and Ki-67 staining (Figure 2F). BrdU is a thymidine analogue often used to track proliferative cells and their descendants along the crypt–villus axis (Nowakowski et al., 1989; Gratzner, 1982). BrdU is incorporated into the newly synthesized DNA of dividing cells during S-phase and transmitted to daughter cells, regardless of whether they proliferate. If the exogenous administration of this molecule is discontinued, the cell label content is diluted by each cell division and is no longer detected after 4–5 generations (Wilson et al., 2008). To simulate the BrdU chase experiment after a single BrdU pulse, we assumed that any cell in S-phase incorporated BrdU permanently into its DNA for the first 120 min after injection of BrdU and BrdU cell content was diluted upon cell division such that after five cell divisions, BrdU was not detectable. See Appendix 1 for a complete description. The BrdU chase simulation showed that the observed initial distribution of cells in S-phase as well as division, differentiation, and migration of BrdU-positive cells over time was replicated by our model (Figure 2E).

Ki-67 is a protein produced by actively proliferating cells during the S-, G2-, and M-phases of the division cycle (Sobecki et al., 2017). Due to the time required for this protein to be catabolized (Miller et al., 2018), Ki-67 is also detected in quiescent or non-proliferative cells after exiting the cycle (Miller et al., 2018) and during G1 in continuously cycling cells (Sobecki et al., 2017). Our simulations assumed that Ki-67 is detected in continuously cycling cells, cells re-entering the cycle after arrest except during G1, as well as in differentiated cells that were cycling within the past 6 hr and recently drug-arrested cells. See Appendix 1 for a complete description. Similarly, we observed that the ABM-simulated spatial distribution along the crypt of Ki-67-positive cells recapitulated observations in mouse ileum (Figure 2F).

In summary, proliferative cells in the ABM respond to mechanical cues by adjusting the cell cycle protein network to dynamically change the duration of the cycle while migrating along the crypt. With this feature, the model replicates spatiotemporal patterns of cell proliferation, differentiation, and migration observed in mouse experiments.

Cell plasticity/de-differentiation enables crypt regeneration following damage of the stem cell niche

Marker-based lineage-tracing studies have demonstrated numerous potential sources available for intestinal stem cell regeneration (Hageman et al., 2020). In line with these studies, our model assumes that cell fate decisions are reversible and both secretory and absorptive cells are able to revert into stem cells when regaining sufficient Wnt and Notch signals.

To investigate the potential of the ABM to describe and explore cell plasticity dynamics, we simulated the repeated ablation of intestinal stem cells resembling a previously published study (Tan et al., 2021). Following the experimental setup in that study, we simulated the diphtheria toxin receptor-mediated conditional targeted ablation of stem cells for four consecutive days considering that ablation was completed after the first 24 hr (Saito et al., 2001) and persistently inducing stem cell death during the remaining days of treatment (Figure 3A–C). Our simulations showed that 6 hr after the last induction, stem cells were not detected, Paneth cells decreased by 75–100% (Figure 3B), and the villus length was reduced by about 10–20% (Figure 3C) which was similar to the reported experimental findings (Tan et al., 2021). Simulated proliferative absorptive progenitors were indirectly affected by stem cell ablation and their decrease was followed by a reduction in mature enterocytes. The progenitors recovered after treatment interruption to later reach values above baseline when responding to the negative feedback signalling from mature enterocytes (Figure 3A). In our simulations, enhanced crypt proliferation was not accompanied by simultaneous villus recovery, which started later. Tan et al., 2021 reported similar results with increased crypt proliferation replenishing first the crypt and not contributing immediately to villus recovery. See Video 1 to visualize the response of the crypt.

Simulated cell dynamics in the epithelium subjected to continuous ablation of stem cells for four consecutive days (grey block) resembling a previously published experiment (A–C) (Tan et al., 2021).

Analysis time denotes 6 hr after ablation interruption for comparison with reported results (Tan et al., 2021). All cell lineages are recorded during treatment and few days after recovery of the simulated crypt for comparison with homeostasis. A simulated 3D image of a crypt in homeostasis can be found in Figure 4A. (A) shows the total number of cells, absorptive progenitors, and enterocytes in the crypt; (B) shows the number of Paneth, stem cells, and uncommitted progenitors, mostly found in the niche; and (C) shows villus cells. (D) Relative frequency of crypt cells moving towards the villus (darker colour) and towards the crypt base, that is, retrograde motion (lighter colour), in homeostasis (blue) and during stem cell ablation (red) at each cell position, showing increased retrograde cellular motion in the niche following stem cell ablation. (E) Leftmost: trajectories (cell position on crypt–villus longitudinal axis vs time) of the progeny of one stem cell, with both daughters leaving the niche and giving rise to a cascade of absorptive and secretory cells that eventually leave the crypt. Rightmost: trajectories of the progeny of an absorptive progenitor dedifferentiating into a stem cell during recovery after stem cell ablation.

Video 1
Simulated cell dynamics in the epithelium subjected to continuous ablation of stem cells for four consecutive days resembling a previously published experiment (Tan et al., 2021).

Plots depict changes in the number of cells in the crypt and villus during the simulation. Colour code of cell types is included below plots.

We next studied the type of cells that were dedifferentiating during the simulated repeated ablation of stem cells and found that in agreement with experimental reports, Paneth cells (Yu et al., 2018), absorptive progenitors (Tetteh et al., 2016), and quiescent stem cells located just above the stem cell niche at the fourth cell position from the crypt base (Tian et al., 2011) dedifferentiated into stem cells. Specifically, from all dedifferentiated cells, about 60% were Paneth cells, 30% absorptive progenitors, and 10% secretory progenitors, which are considered quiescent stem cells as previously suggested (Buczacki et al., 2013). Furthermore, we used our model to explore the retrograde motion, reported using intravital microscopy (Azkanaz et al., 2022), of cells returning to the niche to de-differentiate into stem cells. For cells outside the niche, movement is retrograde when its velocity is negative in the z direction, that is, they move towards the niche across the longitudinal crypt–villus axis. For cells in the hemispherical niche, we consider a cell to move forwards, towards the villus, or backwards, towards the crypt base, if the rate of change of its polar angle is positive or negative, respectively. This implies that cells can be recorded to move backwards despite being located at the crypt base. We observed that the frequency of retrograde, or backward, movements is relatively high at low positions in a crypt in homeostasis (Figure 3D) and increases further after stem cell ablation, reflecting increased retrograde cellular motion as cells repopulate the niche. While in homeostasis the progeny of a stem cell generally differentiates into a cascade of absorptive and secretory progenitors that migrate towards the villus and eventually leave the crypt (Figure 3E). Following the interruption of stem cell ablation, during recovery absorptive progenitors return to the niche and dedifferentiate to regenerate multiple stem and Paneth cells as well as progenitors (Figure 3E).

Taken together, our model recapitulates cellular reprogramming of both multipotent precursors and committed progeny in the crypt and replicates the reported crypt injury dynamics following persistent ablation of stem cells (Tan et al., 2021).

Disturbance of cell cycle proteins spans across scales to impact on crypt and villus organization

The model of Csikász-Nagy et al., 2006 enables the simulation of the disruption of the main proteins governing the cell cycle in each single proliferative cell of the ABM. CDKs play important roles in the control of cell division (Malumbres, 2014), and the development of CDK inhibitors for cancer treatment is an active field of research (Zhang et al., 2021).

To explore the effect of the disruption of the cell cycle on epithelial integrity, we simulated the inhibition of CDK1 for 6 hr, every 12 hr for four consecutive days, resembling epithelial toxicity of a theoretical drug. CDK1 is reported to be the only CDK essential for the cell cycle in mammals (Santamaría et al., 2007). CDK1 triggers the initiation of cytokinesis by inducing the nuclear localization of mitotic cyclins A and B (Pesin and Orr-Weaver, 2008), and its inhibition has been proposed as a cancer therapy with potentially higher efficacy than the inactivation of other CDKs (Diril et al., 2012). To mimic CDK1 inhibition, we added a term to the CycA/CDK1,2 and CycB/CDK1 differential equations of the Csikasz-Nagy model (Csikász-Nagy et al., 2006) that strongly reduces the production of both CycA/CDK1,2 and CycB/CDK1 during the CDK1 inhibition period (Figure 4A–E; Appendix 1).

Simulation of CDK1 inhibition for 6 hr, every 12 hr for four consecutive days in the agent-based model (ABM) and impact on the cell cycle and crypt and villus organization.

All cell lineages are recorded during treatment and a few days after recovery of the simulated crypt for comparison with homeostasis.

(A) A simulated 3D image of a crypt in homeostasis (left) and a crypt subjected to CDK1 inhibition (right). Following CDK1 inhibition, the simulated crypt exhibits apoptotic cells and oversized cells unable to correctly complete the cell cycle and eventually undergo cell cycle arrest. Colour code provided here for apoptotic and arrested cells and in Figure 1A for the rest of cells. (B) Flowchart showing the regular progression through the cell cycle (green path) disturbed by CDK1 inactivation. A disorderly restart of the cycle, leading to enlarged cells, is observed when CDK1 inhibition prevents cells from entering (yellow path) or completing M-phase (orange path) by early reduction of cyclin B, with premature restart of G1 (orange path). Cells in M-phase subjected to greater reduction of cyclins A and B that completely disrupts the protein network undergo mitotic death (red path). (C) Cell cycle protein dynamics in homeostasis. (D, E) Altered cell cycle protein profile by CDK1 inhibition, resulting in premature restart of G1 and arrest of enlarged cell (D) and in disruption of the protein network and cell death (E). Protein concentrations given in arbitrary units (A.U.). Cell dynamics in simulated crypts (F) and villi (G) during CDK1 inhibition period and recovery. The dynamics of all cell lineages are reported in Appendix 1—figure 1. Discontinuous bars denote the beginning of CDK1 inhibition period.

It has been experimentally demonstrated that the selective inhibition of CDK1 activity in cells programmed to endoreduplicate (i.e. cells that can duplicate their genome in the absence of intervening mitosis) leads to the formation of stable nonproliferating giant cells, whereas the same treatment triggers apoptosis in cells that are not developmentally programmed to endoreduplicate (Ullah et al., 2008). Although endoreduplication is not expected in crypt cells, enlarged polynucleated cells have been reported to remain in the epithelium without dying in a recent light-sheet organoid imaging study tracking the progeny of a cell after cytokinesis failure induced by the inhibition of LATS1 (de Medeiros et al., 2022), which is phosphorylated by CDK1 during mitosis (Furth and Aylon, 2017). Thus, we chose to replicate this phenotype to show the capacity of our model to predict possible complex responses in the intestine. Following CDK1 inhibition, we detected oversized cells in the ABM (Figure 4A). The inhibition of the activation of cyclins A and B altered the modelled protein profiles, disturbing progression through G2 and M-phase and preventing the cell mass from dividing before reinitiating a new cycle (Figure 4B). Thus, a cell could either be (i) unaffected if it was at the early stages of the cycle (Figure 4C); or (ii) restart the cell cycle if CDK1 was inhibited while the cell was at the end of G2 and unable to enter M-phase or in M-phase and unable to complete cytokinesis. In this case, the inhibition of cyclins A and B led to an early increase in cyclin E and the premature restart of G1 with the generation of oversized cells, which are ultimately arrested (Figure 4D); or (iii) cells in M-phase can undergo mitotic death if the reduction of cyclins A and B severely disrupts the protein network (Figure 4E). Hence, the failure to culminate M-phase resulted in cell death or generation of oversized, nonproliferating cells, which led to a reduction of the crypt overall cell number (Figure 4F) and the turnover of villus cells (Figure 4G). Appendix 1—figure 1 shows the response of all cell lineages to CDK1 inhibition, and Video 2 shows the 3-D visualization of the crypt during this treatment.

Video 2
Simulated cell dynamics in the epithelium subjected to CDK1 inhibition for 4 d.

Plots depict changes in the number of cells in the crypt and villus during the simulation. Colour code of cell types is included below plots.

Altogether our ABM enables the simulation of how disruptions of the cell cycle protein network span across scales to generate complex phenotypes, such as giant cells, and impact on the integrity of the crypt and villus structure.

A practical application of the ABM to describe 5-fluorouracil (5-FU)-induced epithelial injury at multiple scales

5-FU is a well-studied and commonly administered cancer drug (Longley et al., 2003) with reported high incidence of gastrointestinal adverse effects in treated patients (Stein et al., 2010). 5-FU is a pyrimidine antimetabolite cytotoxin which has multiple mechanisms of action upon conversion to several nucleotides that induce DNA and RNA damage (Longley et al., 2003). Antimetabolites resemble nucleotides and nucleotide precursors that inhibit nucleotide metabolism pathways, and hence DNA synthesis, as well as impair the replication fork progression after being incorporated into the DNA (Helleday et al., 2008).

To explore the performance of our ABM to predict epithelial injury, we used results from experiments in mice dosed with 50 and 20 mg/kg of 5-FU every 12 hr for 4 d to achieve drug exposures similar to those observed in patients (Jardi et al., 2023). 5-FU pharmacokinetics is metabolized into three active metabolites FUTP, FdUMP, and FdUTP (Longley et al., 2003). Based on previous reports, we assumed that FUTP is incorporated into RNA of proliferative cells, leading to global changes in cell cycle proteins (Pritchard et al., 1997), while FdUTP is incorporated into DNA (Longley et al., 2003) during S-phase, resulting in the accumulation of damaged DNA. In our model, DNA and/or RNA damage can be repaired or lead to cell arrest or apoptosis (Figure 5A). We did not implement the inhibition of thymidylate synthase (TS) by FdUMP because the impact of this mechanism on intestinal toxicity is not completely understood (Pritchard et al., 1997). A previously published 5-FU PK model (Gall et al., 2023) was integrated into the ABM to describe the dynamic profile of the concentration of 5-FU and its metabolites in plasma and GI epithelium after dosing (Figure 5B).

Modelling 5-fluorouracil (5-FU) (50 mg/kg twice a day for 4 d) induced injury at several scales in mouse small intestinal epithelium.

(A) Diagram showing the implemented mechanism in the agent-based model (ABM) to describe DNA and RNA damage and cell cycle disruption driven by 5-FU metabolites. Cells trigger the apoptotic pathway if relatively high levels of RNA and/or DNA damage are detected at the cycle checkpoints. Lower levels of DNA damage induced P21 activation, which, together with RNA damage, slow down and could eventually arrest the cycle. (B) Predicted concentration (ng/ml) of 5-FU, FUTP, and FdUTP in plasma in mouse (pharmacokinetics model of 5-FU described in Gall et al., 2023). (C, D) Cell cycle protein dynamics and fate decision when 5-FU challenge starts (C) prior to or at the beginning of S-phase, leading to DNA damage and cell death at the G2-M-phase checkpoint, and (D) at the end of S-phase, resulting in not enough DNA damage, the cell finishes the cycle. (E) Predicted (solid line) and observed (dashed line) proportions of Ki-67-positive cells along the crypt axis at 6 hr, 1 d, 4 d, and 6 d during the 5-FU treatment period. Shadows depict the 95% confidence interval of our simulated staining results assuming that the proportion of staining cells has a beta distribution and estimating its error from experimental data. (F, G) Predicted (lines) and observed (symbols) number of cells in the crypt (F) and villus (G). Vertical bars represent dosing times. Symbols represent cell counts from individual mice.

Figure 5C shows the cell cycle protein dynamics and fate decision when 5-FU challenge took place at the beginning of S-phase and led to the accumulation of relatively high levels of DNA damage which triggered cell death at the G2-M-phase checkpoint. When the challenged cell was at the end of S-phase, the accumulated levels of DNA damage were not high enough to be detected at the G2-M-phase checkpoint and the cell finished the cycle and restarted a new cycle at a slower rate due to concurrent RNA damage and relatively low level of DNA damage (Figure 5D).

Figure 5E shows that predicted and observed Ki-67-positive cells declined gradually over time at all positions in the crypt during the 5-FU high-dose treatment. However, the numbers recovered, reaching values above baseline, 2 d after the interruption of 5-FU administration. The increased rebound of the proliferative crypt compartment after treatment was captured in our ABM by the implemented BMP-mediated feedback mechanism from mature enterocytes to proliferative cells (see Appendix 1, Section 1.7.4). For this treatment, both simulated and observed total number of cell,s in the crypt followed the same pattern as the proliferative compartment (Figure 5F), while the decline in villus cells started later and took longer to achieve full recovery (Figure 5G). Appendix 1—figure 2A and B shows the response of all cell lineages during this treatment, and Video 3 shows the 3-D visualization of the simulated crypt and changes in signalling pathways and cell composition during the high-dose 5-FU challenge. The low dose of 5-FU had a minor impact on crypt proliferation and villus integrity, which was also recapitulated by the model (Appendix 1—figure 2C–E).

Video 3
Simulated cell and signalling molecular dynamics in the epithelium following the administration of 50 mg/kg of 5-fluorouracil (5-FU) twice a day for 4 d in mouse.

Plots depict changes in signal abundance across the crypt longitudinal axis (z), in the number of cells in the crypt and villus, and concentration of 5-FU and metabolites during the simulation. Signals expressed in arbitrary units (A.U.). Colour code of cell types is included below plots.

Overall, the ABM recapitulates DNA and RNA damage, resulting in cell cycle disruption associated with 5-FU administration and describes the propagation of the injury across scales to disturb epithelial integrity. The loss of epithelial barrier integrity is widely accepted to be the triggering event of chemotherapy-induced diarrhoea (McQuade et al., 2016) which is reported in mice at the doses used in this study (Jardi et al., 2023) as well as observed in patients undergoing equivalent treatments (Morawska et al., 2018).

Discussion

We have built a multi-scale ABM of the small intestinal crypt with self-organizing, stable behaviour that emerges from the dynamic interaction of the Wnt, Notch, BMP, and ZNRF3/RNF43 pathways orchestrating cellular fate and feedback regulatory loops and includes contact inhibition of proliferation, RNA and DNA metabolism, and the cell cycle protein interaction network regulating progression across division stages.

In our model, the stability of the niche is achieved by a negative feedback mechanism from stem cells to Wnt respondent cells that resembles the reported turnover of Wnt receptors by ZNRF3/RNF43 ligands secreted by stem cells (Hao et al., 2012; Koo et al., 2012; Clevers, 2013b; Clevers and Bevins, 2013c). Wnt signals generated from mesenchymal cells and Paneth cells at the bottom of the crypt are tethered to receptive cells and divided between daughter cells upon division, which forms a decreasing Wnt gradient towards the villi that stimulates cell proliferation and ensures stemness maintenance (Farin et al., 2016; Sato et al., 2011). The model also implements the BMP signalling counter-gradient along the crypt–villus axis by resembling the production of diffusive BMP signals by mesenchymal telocytes abundant at the villus base as well as the activity of BMP antagonist molecules secreted by trophocytes located just below crypts (McCarthy et al., 2020). This BMP signalling gradient forms an additional negative feedback mechanism that regulates the size of the crypt proliferative compartment and recapitulates the modulation of BMP secretion by mesenchymal cells via villus cells-derived hedgehog signalling (Büller et al., 2012; van den Brink et al., 2004).

Another novel feature of our model is the inclusion of the dynamics of the protein network governing the phases of cell division (Csikász-Nagy et al., 2006). Moreover, in our model, the cell cycle protein network responds to environmental mechanical cues by adapting the duration of the cycle phases. Cells in crowded environments subjected to higher mechanical pressure, such as stem cells in the niche, exhibit longer cell cycles (Wright and Alison, 1984; Marshman et al., 2002; Potten et al., 1997) while progenitors in the transit-amplifying compartment adapt their cell cycle protein dynamics to mainly shorten G1-phase (Wright and Alison, 1984; Carroll et al., 2018) and proliferate more rapidly. This model feature recapitulates the widely reported YAP-mediated mechanism of contact inhibition of proliferation under physical compression (Halder et al., 2012; Aragona et al., 2013; Low et al., 2014). Interestingly, it has been reported that stiff matrices initially enhance YAP activity and proliferation of in vitro cultured intestinal stem cells by promoting cellular tension (Gjorevski et al., 2016); however, that study also proposes that the resulting colony growth within a stiff confining environment may give rise to compression YAP inactivation retarding growth and morphogenesis (Gjorevski et al., 2016).

Furthermore, our model considers that the mechanical regulation of the cell cycle interacts with signalling pathways to maintain epithelial homeostasis, but also to trigger cell dedifferentiation if required. Cells with longer cycles accumulate more Wnt and Notch signals, leading to the maintenance of the highly dynamic niche by replacement of Paneth and stem cells. Cells located outside the niche exhibit shorter cycles and cannot effectively accumulate enough Wnt signals to dedifferentiate into stem cells in homeostatic conditions. However, in case of niche perturbation, progenitor cells reaching the niche as well as existing Paneth cells in the niche are able to dedifferentiate into stem cells after regaining enough Wnt signals, which replicates the injury recovery mechanisms observed in the crypt (Hageman et al., 2020; Tetteh et al., 2016). Our model also concurs with experimental results suggesting that Lgr5+ stem cells are essential for intestinal homeostasis and that their persistent ablation compromises epithelial integrity (Tan et al., 2021).

Altogether, our model implements qualitative and quantitative behaviours to better simulate the functional heterogeneity of the intestinal epithelium at multiple scales. One of the important applications of our modelling approach lies in the development of safer oncotherapeutics. The model enables the prediction of intestinal injury associated with efficacious dosing schedules in order to minimize toxicity while maintaining the efficacy of investigational drugs. We demonstrated the application of our model to predict potential intestinal toxicity phenotypes induced by CDK1 inhibition as well as describe the disruption of the epithelium at multiple scales triggered by RNA and DNA damage, leading to the loss of integrity of the intestinal barrier and diarrhoea following 5-FU treatment. The drug-induced perturbation of other cell cycle proteins or signalling pathways, already integrated into the model, is straightforward to simulate with the current version of the model while the resolution of molecular networks can be increased, or new pathways incorporated into the ABM, to describe additional drug mechanisms of action.

While most of the crypt biology understanding integrated in our model derives from mouse epithelial studies, human-derived intestinal organoids and microphysiological systems, now routinely used in research, can provide highly precise information at the single-cell level to inform ABM development. In return, ABMs can help test hypotheses behind organoid responses in health and disease conditions. Our work highlights the importance of novel modelling strategies that are able to integrate the dynamics of processes regulating the functionality of the intestinal epithelium at multiple scales in homeostasis and following perturbations to provide unprecedented insights into the biology of the epithelium with practical application to the development of safer novel drug candidates.

Materials and methods

Mouse experiments

Request a detailed protocol

We used BrdU tracking and Ki-67 immunostaining data from previously published experiments in healthy mice (Parker et al., 2017; Parker et al., 2019) and following 5-FU treatment (Jardi et al., 2023). The samples from this later study (Jardi et al., 2023) were analysed again to count Ki-67-positive cells at each position along the longitudinal crypt axis for 30–50 individual hemi crypt units per tissue section per mouse as previously described (Williams et al., 2016).

ABM development

Request a detailed protocol

A comprehensive description of the model can be found in Appendix 1 and Appendix 1—table 1. The model has been made available through BioModels (MODEL2212120002) (Malik-Sheriff et al., 2020)

Appendix 1

Technical description of the intestinal epithelial ABM

The model primarily focuses on describing the spatiotemporal dynamics of single epithelial cells, interacting physically and biochemically in the mouse intestinal crypt, undergoing division cycles or differentiating into mature epithelial cells. Single cells both generate and respond to signals and mechanical pressure in the crypt–villus geometry to generate a self-organizing tissue.

Below we describe the assumptions and hypotheses that underpin the model, regarding (1) geometry; (2) cell cycle proteins and cellular growth; (3) drug perturbation of the cell cycle proteins: Cdk1 inhibition; (4) DNA and RNA synthesis; (5) drug perturbations of RNA and DNA synthesis: 5-FU-induced RNA and DNA damage; (6) mechanical cell interactions and contact inhibition; (7) biochemical signalling; (8) cell fate: proliferation, differentiation, arrest, and apoptosis; (9) ABM simulation of Ki-67 and BrdU staining; (10) ‘What-if’ analysis; and (11) model implementation and parameterization.

Geometry

To recreate the morphology of the crypt, we chose the common idealized ‘test tube’ crypt geometry of a hemisphere attached to a cylinder, which describes the basement membrane that the cells are attached to. The parameters describing the average morphology of the crypt, that is, the height and circumference of the ‘tube’, in mouse jejunum and ileum are described in Appendix 1—table 1.

Cells on the villus are terminally differentiated and can be assumed to migrate on a conveyor belt at constant velocity (Parker et al., 2017). Given these simple dynamics, to save computational power and time we modelled individual cells on the villus without spatial granularity. Cells that reach the top of the crypt are collected into a villus compartment. Shedding from the villus tip is mimicked by removing the oldest cells when the number of cells exceeds the maximum capacity of the villus, which is described in Appendix 1—table 1. Cells on the villus keep all properties and still age and undergo apoptosis if required, though in homeostatic conditions cells are usually shed into the lumen before becoming senescent.

Appendix 1—table 1
Parameter values of the agent-based model (ABM) model of the mouse small intestinal crypt.
ParameterValueSource/justification
JejunumIleum
Geometry
Ch , maximum fold-increase of the height of the crypt1.2Calibration of ABM to experimental data in drug-injured crypt
q, Hill coefficient of height scaling3.5Calibration of ABM to experimental data in drug-injured crypt
Crypt circumference9 distance units*Calibration of ABM to experimental data in a homeostatic crypt

Due to cell compression, this results in a crypt circumference of 10 cells
Villus length1000 cells (700 for cell ablation)
Biochemical signalling
BMP/Indian hedgehog feedback
z0, z coordinate at the top of the crypt in homeostatic conditions (measured from the edge of stem cell niche in the direction of the longitudinal crypt–villus z-axis)18 distance units*.16 distance units*.Calibration of ABM to experimental data in homeostatic and drug-injured crypt

Due to cell compression, this results in a homeostatic crypt height of 22–26 (jejunum) and 20–24 (ileum) measured in cell positions
z50. z coordinate at which the number of mature enterocytes becomes greater than the number of absorptive progenitors in homeostatic conditions (measured from the edge of stem cell niche in the direction of the longitudinal crypt–villus z-axis)14.5 distance units*.13 distance units*.Calibration of ABM to experimental data in a homeostatic crypt

Due to cell compression, this corresponds to cell position 18 (jejunum) and 16 (ileum)
Eh , homeostatic enterocyte count940 cells920 cells
(644 for cell ablation experiment)
Calibration of ABM to experimental data in homeostatic and drug-injured crypt (Umar, 2010)

Cell ablation experiment scaled using 0.7× reduction in villus length
p, Hill coefficient of Hedgehog/BMP feedback8.54Calibration of ABM to experimental data in a drug-injured crypt
A, level of BMP signals at height ztop64Calibration of ABM to experimental data in homeostatic crypt
B, exponential transformation of the diffusion coefficient of BMP signals8Calibration of ABM to experimental data in homeostatic crypt
Wnt signalling
kWnt152.80Calibration of ABM to experimental data in homeostatic crypt
WntRange, Wnt signalling short-range field around Paneth cells and Wnt-emitting mesenchymal cells where Wnt signals tethered to receptive cells0.35 distance units*Calibration of ABM to experimental data in homeostatic crypt
Mesenchymalniche, number of Wnt emitting mesenchymal cells surrounding the stem cell niche32 cellsAssumed equal to the total number of epithelial cells in the niche in homeostatic conditions (Wright and Alison, 1984; Snippert et al., 2010)
MWnt, maximum number of Wnt signals a cell can have tethered128Calibration of ABM to experimental data in homeostatic crypt

Value chosen to be a power of 2 to facilitate dividing Wnt signals in half upon cellular division
Notch signalling
knotch200Calibration of ABM to experimental data in homeostatic crypt
Notch range0 during homeostasis, increasing to one cell unit with drop in local cell densityCalibration of ABM to experimental data in homeostatic and drug-injured crypt
ZNRF3/RNF43 signalling
Z, maximum signal strength immediately around the emitting cell1Calibration of ABM to experimental data in homeostatic crypt
LZNRF3, length scale of ZNRF3 signalling.1 distance units*Calibration of ABM to experimental data in homeostatic crypt
K, regulates the dependence of the decay rate of cell surface-tethered Wnt molecules on ZNRF3/RNF43 signals19.1Calibration of ABM to experimental data in homeostatic crypt
ZNRF3, ZNRF3 signal strength experienced at the edge of niche3.5Calibration of ABM to experimental data in homeostatic crypt. Maintains homeostatic number of stem cells, ≈14-16 (Snippert et al., 2010)
u, Hill coefficient of ZNRF3 feedback.1Calibration of ABM to experimental data in a drug-injured crypt
Contact inhibition and mechanical parameters
p03.2Calibration of ABM to experimental data in homeostatic crypt
ϵPaneth-Paneth , adhesive constant for Paneth–Paneth interactions0.1111Calibration of ABM to experimental data in homeostatic crypt
ϵ, adhesive constant for all other interactions0.01111Calibration of ABM to experimental data in homeostatic crypt
ν, Poisson ratio0.4Based on published data (Geissler and Hecht, 1981; Mokbel et al., 2020; Mahaffy et al., 2004)
E, Young’s modulus25 kPa (Paneth), 4 kPa (all others)Pin et al., 2015
μfric.π(r2pij2)10000Calibration of ABM to experimental data in homeostatic crypt

Maintains published cell transfer velocity of 1 cell/hr (Potten, 1998)
μ1+2πrcell23000
(multiplied by 10,000 for Paneth cells)
Calibration of ABM to experimental data in homeostatic crypt

Maintains published cell transfer velocity of 1 cell/hr (Potten, 1998)
Cell timescales
tcycleshort10 hrWright and Alison, 1984
tcyclelong21.5 hrSchepers et al., 2011
Paneth cell lifespan54 dIreland et al., 2005
Other cell lifespan6 d (jejunum), 6.5 d (ileum)
(4 d for villus cells in Barker experiment)
Calibration of ABM to experimental data in drug-injured crypt
Apoptosis duration12 hrSundquist et al., 2006
Paneth ->Stem de-differentiation duration48 hr (reversible for first 36 hr)Calibration of ABM to experimental data in homeostatic crypt (Tan et al., 2021)
Absorptive progenitor ->Enterocyte differentiation duration0 hrCalibration of ABM to experimental data in homeostatic crypt
All other differentiations duration4 hrStamataki et al., 2011
Duration of Ki-67 positivity6 hours post-differentiation.
6/8/10/12 hr after drug-induced cell cycle interruption in G1/S/G2/M-phase
Calibration of ABM to experimental data in homeostatic and drug-injured crypt.(Miller et al., 2018)
KBrdU, theoretical maximum level of BrdU5.7Calibration of ABM to experimental data in homeostatic crypt
Cell fate decision parameters
Stem ->Paneth Notch threshold3Calibration of ABM to experimental data in homeostatic crypt
Absorptive/secretory Notch threshold2Calibration of ABM to experimental data in homeostatic and drug-injured crypt
Paneth ->Stem emergency Notch threshold5Calibration of ABM to experimental data in homeostatic crypt.
Qualitatively reproduces (Yu et al., 2018)
Paneth/stem Wnt threshold64Arbitrary parameterization
Cell cycle modifications
cmasscmassshort0.962Calibrated to enable variable division time responding to mechanical cues
cmasslong1.175
cVsicVsishort0.1Calibrated to modify G1-phase duration for variable division time
cVsilong1
KAWee1pKAWee1pshort2.76Calibrated to maintain S-phase duration for variable division time
KAWee1plong0.97
DNA/RNA module
v12Calibration of ABM to experimental data in homeostatic crypt
v25Calibration of ABM to experimental data in drug-injured crypt
s3Calibration of ABM to experimental data in drug-injured crypt
KVsi5Calibration of ABM to experimental data in drug-injured crypt
t2Calibration of ABM to experimental data in drug-injured crypt
Drug simulation parameters
Cdk1 inhibitor
kdrug,CycA70N/A
kdrug,CycB140N/A
Mass threshold4Mass of cell determining cycle arrest for cells unable to undergo mitosis
5-FU simulation
dDNA50.8Calibration of ABM to experimental data in drug-injured crypt
KFdUTP13100Calibration of ABM to experimental data in drug-injured crypt
n5.82Calibration of ABM to experimental data in drug-injured crypt
dRNA30.3Calibration of ABM to experimental data in drug-injured crypt
KFUTP1610Calibration of ABM to experimental data in drug-injured crypt
m1.940Calibration of ABM to experimental data in drug-injured crypt
Computational parameters
dt, timestep for movement0.0001 dN/A
dtcycle, timestep for cell cycle0.00001 dN/A
  1. *

    All distances are normalized such that an average, isolated cell has a diameter of one distance unit. We have assumed cells are deformable and hence lose the spherical shape upon compression so that the cell diameters, in both the z-axis direction (longitudinal crypt–villus axis), as well as in the y–x plane (crypt transversal circumference), are smaller than one unit and result in inequality between the number of cells and the distance units.

Cell cycle proteins and cellular growth

The division cycle of cells is controlled by a network of interacting proteins which include cyclins, cyclin-dependent kinases (CDKs), and a suite of ancillary proteins (Morgan and Morgan, 2007). The discrete events of the cell cycle, such as DNA replication in S-phase and the various stages of mitosis, are regulated by the activity of this protein network, whose components go through a careful, conserved series of peaks and troughs at the correct pace to complete all processes of the cycle. The dynamics of this protein interaction network is simulated in each cell of the ABM and controls cell division and differentiation.

We have used the model of Csikász-Nagy et al., 2006 that recreates the mammalian cell cycle and is available in Biomodels (Le Novère and Csikasz-Nagy, 2006). This model is an extension of the pioneering work of Novak and Tyson that helped reveal the complex nonlinear dynamics of the cell cycle proteins (Novak and Tyson, 1993; Novak et al., 2001; Novák and Tyson, 2004). The Csikasz-Nagy model provides multiple necessary features such as core cell cycle proteins, a mass variable that can be coupled to the volume of the single cells in our ABM and sufficient mechanistic detail to enable a detailed description of drug–cycle interactions. The model compromises 14 variables that describe the dynamics of the concentration of the main cell cycle proteins as oscillations between alternating peaks and troughs. G1-phase is the default opening state, with low levels of cyclins A, B, and E and high level of Wee1. The level of cyclin D grows exponentially throughout the cycle and is halved between daughter cells after mitosis. S-phase begins with the increase in cyclin E and ends when Wee1 drops to reach its trough. G2-phase is characterized by low Wee1 and high cyclin A, ending with the drop of cyclin A. M-phase ends when cyclin B falls and the cell divides and restarts the cycle in G1.

Stem cells have been reported to have a longer division cycle than absorptive progenitor cells (Wright and Alison, 1984; Marshman et al., 2002; Potten et al., 1997). We hypothesize that this is due to contact inhibition mechanisms caused by increased intercellular forces in the crowded, constrained niche. This implies that the duration of the cycle may significantly vary among single cells. To implement cycles of varying duration in our ABM, we describe below a series of required adjustments in the Csikasz-Nagy model that basically involve changes in the duration of the full cycle, the re-adjustment of the length of the cycle phases, primarily G1 and S-phase, and the modulation of the dynamics of the model mass variable.

To change the duration of the cell cycle, tcycle , we rescaled the time coordinate: tτtcyclet, where τ=140.027 h is the original period of the model (Csikász-Nagy et al., 2006) and tcycle is determined by the internal pressure of the cell as detailed below in Section 1.6.

Without further modifications of the Csikasz-Nagy model (Csikász-Nagy et al., 2006), the duration of all cycle phases would be scaled in proportion with changes in tcycle. However, not all phases are proportionally shortened in fast cycling healthy cells (Csikász-Nagy et al., 2006). The shorter cycle duration in absorptive progenitors is likely due to shortening/omission of G1-phase as reported for rapid cycling progenitors (Wright and Alison, 1984; Carroll et al., 2018 ), while the duration of S-phase is less variable (Wright and Alison, 1984) with reported values of 8 hr for mouse ileal epithelium (Wright and Alison, 1984).

Regarding G1-phase, the p27 protein has been reported to regulate the duration of G1 by preventing the activation of cyclin E-Cdk2 which induces DNA replication and defines the beginning of S-phase (Morgan and Morgan, 2007). We hypothesized that fast cycling cells have low levels of p27 which results in earlier DNA replication, bringing forward the start of S-phase and shortening the length of G1. In support of this hypothesis, it has been experimentally demonstrated that inhibiting p27 has no effect on the proliferation of absorptive progenitors (Zheng et al., 2008). In the Csikasz-Nagy model (Csikász-Nagy et al., 2006), the duration of G1 can be modulated through the parameter Vsi , which is the basal production rate of p21/p27 (in the Csikasz-Nagy model, the p21 and p27 proteins are represented by a single variable, here we refer to that model quantity as p21/p27).

Additionally, the end of S-phase is associated with the decrease in Wee1 to basal levels due to Cdc14-mediated phosphorylation of Wee1. In the Csikasz-Nagy model (Csikász-Nagy et al., 2006), this reaction is described by a Goldbeter–Koshland function, which includes the parameter KAWee1p to regulate the level of Cdc14 required for the phosphorylation of Wee1.

Therefore, we modified these two parameters, Vsi and KAWee1p , to ensure that variations of the cycle duration mostly impact on G1 while the length of S-phase remains constant. We assumed that the value of the two parameters scales linearly with the duration of the division cycle, tcycle , between a lower and upper bound, which prevent aberrant behaviour of the cell cycle model in the dynamically changing conditions of the crypt.

Vsi is scaled according to

VsiCVsiVsi,CVsi={CVsishorttcycle<tcycleshort,(CVsilongCVsishort)(tcycletcycleshort)(tcyclelongtcycleshort)+CVsishorttcycleshorttcycletcyclelong,CVsilongtcycle>tcyclelong,

where tcycleshort and tcyclelong denote the average duration of the cycle of fast cycling progenitors and of the slower cycling stem cells, respectively. CVsishort and CVsilong are values calibrated to ensure the correct duration of G1 for the short and long cycles, respectively, and can be found in Appendix 1—table 1.

Similarly, we scale KAWee1 using the function:

KAWee1p={KAWee1pshorttcycle<tcycleshort,(KAWee1plongKAWee1pshort)(tcycletcycleshort)(tcyclelongtcycleshort)+KAWee1pshorttcycleshorttcycletcyclelong,KAWee1plongtcycle>tcyclelong.

Here, KAWee1short and KAWee1long are the values required to maintain constant duration of S-phase in fast and slow cycling cells and can be found in Appendix 1—table 1.

A further refinement required to modify the length of the cycle in the Csikasz-Nagy model comprises the mass variable. This variable doubles its value over the course of a cycle and drives the progression of the cell cycle by changing the production rates of the cycle proteins. The changing production rates affect the balance of the proteins and the duration of the cell cycle phases, which start and end at particular mass values determined by the abovementioned two rates and other parameters in the model. After the mass doubles, mitosis occurs and the mass is halved to its initial value, returning the model to the original state. From here the mass begins to grow again, repeating the cell cycle. The mass of a cell effectively tracks the cell’s progress through the cell cycle.

In our ABM, tcycle changes continuously in each cell and modifies Vsi and KAWee1p as described above, which in turn changes the mass values of the start/end of the cell cycle phases. Without further changes in the model, this would cause the cells to not progress through the cell cycle correctly, with unbalanced phases duration and dividing at unwanted mass values, causing erroneous and unrealistic behaviour in the ABM.

This can be solved by normalizing the mass in the cell cycle model, chosen such that a cell begins at mass=massinit1 and always divides at mass=2. To do this, we first define a normalized mass variable, assumed to be proportional to the volume of the cell:

massVr3mass=2r3rfinal3,

where r is the cell radius that takes values between rinit. and rfinal. When a proliferative cell is created, it is assigned a desired final size, rfinal=23r , where rN(0.35,0.00875) for stem cells and rN(0.5,0.0125) for all other cells. The mean values, 0.5 and 0.35, of the radius of progenitor and stem cells, respectively, were determined for an average, non-proliferative or proliferative progenitor cell to have, without loss of generality, a diameter of 1 while the diameter of an average stem cell is slightly smaller, 0.7. In this way, the model captures the smaller size described for columnar LGR5+ stem cells (Barker et al., 2008), which additionally helps recapitulate the mechanics and cell composition of the niche. The variance of the radius was determined by our implementation of the cell cycle model in the ABM. In our model, the volume of the cell is equated to the cell’s mass parameter of the Csikasz-Nagy model and, hence, the cell final radius determines the duration of the cell cycle as described above. By simulating the cell cycle model, we observed that large values of the standard deviation resulted in some cells progressing through the cycle too quickly and, therefore, failing to complete the cell cycle correctly. This analysis provided an upper limit to the coefficient of variation (CV) = 0.025 to ensure all cells progress regularly through the cycle during homeostasis. This results in values of the standard deviation of the radius of 0.0125 and 0.00875 for progenitor cells and stem cells, respectively. Of note, a cell radius CV of 0.025 corresponds to a cell volume CV of about 0.075 which is not far from the reported experimental CV for cell volume, about 0.11 (Bell and Anderson, 1967).

We then introduce a factor cmass onto the four terms involving the mass variable in the cell cycle model. These terms are the basal production rates of the four cyclins A, B, D, and E, called Vsa, Vsb , CycD0, and Vse, respectively. cmass is given by

cmass={cmassshorttcycle<tcycleshort,(cmasslongcmassshort)(tcycletcycleshort)(tcyclelongtcycleshort)+cmassshorttcycleshorttcycletcyclelong,cmasslongtcycle>tcyclelong.

The values cmassshort and cmasslong are values found by calibration of the cell cycle model to guarantee the cell always divides at mass=2 for the short and long cycle durations.

Moreover, the cell mass is assumed to grow exponentially. A proliferative cell always reaches a final value of mass=2, corresponding to the radius rfinal=23r , during the cycle time, tcycle, so that mass must grow as

dmassdt=ln2massinittcyclemass.

This corresponds to a radial growth rate of

drdt=ln(rfinalrinit)tcycler.

As tcycle changes dynamically through the cell cycle, the growth rate holds only for the instantaneous conditions the cell is experiencing and changes dynamically through the cell’s lifetime. However, in a healthy crypt, extracellular conditions vary slowly, and the value of tcycle and all derived adjustment factors remain relatively unchanged.

We assumed that cells divide symmetrically. Each daughter cell has a starting radius of rinit.=rparent and is assigned with a new randomly generated r value which determines rfinal=23r. If rinit.>rfinal, then we set rfinal=rinit. to prevent values of mass>2. Since cells have a variable maximum size uncorrelated to their birth size, that is, rparentrchild, the initial mass value is not necessarily 1. Longer or shorter G1 phases emerge from the model to adjust the cycle duration in cells that begin with mass<1 or mass>1, respectively.

Proliferative daughter cells continue through its own cell cycle and proceed to grow to its own rfinal=23r. Non-proliferative secretory cells differentiate from stem cells, which are smaller than other cells. To compensate for this, secretory cells grow to reach a radius r, generated as rN(0.5,0.0125), in a time equal to tcycleshort. The other type of non-proliferative cells, enterocytes, derive from absorptive progenitors and remain at rinit.0.5 without increasing size.

These definitions of mass, cell radius, and cell growth were chosen to ensure that cells have a consistent radius and guarantee that the cell cycle model correctly proceeds through all phases in each cell. Due to the varying cycle duration and extracellular conditions, this control is essential to the correct functioning of the cell cycle and overall behaviour of the ABM.

Drug perturbations of the cell cycle model: CDK1 inhibition

We have used the Csikasz-Nagy cell cycle model to implement drug-induced perturbations of the cell cycle proteins, which are common mechanisms of action of oncotherapeutics, in our ABM. For an arbitrary component of the cell cycle model, X, we introduce a term dependent on the drug and X:

dXdtf(Drug,X),

where ⊂ means ‘contains the term’ and Drug represents the cell concentration of the active compound/metabolite which is often described by a pharmacokinetics model. fDrug,X, quantifies the effect of the drug on X. This function can take several forms such as a mass-action term or a Michaelis-Menten or Hill equation. Multiple terms like this can be added concurrently to the proteins described by the Csikasz-Nagy model.

As an example, we have modelled the effects of a Cdk1 inhibition at the single-cell level in our ABM. Cdk1 binding is reported to induce nuclear translocation of cyclins A and B require to initiate mitosis (Pesin and Orr-Weaver, 2008). Accordingly, we have added a mass-action term onto the rate of change of the CycA/Cdk1/2 and CycB/Cdk1 complexes as follows:

dCycAdt=Vsa+(Vdi+kdissa)TriAkassap27CycAVdaCycAkdrug,CycA[Drug]CycA,
dCycBdt=Vsb+(Vdi+kdissb)BCKI+V25pBkassbp27CycB(Vdb+VWee)CycBkdrug,CycB[Drug]CycB,

where CycA and CycB are used to refer to CycA/Cdk1/2 and CycB/Cdk1 to improve readability of the equation. kdrug,CycA and kdrug,CycB are parameters that quantify the drug effect, with values specified in Appendix 1—table 1, and [Drug] denotes a theoretical drug dynamical concentration. For the simulation in Figure 4, we considered a CDK1 inhibitor that was administered every 12 hr for 4 d, with active cytotoxic effects for 6 hr. To model this, [Drug] is given by the formula

[Drug]={1tdosetime<tdose+6,0otherwise,

where tdose0,12,24,36,48,60,72,84 hr. Also, we considered a smaller value for kdrug,CycA than for kdrug,CycB to reflect the fact that CycA represents both CycA/Cdk1 and CycA/Cdk2 and only CycA/Cdk1 is inhibited.

These perturbations of the cell cycle proteins can cause incorrect progression through the cell cycle, whereupon a cell is permanently arrested. A disorderly restart of the cycle, leading to enlarged cells, is observed when CDK1 inhibition prevents cells at the end of G2 from entering M-phase or induces early reduction of cyclins A and B during M-phase, with cells failing to complete cytokinesis and prematurely restarting G1. Cells in M-phase subjected to greater reductions of cyclins A and B, which completely disrupt the protein network, undergo mitotic death.

DNA and RNA synthesis

Since one of the most common means of targeting the cell cycle is to exploit the effect of DNA-damaging drugs (Helleday et al., 2008), we added the dynamics of DNA replication during S-phase and RNA synthesis during the cell cycle.

Replicating DNA is represented by two variables, DNA1 and DNA2, which denote two DNA double helices formed during S-phase. DNAi is an abstraction of the proportion of undamaged DNA, which takes values from 0, representing total DNA disruption, to 1 for the whole undamaged double helix.

At the onset of S-phase, the original DNA double helix, DNA1, unwinds to start the replication of strands and rapidly generates two complete sets of DNA, DNA1 and DNA2. This is represented in the model by

{DNA1,DNA2}={C,0}S-phase onset{C/2,C/2},

Both DNA1 and DNA2 aim to reach DNAi=1: DNA synthesis is assumed to be at a faster rate during S-phase, and outside S-phase DNA synthesis takes place solely for repair at a slower rate. Hence, in healthy cells, these variables obey the following equations and algorithm:

dDNAidt=kDNA,kDNA={0ifDNAi=1,v1if in S-phase and0<DNAi<1,v12if0<DNAi<1.

The DNA replication rate, ν1, is sufficiently fast to ensure DNAi reaches 1 during S-phase in healthy cells. Outside of S-phase, we assumed a twofold slower rate for DNA repair when the cell is not actively replicating its DNA. Values are specified in Appendix 1—table 1.

When the cell divides, the daughter cells are given one DNA double helix each (which are both assigned to DNA1 in the respective daughter cell) to restart the cycle.

RNA levels are represented by a single RNA variable. Similarly, this variable is an abstraction of the proportion of undamaged RNA in the cell, with RNA=1 in a healthy cell and RNA=0 for total RNA disruption. RNA synthesis is assumed to be governed by a simple linear-growth differential equation until its maximum value, RNA=1, and remains at this value unless damage is induced as follows:

dRNAdt=kRNA,kRNA=0if RNA=1,elsekRNA=v2,

with parameter values specified in Appendix 1—table 1.

Along with these equations for DNA and RNA levels, we added DNA and RNA-damage checkpoints to modulate the response of the Csikasz-Nagy cell cycle model to perturbations. We considered both the G1/S and the G2/M checkpoints (Morgan and Morgan, 2007), with cells checking their DNA and RNA levels as they progress from G1 to S-phase, and from G2 to M-phase. If the DNA and/or RNA levels are below the threshold values (see Appendix 1—table 1), the cell undergoes apoptosis. Checkpoint failures can occur upon drug-induced DNA or RNA damage, as explained below.

Drug perturbations of RNA and DNA synthesis: 5-FU-induced RNA and DNA damage

Similar to the cell cycle model, drug effects are represented by adding a negative term to these differential equations:

dDNAidt=kDNA-fDrug,DNAi,
dRNAdt=kRNAf(Drug,RNA),

where f(Drug,X) could be a mass action, Hill equation, or Michaelis–Menten term quantifying the drug-induced RNA or DNA damage.

DNA damage induces increased p21 expression in cells, which prevents progression through the cell cycle and can lead to cell cycle arrest or apoptosis (Abbas and Dutta, 2009). To replicate this, we further modified the p21/p27 term in the Csikasz-Nagy model to respond to the DNA levels of the cell. Recall that Vsi was the production rate of p21/p27 in the model, and we multiplied this by CVsi to moderate the production of p21/p27 (see details in ‘Cell cycle proteins and cell size/growth’ above). Recall that VsiCVsiVsi ; to replicate DNA-damaged induced production of p21, we replace CVsi with a bounded function dependent on the cell’s DNA levels

CVsiKVsi1+(KVsiCVsi1)DNAs.

DNA=DNA1 in G1, and DNA=min(DNA1+DNA2,1) in all other phases and s is a scaling coefficient. In homeostasis, with DNA1+DNA21, this function is equal to CVsi and the cell cycle model proceeds as before. With severe DNA damage, DNA1, the function is approximately equal to KVsi, always >CVsi, that represents the maximum fold increase in the production rate of p21, that is, VsiKVsiVsi. Parameter values can be found in Appendix 1—table 1. When DNA levels are reduced by drug-induced injury, this new function increases the production rate of p21/p27 which slows down the production of cyclins and the progression of the cell cycle, recapitulating a reversible cell cycle arrest for low-to-moderate DNA damage (Shaltiel et al., 2015).

Cell growth is dependent on the correct translation of mRNA into proteins. We hypothesized that RNA damage reduces a cell’s capability of biosynthesis and leads to slower cellular growth (Wurtmann and Wolin, 2009). This is modelled by adding an RNA-dependent factor to the growth rate of cells:

dmassdt=2RNAtRNAt+1ln(2massinit)tcyclemass,

where RNA takes as defined above values between 0 and 1 and t is a scaling coefficient. Parameter values can be found in Appendix 1—table 1. By linking RNA integrity to cellular growth, we allow RNA damage to induce a form of cell cycle arrest, as previously reported (Chernova et al., 1995; Bellacosa and Moss, 2003).

The result of these responses to DNA and RNA damage, in combination with the cell cycle checkpoints, allows the cells in our model to exhibit a progression of responses to increasingly severe DNA and RNA damage. Cells with slightly damaged DNA and/or RNA levels grow and proliferate slowly due to impediment of their cell cycle and/or cellular growth. With moderate DNA and RNA damage, a cell enters an impermanent, reversible cell cycle arrest (characterized by a near-zero growth rate and p21-induced halt of the cell cycle). Upon interruption of the drug-induced insult, these cells will re-enter the cell cycle. In case of severe DNA and/or RNA damage, a cell will undergo DNA/RNA damage-induced apoptosis caused by failing a cell cycle checkpoint. Additionally, drug-induced perturbations may result in incorrect progression through the cell cycle, which causes the cell to enter a permanent arrested state or die as described above. Note that though RNA damage is known to cause cell cycle arrest and apoptosis (Bellacosa and Moss, 2003), the mechanisms are poorly known, so we made the conservative decision to check the level of RNA damage at the same checkpoints as DNA damage.

As an example, we modelled 5-FU-induced RNA and DNA damage in the intestinal epithelium. We considered the two main downstream metabolites of 5-FU, FdUTP and FUTP, causing DNA and RNA damage, respectively (Longley et al., 2003). To do this, we implemented in the ABM a previously published model that describes 5-FU distribution post-dosing in mouse and a reduced version of the 5-FU metabolic pathway (Gall et al., 2023).

Furthermore, we implemented the effect of FdUTP and FUTP on DNA and RNA synthesis, respectively, on each cell of our ABM using a Hill function as follows:

dDNAidt=kDNA-dDNAiDNAiFdUTPnFdUTPn+KFdUTPn,
dRNAdt=kRNA-dRNARNAFUTPmFUTPm+KFUTPm.

Parameter values can be found in Appendix 1—table 1.

The impact of these metabolites on DNA and RNA of each cell of the epithelium resulted in the arrest of the majority of proliferative cells, with a small proportion undergoing apoptosis after failing the G1/S or G2/M checkpoint.

Mechanical cell interactions and contact inhibition

Intestinal stem cells and early progenitor cells compete for limited niche space and, therefore, the ability to retain or regain stemness. Cell proliferation creates a constant battle for space, inducing forces that drive cell migration away from the hard boundary of the stem cell niche towards the top of the crypt and onto the villus.

We assumed intercellular physical forces based on Hertzian contact mechanics with adhesive and frictional forces, similar to those in published reports (Galle et al., 2005; Buske et al., 2011). For the sake of simplicity and differently from previous approaches, we did not include the extra repulsive force opposing the reduction in cell volume caused by cell overlapping and did not consider radial expansion of cells to compensate for the loss of volume in compressed cells.

In our model, cells experience repulsive, adhesive, and frictional forces. Forces result in movement according to Stoke’s flow, where viscous forces dominate inertial forces, such that cell velocity is directly proportional to the resultant forces on the cell. For very shallow overlapping distances (10% of the cells radius), the adhesive force holds the cells together and replicate continuity of a biological tissue, but for greater overlap distances, repulsive forces dominate. Frictional forces help create collective movement by counteracting cell migration in the opposite direction to the general flow of cells.

All distances are expressed in arbitrary units (A.U.) defined such that 1 distance unit is equal to the diameter of an average, isolated cell. Forces are then measured in the resulting units. We have assumed cells are deformable and hence can lose their spherical shape when responding to mechanical forces. Regions with high proliferation result in cell diameters, in both the z-axis direction (longitudinal crypt–villus axis) and the x–y plane (crypt transversal circumference), smaller than 1 unit and, hence, in inequality between the number of cells and the distance units.

Contact repulsion

Cells are assumed to be elastic spheres with intercellular forces derived from Hertzian contact mechanics. The magnitude of the repellent force, Frep.ij between cell i (with position vector xi, radius Ri, Young’s modulus Ei, and Poisson ratio νi) and cell j (with position vector xj, radius Rj, Young’s modulus Ej, and Poisson ratio νj) is described as follows:

Frep.ij=43ER12dij32,1E=1νi2Ei+1νj2Ej,1R=1Ri+1Rj,

where dij=Ri+Rj|xij| is the overlapping distance between cells measured on the line joining the cell centres, with xij=xjxi the displacement vector joining the two cell centres. This repulsive force acts on both cells in opposing directions, pushing them away along the unit vector joining the two cells x^ij:

Frep.ij=Frep.ijx^ij.

The reported value for the Young’s modulus of Paneth cell is relatively large (Pin et al., 2015) and results in a relatively large force acting on neighbouring stem cells which helps to confine them in the niche. In addition, the previously published values of the Poisson ratio indicate that cells are marginally compressible (Geissler and Hecht, 1981).

Adhesive force

All cells in contact experience adhesive forces proportional to the area of contact and the cells inherent adhesiveness, parameterized by ϵ. The magnitude of adhesive force between cell i and j is quantified as follows:

Fadh.ij=2ϵπpij(1pij|xij|),

where |xij| is the distance between cell centres and

pij=Ri2Rj2+|xij|22|xij|.

This force is again directed along x^ij, pulling the cells together: Fadh.ij=Fadh.ijx^ij and its magnitude is derived by assuming the associated energy, Eadh.ij, is proportional to the area of contact between cells i and j, Eadh.ij=ϵAij, where Aij=πRi2-pij2, and differentiating with respect to the distance between the cells.

Two cells in isolation will be at rest when the repulsive and adhesive forces are equal; however, in our simulations, this rarely happens due to the constant proliferation and growth of surrounding cells. In vivo crypts have a highly compressed niche with tightly packed stem cells wedged between Paneth cells. In our model, the repulsive force is parameterized entirely by observed quantities (the Young’s modulus and Poisson ratio), leaving ϵ in the adhesive force as a free parameter. The value of ϵ determines intercell separation at rest. This value was chosen to allow overlapping of Paneth cells at rest of 0.15 distance units, which corresponds to 15% of the diameter of an average Paneth cell. This results in ϵ=0.216 for Paneth–Paneth adhesion. Qualitatively, all other cells are less tightly packed, so all other adhesive forces (including Paneth cells with any other cell type) are assumed to be tenfold weaker with ϵ=0.0216, which produces an overlap of approximately 0.075 cell units. These assumptions facilitate the recapitulation of the tighter packed cells in the niche, resulting in increased mechanical pressure (defined in the following sections) which induces proliferation contact inhibition mechanisms.

Frictional force

Cells that are in contact experience a frictional force proportional to their relative velocity. The force acting upon cell i due to friction with cell j is quantified as follows:

Ffric.ij=μfric.Aij(dxidtdxjdt),

where Aij is the area of contact between cells i and j defined above, and μfric. is a numerical constant calibrated to enforce orderly cell dynamics. This force is comparatively smaller than the other forces but helps the collective motion of cells by opposing cell migration against the common direction.

Cell migration

Under a force, cells move according to Stoke’s flow, where viscous forces are assumed to dominate over inertial effects:

md2xidt2=jFijμdxidtdxidt=1μjFij.

Therefore, the position vector of the i-th cell, xi , is updated according to

Δxi=1μjFijΔt,

where Fij=Frep.ij+Fadh.ij+Ffric.ij is the resultant of all forces on cell i due to cell j.

The parameter μ links the forces to cellular motion. The value of this parameter is estimated to recapitulate the transfer velocity in the crypt–villus junction measured in in vivo experiments to be approximately one cell position per hour in mice (Potten, 1998). However, cell motion response to these forces may vary for different cell types. It has been reported that Paneth cells persist in the stem cell niche at the crypt base for relatively long periods of up to 57 d in mice (Ireland et al., 2005; Roth et al., 2012) and exhibit elevated β4 -integrin expression anchoring them to the mesenchyme (Langlands et al., 2016). Additionally, Paneth cells are larger and stiffer than the comparatively malleable stem cells which suggest that they require greater forces to be displaced. In our model, we used μ to replicate this behaviour and recreate drag effects of the basal membrane/mesenchyme. We implemented a value of μ for Paneth cells 10,000-fold greater than for other cells, effectively making Paneth cells difficult to move by other cells but allowing them to slowly move one another to form an orderly niche over longer timescales.

Internal pressure and contact inhibition

The forces described above are used to calculate the internal pressure experienced by cells, which varies according to the cell-intrinsic properties and local environment, that is, a stem cell in the crowded niche has higher internal pressure. Cell pressure is used to recapitulate contact inhibition by modulating the duration of the division cycle which increases when cells are densely squeezed together and decreases if cell density falls to enable, for instance, fast recovery from injury.

A cell feels internal stress from the surrounding cells, and this is used to simulate contact inhibition. To do this we use the concept of virial stress outlined in Van Liedekerke et al., 2015. The stress tensor for cell i, σi, is defined as follows:

σi=1VcelljFijrij,

where rij is the vector from the centre of the cell i to the plane of contact with cell j, always assumed to be on the surface of cell i, and ⊗ is the tensor/outer product combining two vectors into a ‘matrix’. Using this stress tensor, we extract the pressure in the conventional manner:

p=13tr(σ).

As all our forces are normal to the plane of contact, this reduces to

p=14πRi2j|Fij|.

This provides a rough, first-order approximation to the pressure experienced at the centre of the cell that is straightforward to compute and essential to implement contact inhibition in proliferative cells. Note that we do not consider the hydrostatic pressure induced by cell compression.

On the other hand, physical compression has been reported to lead to YAP inactivation, retarding growth and morphogenesis in the GI epithelium (Halder et al., 2012; Aragona et al., 2013; Low et al., 2014). We used our estimate of pressure to implement this contact proliferation inhibition mechanism responding to environmental mechanical cues and described the increase in the cell cycle duration, tcycle., as pressure, p, increases using a scaled logistic function as follows:

tcycle=g1+e2p0-p+tcycleshort.

Here, p0 is the average pressure experienced by cells in the niche; tcycleshort is the average division time of absorptive progenitors; and g=2tcyclelong-tcycleshort, where tcyclelong denotes the longer division time of a stem cell in average niche conditions.

This function captures the variation of the duration of the division cycle from a minimum to a maximum value in highly compressed cells which leads to longer division times in the tightly constrained stem cell niche of the crypt, while the cycle is shorter in the less compressed transit-amplifying zone, in agreement with experimental reports (Wright and Alison, 1984; Bach et al., 2000; Schepers et al., 2011).

Biochemical signalling

Next, we detail how the cells interact with one another, communicating the local composition of the crypt to maintain homeostasis through simulated biochemical signalling.

To achieve stable crypt cell composition and structure, we have implemented five signalling mechanisms including Wnt, Notch, and BMP pathways which have been demonstrated to be essential for morphogenesis and homeostasis of the intestinal crypt (Gehart and Clevers, 2019; Fevr et al., 2007; VanDussen et al., 2012; Pellegrinet et al., 2011; He et al., 2004). We have modelled contact proliferation inhibition mediated by the YAP-Hippo signalling pathway responding to mechanical forces (Halder et al., 2012; Aragona et al., 2013; Low et al., 2014) as described above and following experimental evidence (Hao et al., 2012; Koo et al., 2012; Farin et al., 2016), implemented a ZNRF3/RNF43-like mediated feedback mechanism between Paneth and stem cells.

These minimal signalling mechanisms were chosen because a full understanding of the protein interaction networks is still a topic of active research. However, even with our conservative assumptions, we implicitly introduce crosstalk between the different signalling pathways. For example, the nature of cell fate decisions leads to interaction between Wnt and Notch levels, and changes in the duration of the cell cycle caused by contact inhibition regulate the ability of a cell to accumulate signalling molecules.

Wnt signalling

The Wnt pathway is the primary pathway associated with stem cell maintenance and differentiation in the intestinal crypt as well as in many other tissues (Fevr et al., 2007; van der Flier and Clevers, 2009; Nusse and Clevers, 2017). Two sources of Wnt signals have been described in the mouse crypt: Paneth cells (Sato et al., 2011) and specific mesenchymal cells surrounding the stem cell niche at the crypt base (Stzepourginski et al., 2017).

We did not consider the dynamics of the canonical Wnt signalling molecular cascade but directly implemented downstream cellular responses to Wnt levels. We modelled Wnt signalling as a short-range field around Paneth cells and Wnt-emitting mesenchymal cells at the bottom of the crypt, acting within a distance WntRange from the surface of these cells (see Appendix 1—table 1 for value). Receptive cells within this range tether Wnt signals to their surface as previously reported (Farin et al., 2016; Clevers and Nusse, 2012). This is described by the following equation:

dWntdt={kWntIncWntdWnt(ZNRF3)WntWnt<MWnt,dWnt(ZNRF3)WntWntMWnt.

The variable ‘Wnt’ is an abstraction of the total number of Wnt ligands tethered to the surface of the cell. kWnt is the rate of Wnt signal tethering by a receptive cell and dWnt(ZNRF3) is the decay rate of Wnt signal tethered molecules. dWnt depends on the turnover of Wnt receptors assumed to be regulated by RNF43 and ZNRF3 ligands produced by stem cells, which forms a Wnt-mediated negative feedback loop as described below. MWnt describes the maximum number of Wnt signals a cell can have tethered and its value is chosen to be a power of 2 to facilitate dividing Wnt signals in half upon cellular division. IncWnt is the total amount of Wnt signal sources within range of the cell and is calculated as follows:

IncWnt=Paneth in range+MesenchymalnicheCellsniche.

Mesenchymalniche represents the number of Wnt-emitting mesenchymal cells surrounding the niche, which we assume is equal to the total number of epithelial cells in the niche in homeostatic conditions (Wright and Alison, 1984; Snippert et al., 2010). Additional Wnt production by Paneth cells is required to support the homeostatic number of stem cells in homeostasis. In the presented modelling scenarios, we assumed constant exogeneous Wnt source, that is, constant Mesenchymalniche, shared by all cells in the niche and enhancing niche recovery after damage. For instance, with lower number of cells in the niche, the survival cells will receive stronger mesenchymal Wnt signalling that enhances proliferation and recovery after perturbations. We assumed that surface-tethered signals are equally distributed between daughter cells upon cell division (Gehart and Clevers, 2019; Farin et al., 2016), so that cells eventually lose Wnt signals and their capacity to proliferate if not within the range of a Wnt source. These assumptions are partly supported by observed in vivo and in vitro behaviour, where the mesenchymal and Paneth cell-derived Wnt sources are mutually redundant (Farin et al., 2012).

ZNRF3/RNF43 signalling

In our model, Paneth cells enhance their own production by generating high Wnt local environments (van Es et al., 2005). In addition, due to their high Young’s modulus, Paneth cells create a region of high intercellular forces on neighbouring cells which leads to prolonged division times with greater opportunity for Wnt accumulation. This, in turn, expands the niche region with high Wnt and high cell pressure, promoting further differentiation into stem and Paneth cells. Therefore, without a negative feedback mechanism in our model, these features would result the expansion of the niche with stem and Paneth cells occupying the entire crypt. Additionally, two recent studies have demonstrated the existence of a negative feedback loop mediated by RNF43 and ZNRF3 ligands produced by stem cells (Hao et al., 2012; Koo et al., 2012). These studies proposed that RNF43 and ZNRF3 inhibit Wnt signalling by promoting the turnover of Wnt receptors such as Frizzled and LRP5 (de Lau et al., 2011), and showed that simultaneous deletion of these two receptors results in the formation of adenomas comprising mostly stem and Paneth cells (Koo et al., 2012).

We assumed that ZNRF3/RNF43 (henceforth called ZNRF3 for simplicity) is a diffusing, decaying signal secreted by stem cells. Without explicit knowledge of the chemical and physical properties of ZNRF3 signalling, this process is assumed to immediately reach steady state at the timescale of cellular decisions. Therefore, the ZNRF3 signal strength, ZNRF3(r), received by a cell at position r from a stem cell located at position R, is described by the diffusion equation as follows Crank, 1975:

ZNRF3(r)=Ze|Rr|LZNRF3,

where Z represents the maximum signal strength immediately around the emitting cell and LZNRF3 determines the spatial scale of diffusion, which we assume is equal to the length of a cell in order to maintain high signalling levels primarily in the niche.

The total ZNRF3 signalling received by a cell at position r is calculated, therefore, as the sum of the signal received from all stem cells:

ZNRF3(r)total=stem cellsZe|Rir|LZNRF3,

where Ri is the position of the i-th stem cell.

The strength of ZNRF3 signalling received by a cell is proportional to the number of stem cells in the immediate vicinity of the cell: in typical, homeostatic conditions, ZNRF3(r)total is high in the niche, falling off exponentially as a cell moves towards the villus.

The ZNRF3 signalling level detected by a cell, located at position r, regulates the decay rate of its surface-tethered Wnt molecules, dWnt as follows:

dWnt(ZNRF3)=K1+(ZNRF3ZNRF3(r)total)u,

where u is a scaling coefficient, and K and ZNRF3 are constants calibrated to maintain the size of niche at its homeostatic level. In particular, ZNRF3 is determined by the homeostatic number of stem cells in the niche (Snippert et al., 2010), while K was calibrated to produce a Wnt decay rate high enough to prevent Wnt values ≥64 in cells located at the edge of the niche when the number of stem cells is excessive such that ZNRF3(r)totalZNRF3. These considerations prevent the expansion of the niche by preventing cells from differentiating into the Paneth or stem cell fate (which requires Wnt64) when a cell is outside the niche.

With this implementation of ZNRF3-mediated negative feedback, the Wnt decay rate within the niche is high but is compensated by the abundant Wnt supply from mesenchymal and Paneth sources, while the Wnt decay rate rapidly drops to zero outside the niche. This means that the degradation of Wnt outside the niche has little impact on a healthy crypt and the Wnt gradient in our model is mainly generated by the halving of the surface bound Wnt signals between daughter cells upon division. Growth and proliferation derived forces drive migration of cells towards the villus while the amount of tethered Wnt decreases after each division. This recreates the observed (Farin et al., 2016) decreasing gradient of Wnt signals moving up the crypt (Figure 1), with the highest values in the niche, intermediary values in the transit-amplifying zone, and low levels in the upper crypt region of differentiated enterocytes.

The stem cell-mediated negative feedback loop regulating Wnt signalling, together with the differentiation rules described below, ensures the maintenance of the niche size and crypt composition in homeostasis. In addition, it also facilitates crypt recovery as stem cells in low numbers are able to reach greater surface-tethered Wnt levels to pass to their offspring which, in turn, can more readily acquire the required amount of Wnt to become stem cells.

Notch signalling

Active Notch signalling requires direct membrane contact between two cells, one expressing Notch ligands and the other Notch receptors (Gehart and Clevers, 2019; Pellegrinet et al., 2011; Baron, 2003; Sancho et al., 2015). In the intestinal epithelium, Notch ligands present in secretory cells bond to transmembrane notch receptors of stem cells to induce a transcriptional cascade which blocks differentiation of stem cells into the secretory lineage in a process known as lateral inhibition and leads to checkerboard/on-off pattern of Paneth and stem cells in the niche (VanDussen et al., 2012; Chen et al., 2017). With these considerations, Notch signalling, Notch, is implemented in each cell according to the following equation:

dNotchdt=kNotch(IncNotchNotch),

where IncNotch is the amount of incoming notch ligands to the cell which we assumed is equal to the number of ligands expressing cells in contact with the cell. At steady state, a cell’s Notch value corresponds to the number of incoming Notch ligands the cell is receiving: for example, a stem cell receiving Notch from one single neighbouring cell reaches equilibrium with Notch=1. The factor kNotch denotes the rate of Notch accumulation and has a relatively high value to ensure that the equilibrium is reached before the fate-commitment point at the end of G1. As described in the cell cycle section, the duration of G1 changes with the length of the overall division cycle: shorter cycles have a shorter G1 phase, shortening the time the cell has to receive Notch signals before deciding whether to differentiate or divide. Additionally, kNotch is also the decay rate of the cell Notch signalling and this relatively fast rate means that Notch must be constantly supplied for a stem cell to maintain stemness.

A reduction in cell density (e.g. by ablation of cells) can introduce gaps in the simulated epithelial tissue. In real tissues, these gaps would be covered by expansion-flattening of surviving cells to restore epithelial integrity and contact to neighbouring cells. These new contacts would allow cells to exchange Notch ligands. In our model, we do not explicitly consider the expansion of cells to fill gaps in the epithelium; however, we simulate this effect by allowing a cell to pass Notch signals to receiving cells within a larger range (one cell diameter) following a drop in local density. This allows our model to recreate the correct recovery response following ablation of cells.

BMP signalling

The Wnt gradient in the crypt is opposed by a gradient of BMP generated by mesenchymal telocytes, which are especially abundant at the villus base and provide a BMP reservoir, and by the recently identified trophocytes located just below crypts and secreting the BMP antagonist Gremlin1 (McCarthy et al., 2020). BMP signals inhibit cell proliferation and promote terminal differentiation (Qi et al., 2017). Large levels of BMP at the crypt–villus junction prevent proliferative cells from reaching the villus (Beumer et al., 2022). BMP signalling has been reported to be modulated by matured epithelial cells on the villus via hedgehog signalling (Büller et al., 2012; van den Brink et al., 2004) such that a decrease in villus cells decreases BMP signalling in the crypt, which enhances proliferation and expedites villus regeneration.

We propose a simple model that assumes that enterocytes, E, secrete diffusing signals, which could be interpreted as Indian hedgehog, to regulate BMP secretion by mesenchymal cells. The explicit pathways and associated timescales involved in BMP signalling are unknown; therefore, similar to our implementation of ZNRF3 signalling, this process is assumed to instantaneously reach steady state at the timescale of cellular decisions. As before, we assume that BMP is a diffusing, decaying signal in steady state (Crank, 1975) described by

BMP(z,E)=f(E)eln(B)(ztop(E)z)(ztop(E)z50),

where z is the position coordinate corresponding to the crypt–villus longitudinal axis; in our model z0 for cells located in the stem cell niche while z>0 for crypt cells outside the niche; ztopE is the value of z at the top of the crypt, which depends on the number of enterocytes on the villus; B is the exponential transformation of the diffusion coefficient. To facilitate the use of the model for different species, the z coordinate is standardized using z50, which is the crypt axis position at which the number of mature enterocytes becomes greater than the number of absorptive progenitors. As mentioned above, mesenchymal cells surrounding the niche secrete BMP antagonists (McCarthy et al., 2020), and we assumed that BMP signalling is effectively blocked in the niche such that BMP=0, which is approximately true for the above formula. fE describes the relationship between the number of enterocytes and maximum BMP signal intensity using an increasing Hill function:

f(E)=2A1+(EhE)p,

where Eh is the homeostatic number of enterocytes determined by in vivo experiments, p is the Hill interaction coefficient, and A denotes the level of BMP signals at position ztop. In our model, absorptive progenitors differentiate into enterocytes when BMP>Wnt, representing that the anti-proliferative BMP signalling received by the cell is sufficient to overcome the proliferative effect of Wnt (He et al., 2004). We achieved a homeostatic crypt cell composition with values of A and B that allow progenitors cells to divide in a healthy crypt at least three times before differentiating. Differentiation occurs when the Wnt content of a cell, at position z, reaches values below BMPz when migrating towards the villus.

In addition, these equations describe a frequently reported feedback response to villus injury consisting of enhanced proliferation within hypertrophic crypts (Büller et al., 2012; Pont and Yan, 2018; Sprangers et al., 2021). In our model, when the number of enterocytes on the villus falls below the homeostatic level, the production of BMP signals decreases and makes it possible for absorptive progenitors to divide more times and reach higher positions in the crypt before becoming terminally differentiated. Concurrently, the height of the crypt must increase to provide sufficient space for the extra proliferative cells. We modelled the enlargement of the crypt height responding to villus injury by varying the maximum z-coordinate of the crypt, ztop, using a decreasing Hill function as follows:

ztop(E)=Chz01+(Ch1)(EEh)q,

where z0 is the calibrated homeostatic value of ztop, Ch is the maximum fold increase in the height of the crypt, and q the Hill interaction coefficient. We do not consider cases in which the number of enterocytes on the villus increases above homeostatic levels, such that if E>Eh then ztopE=z0.

The standard manner to report the height of cells along the crypt–villus axis is in terms of cell positions, which is related to but not equal to ztop. This is because cell positions are counted from the bottom of the niche (and we defined z=0 to be the top of the niche), and that in our model the cells are squeezed together, causing the height of the crypt measured in cell positions to be larger than ztop.

Cell fate: Proliferation, differentiation, arrest, and apoptosis

In the sections above, we have outlined the dynamics of signalling pathways, cell cycle proteins, and mechanical forces. These processes interact with each other to maintain epithelial homeostasis by precisely tuning cell proliferation, differentiation, and migration within the crypt geometry.

An overall picture integrating the rules governing cell fate decision is described in Figure 1. Wnt levels ≥64 A.U. are required for stemness maintenance. For a stem cell, lateral inhibition is repressed when Notch < 3 A.U., equivalent to less than three secretory cells in the local neighbourhood. If Notch is repressed (<3 A.U.) and Wnt > 64 A.U., stem cells differentiate into Paneth cells. Paneth cells generate Wnt signals which enhance the production of stem cells and of Paneth cells themselves. Niche expansion is modulated by the ZNRF3/RNF43-mediated negative feedback mechanism (Hao et al., 2012; Koo et al., 2012; Farin et al., 2016) that makes Wnt > 64 unobtainable after reaching the homeostatic number of stem cells. Furthermore, the duration of the division cycle is dependent on local forces experienced by the cell. Cells under high mechanical pressure (in the niche) are subjected to YAP-Hippo-regulated contact inhibition and with longer cycles accumulate more Wnt and Notch signals. On the other hand, cells located outside the niche exhibit shorter cycles and cannot effectively accumulate enough Wnt signals to become stem or Paneth cells.

Stem cells with decreased levels of Wnt signalling (<64), usually located outside the niche, differentiate into absorptive proliferating progenitors if Notch signalling is active or into secretory progenitors if Notch signals <2 A.U. This lower Notch threshold value is required to maintain the correct balance of absorptive and secretory cells outside the niche in the absence of large numbers of Notch secreting Paneth cells. All cells migrate towards the crypt mouth driven by proliferation forces. During this migration, the Wnt content in absorptive progenitors is halved in each division and, away from Wnt sources, progressively decreases, while BMP signals increase, towards the villus. In our model, differentiation into enterocytes occurs when progenitors encounter a BMP signal level higher than their Wnt signal content. For instance, in the ileal crypt in homeostasis, this occurs approximately at cell position 16 from the crypt base, where progenitors migrating from the stem cell niche reach a reduced content of Wnt signals of about 8 A.U. On the other hand, the BMP signalling level has a maximum value of 64 at approximately cell position 23 from the crypt base, where BMP signals are generated by mature enterocytes. These BMP signals diffuse towards the crypt base and, hence, decrease exponentially to reach values of 8 A.U. at approximately position 16, which enables differentiation into enterocytes. Epithelial injuries resulting in a decreased number of enterocytes reduce BMP signal production and its diffusion range which results in the enlargement of the proliferation compartment as cells encounter the required level of BMP signals for differentiation only at higher positions in the crypt.

All fate decisions are assumed to be made at the restriction point which in our model is located at the end of G1 (Blomen and Boonstra, 2007). At the restriction point, cells assess their internal Wnt and Notch levels, and if these values fulfil the criteria to differentiate, they enter a quiescent state or G0, otherwise they proceed to S-phase and become irreversibly committed to complete the cell cycle of variable duration depending on local forces. This quiescent state lasts for 4 hr for all differentiating cells, except for absorptive progenitors, which differentiate straightaway into enterocytes. In accordance with Stamataki et al., 2011, a secretory progenitor requires an additional 4 hr to fully mature into a goblet or enteroendocrine cell.

Therefore, quiescent stem cells located above the fourth cell position from the crypt base (Gehart and Clevers, 2019; Potten et al., 1978; Sangiorgi and Capecchi, 2008) emerge naturally in the model as stem cells migrate outside the niche and pause the cycle to give rise to non-proliferative secretory progenitors, which have been identified with quiescent stem cells (Buczacki et al., 2013; Clevers, 2013a). Features and behaviours of these cells could be expanded if of interest for the model application.

Cell fate decisions are reversible; a stem cell that leaves the niche and differentiates into a progenitor cell can relatively quickly become a stem cell again if regaining enough Wnt signals by being pushed back into the niche. This plasticity extends to all cells: all progenitors and fully differentiated cells can revert to stem cells when exposed to sufficient levels of Wnt and Notch signals, replicating injury recovery mechanisms observed in the crypt (Hageman et al., 2020; Tetteh et al., 2016). We have assumed that all cells, except Paneth cells, need to acquire and maintain high levels of Wnt signals (>64) over 4 hr to complete the process. Dedifferentiating cells shrink to their new smaller size during the process if required.

Notch signalling mediates the process of Paneth cell de-differentiation into stem cells to regenerate the niche as previously reported (Mei et al., 2020; Yu et al., 2018). Paneth cells not supplying Notch ligands for 12 hr to recipient cells dedifferentiate into stem cells in a process that takes 36 hr to complete in agreement with published findings (Yu et al., 2018).

Additionally, Paneth cells in low Wnt conditions (e.g. a Paneth cell that is forced out of the niche) for 48 hr will also dedifferentiate into a stem cell, which with low Wnt content rapidly becomes a secretory or absorptive progenitor.

Additionally, injured proliferative cells can experience cell cycle arrest and apoptosis, induced by drug injury or by natural senescence. In arrested and apoptotic proliferative cells, the production rates of the cyclins (Vsa, Vsb, CycD0, and Vse) are set to 0 to interrupt the cell cycle. We assumed that cells remain arrested until they are shed from the villus tip or reach the end of their lifespan and become apoptotic. Apoptotic cells shrink and die with a negative linear rate of

drdt=rapop.tapop.,

where rapop. is the radius at the onset of apoptosis and tapop. is the time for the completion of apoptosis.

ABM simulation of Ki-67 and BrdU staining

This section discusses the implementation of the Ki-67 and BrdU staining simulations, which can be found in Figures 3 and 5, and is discussed in the ‘Results’ section.

For the Ki-67 staining simulation, we considered that a cell is Ki-67 positive if it is going through S-, G2-, or M-phase of the division cycle. Daughter cells are considered Ki-67 positive, regardless of their fate, during the first 6 hr after cell division. This assumption recapitulates the time reported for the Ki-67 protein to decay below detectable levels after exiting the cycle (Miller et al., 2018) and the detection of Ki-67 in G1 in continuously cycling cells (Sobecki et al., 2017). Similarly, cells are assumed to remain Ki-67 positive for 6–12 hr after drug-induced cell cycle interruption depending on the phase the cell was in upon interruption in our simulations, which recapitulates a previously published report (Miller et al., 2018), where cells exhibited greater Ki-67 levels in later cell cycle phases. In particular, cells arrested during G1, S, G2, and M phases are Ki-67 positive for 6, 8, 10, and 12 hr after arrest, respectively.

For the BrdU staining simulation, we assumed that cells become BrdU positive by effectively incorporating BrdU into their DNA when they are in S-phase or enter S-phase during the BrdU exposure window, which is considered to last 2 hr after BrdU administration in agreement with previous experimental reports (Parker et al., 2017). The initial level of BrdU after dosing in each cell is quantified by

BrdU=Round(KBrdUTBrdU2)

where Round(X) is the function that rounds X to the nearest integer, KBrdU is the theoretical maximum level of BrdU a cell can incorporate, and TBrdU is the remaining BrdU exposure time. For cells that enter S-phase after BrdU administration, TBrdU is equal to the remaining BrdU exposure time, while for cells already in S-phase at the time of BrdU dosing, TBrdU is equal to the exposure window or, alternatively, to the remaining duration of S-phase if this is shorter than the exposure window. Furthermore, we considered that a cell is BrdU positive if its BrdU level is >0 and, if dividing, the two daughter cells are given a BrdU value of BrdUdaughter = BrdUparent-1. This consideration recapitulates experimental reports indicating that the BrdU cell content is diluted in each division and it is no longer detected after 4–5 generations (Wilson et al., 2008).

The spatial data from the Ki-67 and BrdU staining experiments comprises the proportion of positive cells at each cell position by aggregating spatial counts from 20 to 50 one-dimensional longitudinal strips running from the crypt base to the villus (Parker et al., 2017; Williams et al., 2016). Therefore, cell position is reported in a one-dimensional space and measured as the cell count from the base of the crypt to the cell itself. To match these observations, we have implemented an algorithm that slices longitudinally the simulated crypts to generate 100 one-dimensional strips which are aggregated to estimate the proportion of stained cells at each position.

Furthermore, we estimated a 95% confidence interval, based on experimental error, around the simulated spatial profiles of Ki-67- and BrdU-positive cells by assuming that the proportion of stained cells follows a beta distribution with parameters α and β, PBeta(α,β). These parameters are estimated as follows:

α^=p21pe21p,β^=p(1p)1pe21p,

where p is the simulated proportion and e its standard error. We used an estimate of the standard error, e, derived from experimental data. We first studied the relationship between the mean value and the standard deviation of the proportion in three replicated control samples. The experimental data suggest that the error is lower for extreme values of p, that is, around 0 or 1, and larger for values of p around 0.5 (Appendix 1—figure 3). Thus, we described this relationship with a quadratic expression:

e(p)=a0+a1p+a2p2,

where a0, a1 and a2 are the coefficients determined from the replicated control samples and their values are displayed in Appendix 1—figure 3.

Appendix 1—figure 1
Simulated cell lineages in the crypt (A) and villus (B) during a 4-day CDK1 inhibition treatment for 6 hr, every 12 hr for four consecutive days and following recovery in the agent-based model (ABM) as described in Figure 4.
Appendix 1—figure 2
Simulated (lines) and observed (symbols) number of cells in the crypt.

(A) and villus (B) during the administration of 50 mg/kg of 5-fluorouracil (5-FU) twice a day for 4 d and following recovery as described in Figure 5. (C) Predicted (continuous lines) and observed (dashed line) proportions of Ki-67-positive cells along the crypt axis at 6 hr, 1 d, 4 d, and 6 d following the administration of 20 mg/kg of 5-FU twice a day for 4 d. Shadows depict the 95% confidence interval of our simulated staining results assuming that the proportion of staining cells has a beta distribution and estimating its error from experimental data. (D, E) Predicted (lines) and observed (symbols) number of cells in the crypt (D) and villus (E) during the administration of 20 mg/kg of 5-FU twice a day for 4 d and subsequent recovery. Vertical bars represent dosing times. Symbols represent cell counts from individual mice.

Appendix 1—figure 3
Relationship between the mean and the standard deviation of the proportion of Ki-67- (A) and BrdU-positive (B) cells observed at several crypt positions in three replicated control experiments.

What-if analysis

We investigated the effect on the simulated crypt of increasing and decreasing the strength of the main signalling pathways, Wnt, BMP, and ZNRF3/RNF43 signalling, and modifying the Notch thresholds. For each alternative parameterization, except when decreasing ZNRF3/RNF43 signalling, the simulation was run for 30 d to ensure stability was reached with the new parameter set and the final 10 d were included in the analysis. When decreasing ZNRF3/RNF43 signalling, we simulated 60 d to demonstrate the expansion of the niche and analysed the final 10 d. The reference parameter set used as baseline was the ileal mouse crypt parameter set reported in Appendix 1—table 1. In all cases, we only consider modifications of one signalling mechanism at a time.

To study alternative Wnt signalling scenarios, we used the WntRange parameter (Appendix 1—table 1) to double and halve the spreading area of Wnt signals emitted by Paneth cells while we maintained the original WntRange value for Wnt-emitting mesenchymal cells at the bottom of the crypt (Section 1.7.1; Appendix 1—figure 4A–F). When WntRange was doubled, we observed increased number of stem and Paneth cells in a noticeably enlarged niche (Appendix 1—figure 4C and D), with cells choosing the stem cell fate instead of differentiating into absorptive progenitors. On the other hand, decreasing Wnt signalling, by halving WntRange in Paneth cells but maintaining its homeostatic value in mesenchymal cells, resulted in no apparent changes in niche cell composition (Appendix 1—figure 4E and F), which resembled published experimental results of persisting functional stem cells after Paneth cell ablation (Durand et al., 2012).

Appendix 1—figure 4
Simulated ileal mouse crypt with modified Wnt signalling.

Three-dimensional image (A) and cell composition (B) by position in homeostasis and likewise after doubling (C, D) and halving (E, F) Paneth cell generated Wnt signals, maintaining homeostatic levels of Wnt of mesenchymal sources.

The ZNRF3/RNF43-mediated negative feedback mechanism regulates the size of the niche by modulating Wnt signalling. We simulated increasing and decreasing the strength of the ZNRF3/RNF43 by doubling and halving, respectively, the parameter Z described in Section 1.7.2; Appendix 1—figure 5A–F. Following the increase in the intensity of ZNRF3/RNF43 signalling, we observed a decrease in the number of stem and Paneth cells, together with relatively minor changes in the transit-amplifying region (Appendix 1—figure 5C and D). On the other hand, when decreasing ZNRF3/RNF43 signalling levels, the niche expanded, resulting in a crypt dominated by Paneth and stem cells (Appendix 1—figure 5E and F) which replicates reported experimental phenotypes (Koo et al., 2012).

Appendix 1—figure 5
Simulated ileal mouse crypt with modified ZNRF3/RNF43 signalling.

Three-dimensional image (A) and cell composition (B) by position in homeostasis and likewise after doubling (C, D) and halving (E, F) ZNRF3/RNF43 signalling strength.

To modify Notch signalling, we increased and decreased by 1 A.U. the Notch threshold required for lateral inhibition (Appendix 1—figure 6A–F). This Notch signalling threshold determines the number of contacting Notch-secreting cells (secretory lineage) required to inhibit the differentiation of stem cells into the secretory lineage. Thus, increasing this Notch threshold enhances the production of secretory cells leading to the increase in Paneth, goblet, and enteroendocrine cells (Appendix 1—figure 6C and D). Alternatively, decreasing the Notch threshold enhances differentiation into the absorptive lineage, reducing the number of Paneth and secretory cells (Appendix 1—figure 6E and F).

Appendix 1—figure 6
Simulated ileal mouse crypt with modified Notch signalling.

Three-dimensional image (A) and cell composition (B) by position in homeostasis and likewise after increasing (C, D) and decreasing (E, F) the Notch threshold by 1 A.U.

We modified the range of diffusion of BMP signals by doubling and halving the parameter A (Appendix 1—figure 7A–F) which denotes the amount of diffusing BMP signals, and hence affects the diffusion range, towards the base of the crypt (Section 1.7.4). When we increased the BMP signalling range, enterocytes differentiated at lower crypt positions, effectively reducing the transit-amplifying zone (Appendix 1—figure 7A and B). Decreasing BMP signalling strength by halving A resulted in the increase in proliferative absorptive progenitors, which reach higher positions in the crypt (Appendix 1—figure 7C and D). The niche was largely unaffected in both cases (Appendix 1—figure 7E and F).

Appendix 1—figure 7
Simulated ileal mouse crypt with modified BMP signalling.

Three-dimensional image (A) and cell composition (B) by position in homeostasis and likewise after a fourfold increase (C, D) and a fourfold decrease (E, F, respectively) of BMP signalling strength.

Model implementation and parameterization

The model is implemented using the Julia programming language. The mechanical forces, cellular motion, and biochemical signalling are simulated with a fixed timestep of dt=0.0001 d, while the proteins of the cell cycle model are simulated with a timestep of 0.00001 d. Parameter values and means used for their identification are detailed in Appendix 1—table 1.

Data availability

The current manuscript is a computational study. No data have been generated for this manuscript. Modelling code is uploaded as Source code 1.

References

  1. Book
    1. Crank J
    (1975)
    The Mathematics of Diffusion
    Clarendon Press.
  2. Book
    1. Morgan D
    2. Morgan DO
    (2007)
    The Cell Cycle: Principles of Control: OUP/New
    Science Press.
    1. Potten CS
    (1998) Stem cells in gastrointestinal epithelium: numbers, characteristics and death
    Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 353:821–830.
    https://doi.org/10.1098/rstb.1998.0246
    1. Sundquist T
    2. Moravec R
    3. Niles A
    4. O’Brien M
    5. Riss T
    (2006)
    Timing your apoptosis pathway
    Cell Notes 2006:18–21.
  3. Book
    1. Wright NA
    2. Alison MR
    (1984)
    The Biology of Epithelial Cell Populations
    Oxford: Clarendon Press.
    1. Zhang M
    2. Zhang L
    3. Hei R
    4. Li X
    5. Cai H
    6. Wu X
    7. Zheng Q
    8. Cai C
    (2021)
    CDK inhibitors in cancer therapy, an overview of recent development
    American Journal of Cancer Research 11:1913–1935.

Decision letter

  1. Mariana Gómez-Schiavon
    Reviewing Editor; Universidad Nacional Autónoma de México, Mexico
  2. Aleksandra M Walczak
    Senior Editor; École Normale Supérieure - PSL, France
  3. Mariana Gómez-Schiavon
    Reviewer; Universidad Nacional Autónoma de México, Mexico

Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "Homeostasis, injury and recovery dynamics at multiple scales in a self-organizing intestinal crypt" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, including Mariana Gómez-Schiavon as Reviewing Editor and Reviewer #1, and the evaluation has been overseen by Aleksandra Walczak as the Senior Editor.

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Recommendations for the authors:

– We recommend revising the writing in the introduction and the abstract to clarify the background, motivation, impact, and context of the research. We suggest that the authors explicitly state why studying the intestinal crypt is broadly important to the scientific community and why a mathematical model of it is necessary for progress in the development of oncotherapies. Further, we recommend highlighting the impact of their model of the intestinal crypt on the development of oncotherapies. Most of this information is included in the Discussion section, but not in the introduction or the abstract.

– ABM must be defined in the main text, not just the abstract.

– The authors claimed to demonstrate the model's application to epithelial research and drug discovery. They evaluate some known drugs, but we don't believe the term "discovery" is fully justified in the manuscript.

– We suggest that the authors specify how the model they developed compares to others available in the literature (references (15-18)) and its specific contribution. The authors mentioned that several agent-based models have been proposed to address the question of interest, but they failed to state what is the gap that their model is filling and in general how it differentiates from previous work. Given the comment on line 74, we can assume the novelty of the developed model lies in not including any "externally imposed behaviors", but this needs to be specified in the introduction.

– The importance of negative feedback loops in the signaling pathways is stated by the authors. Nevertheless, a clear diagram of the responsible regulatory networks is lacking. Such a diagram would allow the reader to clearly judge the presence and details of this or any other feedback in the networks in the context of the proposed conclusions. Additionally, this diagram would allow a clearer interpretation of how the multiple signaling pathways connect in the proposed model.

– As a suggestion for improving the readability of the manuscript, the third paragraph of the introduction (lines 46-55) seems disconnected from the rest of the text. Some rephrasing could help with the flux of the introduction.

– In lines 98 and 99, the authors claim "Specifically, in our model, high Wnt, and Notch signalling environments are required to maintain stemness." Is this statement supported experimentally?

– Often, it is not clear which statements are derived from experimental evidence versus model results. For example, in the sentence on lines 138-140, "a decreased number of enterocytes results in reduced production of BMP, which enables progenitor cells to divide and migrate further up the crypt before meeting BMP levels higher than the differentiation threshold".

– Regarding Figure 1B: The mathematical equations are not necessary. They're not in the main text of the paper. What do the green checkmark and red X mean?

– Regarding Figure 1C: Improperly labeled axes. What are numbers 0, 2, 3, 4, and 64 on the axes? What is p = 0.1 and p = 0.9? Some of these numbers are explained in the supplementary, but not in the main text of the paper.

– Regarding Figure 1D: What types of cells are in the light blue and dark green colors? Where did the authors obtain the data to generate this panel?

– Regarding Figure 2A-D: We find Figure 2A-D not particularly informative. We suggest the authors reevaluate its purpose and redesign it accordingly. For instance, if the purpose is to compare the cell cycle behavior between long and short cycles, panels B and D must be truly comparable, and the x-axis for both figures should be in equal proportions respective to the cell cycle lengths. If the purpose is to show the dependency of p27 over the cell cycle length, the values of p27 need to be shown as well, ideally including a diagram highlighting the interaction between p27 and the rest of the system. Additionally, each panel should have its own caption to increase clarity. Alternatively, the authors could plot the comparison between stem cells and transit-amplifying cells in the same panel and with fewer biochemical species.

– Regarding Figure 2E,F: The simulation and the data don't match. The conclusions drawn from these panels are questionable.

– We suggest rearranging the paragraph on lines 175-186: the purpose (why) of BrdU tracking and Ki-67 staining needs to be specified before explaining how it was done. In its current form, the paragraph can be very confusing for readers not familiar with the experimental technique.

– Some visualization aid (potentially in Figure 1A) would be helpful to understand retrograde movement defined as "when its velocity is negative in the z direction".

– Regarding Figure 3,4: It would be helpful to see side-by-side figures of healthy/homeostatic crypt organization vs the results of the ablation of stem cells or of CDK1 inhibition. Consider including adjacent panels or plotting on top of each other.

– Regarding Figure 3B: The red line for enterocytes appears in the top right corner.

– Regarding Figure 3D: Please use different colors to show directionality. Additionally, red and blue are used to indicate types of cells up to this point. And in general, it is really hard to interpret Figure 3D. Separating the bars (homeostasis vs cell ablation) might help.

– Figure 3E is never referenced in the main text, which is in general bad practice. As its goal is to be a comparison point with Figure 3F, we suggest putting both plots in one panel, each with a clear title.

– To be able to appreciate Figure 4A, we need a reference point, i.e. the case without CDK1 inhibition. The authors show this in a video, but we suggest including it as a figure here too.

– Regarding Figure 4C: Labels describing the cell cycle progress (as in Figure 2B,D) would be helpful.

– Regarding Figure 5: Each panel should have its own caption to increase clarity.

– Regarding Figure 5C: Recovery is after 144 hrs (yellow). That is more than two days (main text). Also, data vs experiments in the purple lines looks qualitatively different.

– Figure 5D: After 6, 24, 96, and 144 hrs, the proportion of proliferative cells is roughly the same. There is no recovery. Please explain why.

– Figure 5E-H: Not referenced in the main text of the paper. Probably not necessary to include.

– With exception of Figure 1, the quality of the figures is poor and can be easily improved. For instance, colors for plots Figure 2A-D can be better chosen to facilitate distinction between the lines (ideally considering colorblindness). We appreciate the attempt to keep the color associated with each cell type consistent through the plots, nevertheless in this case that just complicates the readability of the figures: there are too many cell types to keep track of, and some colors are too similar to each other, some colors are hard to see, there are too many lines in individual panels, and some of the designated colors are used in some other contexts. We suggest simplifying Figure 2A-C and similar plots to only show relevant lines, each with its own legend, and always starting with the y-axis in zero (note: the y-axis range in Figure 2A can be misleading, as it doesn't start in zero, but it isn't obvious at a first sight). Colors can then be chosen to help visualization in Figures 2E-F. (If some lines are excluded for readability in the panels, they can also be shown as figure supplements.)

– Clarify: Line 133 – How do you know what an optimal crypt cell composition is?

– Clarify: Line 149 – Why do you use the Csikasz-Nagy model? Clarify the main text and give some background. Additionally, what is the advantage of this model compared to others models that were described before? Could a similar model be used or was another model tested?

– Clarify: Line 176 – BrdU is field-specific knowledge. Please explain what it does to a more general audience.

– Clarify: Line 284 – What is partially published data?

– The full model could not be evaluated, as it is not available yet in BioModels (even if it appears as submitted). The authors should provide access to the reviewers for a full assessment.

– We suggest slowing down the cell movements/dynamics in all the videos. It is easy to follow the general dynamic of the crypt but a bit difficult to follow how single cells behave.

– A reminder in the video of which color represents each cell type in the crypt would help a lot to follow it.

– Clarify: Line 115 in the Supplementary – Why do the normal distributions have those means and variances? They're very specific.

– Clarify: Please clarify the model organism throughout the paper. Both mouse and human are referenced.

– Discussion: The potential applications of the model should be expanded in the discussion and/or response to reviewers. How exactly can this model be used in different intestinal diseases such as chronic inflammation?

– Discussion: In the last years 3D organoid technology has been a major advance in studying intestinal biology. Would it be possible to apply the model to study organoids' responses in vitro? And study the single-cell response using organoids?

– Discussion: It is well known that small intestinal crypts have a different spatial organization than a colonic crypt, and furthermore, in the small intestinal tissue there are differences among duodenum, ileum, and jejunum. Can this be predicted with the model?

– Discussion: Can the authors expand on how different drugs' actions can be predicted using the same model (e.g. 5 Fu compared to oxaliplatin)?

– Discussion: How close is the model to predicting how many cycles of chemotherapy should be given to a patient to eliminate cancer cells and not cause severe toxicity?

– Discussion: Could a knock-out mouse model of a specific pathway be used to validate the model?

https://doi.org/10.7554/eLife.85478.sa1

Author response

Recommendations for the authors:

– We recommend revising the writing in the introduction and the abstract to clarify the background, motivation, impact, and context of the research. We suggest that the authors explicitly state why studying the intestinal crypt is broadly important to the scientific community and why a mathematical model of it is necessary for progress in the development of oncotherapies. Further, we recommend highlighting the impact of their model of the intestinal crypt on the development of oncotherapies. Most of this information is included in the Discussion section, but not in the introduction or the abstract.

We have revised the abstract which now contains the following paragraph:

“The maintenance of the functional integrity of the intestinal epithelium requires a tight coordination between cell production, migration and shedding along the crypt-villus axis. Dysregulation of these processes may result in loss of the intestinal barrier and disease. With the aim of generating a more complete and integrated understanding of how the epithelium maintains homeostasis and recovers after injury, we have built a multi-scale agent-based model (ABM) of the mouse intestinal epithelium. We demonstrate that a stable self-organizing behaviour in the crypt emerges from the dynamic interaction of multiple signalling pathways, such as Wnt, Notch, BMP, RNF43/ZNRF3 and YAP-Hippo pathways, which regulate proliferation and differentiation, respond to environmental mechanical cues, form feedback mechanisms and modulate the dynamics of the cell cycle protein network.

The model recapitulates the crypt phenotype reported after persistent stem cell ablation and after the inhibition of the CDK1 cycle protein. Moreover, we simulated 5-fluorouracil (5-FU)-induced toxicity at multiple scales starting from DNA and RNA damage, which disrupts the cell cycle, cell signalling, proliferation, differentiation and migration and leads to loss of barrier integrity. During recovery, our in-silico crypt regenerates its structure in a self-organizing, dynamic fashion driven by dedifferentiation and enhanced by negative feedback loops. Thus, the model enables the simulation of xenobiotic-, in particular chemotherapy-, induced mechanisms of intestinal toxicity and epithelial recovery.

Overall, we present a systems model able to simulate the disruption of molecular events and its impact across multiple levels of epithelial organization and demonstrate its application to epithelial research and drug development.”

And we have rewritten the introduction which now includes the following paragraph

“The imbalance of this tightly orchestrated system contributes to pathological conditions, including microbial infections, intestinal inflammatory disorders, extra-intestinal autoimmune diseases, and metabolic disorders (Chelakkot et al, 2018). In addition, critically ill patients and patients receiving chemotherapy/radiotherapy often show severely compromised intestinal barrier integrity (Chelakkot et al., 2018). For instance, oncotherapeutics-induced gastrointestinal toxicity is frequently a life-threatening condition that leads to dose reduction, delay and cessation of treatment and presents a constant challenge for the development of efficient and tolerable cancer treatments (McQuade et al, 2016; Saltz et al, 2000; Saltz et al, 2001; Stein et al, 2010). This intestinal toxicity often results from the interaction of the drug with its intended molecular target such as cell cycle proteins (Zhang et al, 2021) or the disruption of the cycle through DNA-damage (Helleday et al, 2008). Multiscale models integrating our knowledge on how the epithelium maintains homeostasis and responds to injury can contribute to understand epithelial biology and to quantify the risk of intestinal toxicity during drug development.”

– ABM must be defined in the main text, not just the abstract.

We have amended the manuscript as suggested.

– The authors claimed to demonstrate the model's application to epithelial research and drug discovery. They evaluate some known drugs, but we don't believe the term "discovery" is fully justified in the manuscript.

All references to ‘drug discovery’ have been changed to ‘drug development’.

– We suggest that the authors specify how the model they developed compares to others available in the literature (references (15-18)) and its specific contribution. The authors mentioned that several agent-based models have been proposed to address the question of interest, but they failed to state what is the gap that their model is filling and in general how it differentiates from previous work. Given the comment on line 74, we can assume the novelty of the developed model lies in not including any "externally imposed behaviors", but this needs to be specified in the introduction.

A paragraph describing previous models has been added to the introduction and it is attached below. The main novel features of our model are described in detail in the discussion:

“Several agent-based models (ABMs) have been proposed to describe the complexity and dynamic nature of the intestinal crypt. Early models were used as in silico platforms to study the dynamics and cellular organisation of the crypt. For instance, one of the pioneering ABMs was used to study the distribution and organisation of labelling and mitotic indices (Meineke et al, 2001). This model comprises a fixed ring of Paneth cells beneath a row of stem cells, which divide asymmetrically to produce a stem cell and a transit-amplifying cell that terminally differentiates after a fixed number of divisions. Some subsequent models are lattice-free, recapitulate neutral drift of equipotent stem cells and describe proliferation and cell fate regulated by a fixed Wnt signalling spatial gradient, which is defined by the distance from the crypt base, with proliferating cells progressing through discrete phases of the cell cycle and showing variable duration of the G1 phase (Pitt-Francis et al, 2009). Further model refinements can be seen in the model of Buske et al. (2011), with stochastic cell growth and division time (Buske et al, 2011), Wnt levels defined by the fixed local curvature of the crypt and lateral inhibition driven by Notch signalling. Here, we present a lattice-free agent-based model that describes the spatiotemporal dynamics of single cells in the small intestinal crypt driven by the interaction of surface tethered Wnt signals, cell-cell Notch signalling, BMP diffusive signals, RNF43/ZNRF3-mediated feedback mechanisms and the cycle protein network responding to the crypt mechanical environment. We show that our computational model enables the simulation of the ablation and recovery of the stem cell niche as well as of how drug-induced molecular perturbations trigger a cascade of disruptive events spanning from the cell cycle to single cell arrest and/or apoptosis, altered cell migration and turnover and ultimately loss of epithelial integrity.”

Also we have rephrased the sentence “not including externally imposed behaviours”, which now reads as:

“with behaviours (proliferation, differentiation, fate decision, migration, etc.) determined largely by endogenous intracellular and intercellular interactions.”

– The importance of negative feedback loops in the signaling pathways is stated by the authors. Nevertheless, a clear diagram of the responsible regulatory networks is lacking. Such a diagram would allow the reader to clearly judge the presence and details of this or any other feedback in the networks in the context of the proposed conclusions. Additionally, this diagram would allow a clearer interpretation of how the multiple signaling pathways connect in the proposed model.

We have integrated a diagram showing the feedback regulatory loops into Figure 1A.

– As a suggestion for improving the readability of the manuscript, the third paragraph of the introduction (lines 46-55) seems disconnected from the rest of the text. Some rephrasing could help with the flux of the introduction.

We have re-worded this paragraph for a more logical flow, primarily by putting gut-relevant information first to establish why we focus on oncotherapeutics. The paragraph now reads:

“The imbalance of this tightly orchestrated system contributes to pathological conditions, including microbial infections, intestinal inflammatory disorders, extra-intestinal autoimmune diseases, and metabolic disorders (Chelakkot et al., 2018). In addition, critically ill patients and patients receiving chemotherapy/radiotherapy often show severely compromised intestinal barrier integrity (Chelakkot et al., 2018). For instance, oncotherapeutics-induced gastrointestinal toxicity is frequently a life-threatening condition that leads to dose reduction, delay and cessation of treatment and presents a constant challenge for the development of efficient and tolerable cancer treatments (McQuade et al., 2016; Saltz et al., 2000; Saltz et al., 2001; Stein et al., 2010). This intestinal toxicity often results from the interaction of the drug with its intended molecular target such as cell cycle proteins (Zhang et al., 2021) or the disruption of the cycle through DNA-damage (Helleday et al., 2008). Multiscale models integrating our knowledge on how the epithelium maintains homeostasis and responds to injury can contribute to understand epithelial biology and to quantify the risk of intestinal toxicity during drug development.”

– In lines 98 and 99, the authors claim "Specifically, in our model, high Wnt, and Notch signalling environments are required to maintain stemness." Is this statement supported experimentally?

A supporting reference with experimental evidence has been added to this sentence.

Tian et al. Nature 2011 Vol. 478 Issue 7368 Pages 255-9.

– Often, it is not clear which statements are derived from experimental evidence versus model results. For example, in the sentence on lines 138-140, "a decreased number of enterocytes results in reduced production of BMP, which enables progenitor cells to divide and migrate further up the crypt before meeting BMP levels higher than the differentiation threshold".

We have added “in our model” here and throughout the manuscript to clarify modelling assumptions/hypotheses vs experimental data when needed.

– Regarding Figure 1B: The mathematical equations are not necessary. They're not in the main text of the paper. What do the green checkmark and red X mean?

We have revised Figure 1B and equations have been removed. The ‘Green checkmark and red X’ marks are omitted from the new diagram explaining Notch inhibition. Also new diagrams explain the Wnt, BMP and ZNRF3 feedback loops.

– Regarding Figure 1C: Improperly labeled axes. What are numbers 0, 2, 3, 4, and 64 on the axes? What is p = 0.1 and p = 0.9? Some of these numbers are explained in the supplementary, but not in the main text of the paper.

We have removed all numbers from the Notch and Wnt signalling axes and indicated the direction of signalling increase with “Low” and “High”. All parameter values can be found in the Appendix and Table 1. This includes references to probability of fate differentiation, p, and threshold values for Wnt, BMP and Notch axis.

Wnt/BMP axis: replaced with description of values: 64 is replaced with ‘Fate Determinant Wnt Level’ and 4 with ‘BMP and Wnt level for Differentiation’.

Notch axis: removed numerical values, replaced with ‘Low’ and ‘High’ Notch.

The Figure caption has been amended as:

“C) Cell fate determination. High Wnt signalling and activation of Notch are required to maintain stemness. Low Notch signalling determines differentiation into secretory fates, including Paneth cells in high Wnt signalling regions, or goblet/enteroendocrine progenitors in low Wnt regions. Absorptive progenitors develop from stem cells in low Wnt conditions and divide 3-5 times, before becoming terminally differentiated when Wnt signal levels are decreased and cells find sufficient BMP signals;”

– Regarding Figure 1D: What types of cells are in the light blue and dark green colors? Where did the authors obtain the data to generate this panel?

Colours (of enteroendocrine = light blue, secretory progenitor = dark green) in 1A are now amended to be identical to 1D. The same colour code is now used for all plots.

Panel 1D shows the average cell composition, per cell position, of a simulated crypt in homeostasis. The data is generated by the model and shows a good agreement with reported experimental data in (Buske et al., 2011).

We have clarified the Figure caption that now reads:

D) Average composition of a simulated healthy/homeostatic crypt (over 100 simulated days), showing the relative proportion of cells at each position.

and the manuscript text now reads:

“… our model describes single cells that generate and respond to signals and mechanical pressures in the crypt-villus geometry to give rise to a self-organizing crypt which has stable spatial cell composition over time (Figure 1D) and reproduces reported experimental data (Buske et al., 2011).”

– Regarding Figure 2A-D: We find Figure 2A-D not particularly informative. We suggest the authors reevaluate its purpose and redesign it accordingly. For instance, if the purpose is to compare the cell cycle behavior between long and short cycles, panels B and D must be truly comparable, and the x-axis for both figures should be in equal proportions respective to the cell cycle lengths. If the purpose is to show the dependency of p27 over the cell cycle length, the values of p27 need to be shown as well, ideally including a diagram highlighting the interaction between p27 and the rest of the system. Additionally, each panel should have its own caption to increase clarity. Alternatively, the authors could plot the comparison between stem cells and transit-amplifying cells in the same panel and with fewer biochemical species.

We reduced 4 plots, A-D, to 2 plots (A and B). The new plots (A and B) show the dynamics of the main cycle proteins, including p27, in cells with long (figure A) and short (figure B) cycle in a comparable manner, with equal x-axis. We highlighted in the figure caption that a diagram of the protein interaction network, including the interaction of p27 with the rest of the proteins, can be found in the original paper of Csikasz-Nagy, where this model was first described. We did not bring here the diagram due to its size and complexity.

– Regarding Figure 2E,F: The simulation and the data don't match. The conclusions drawn from these panels are questionable.

Figure 2E-F are now Figures 2C-D. As per response to one of the main comments, we have performed several adjustments that include, refinement of the counting algorithm and of the Ki67 and BrdU staining simulations and addition of an estimation of the experimental error to the simulated responses. A description of these changes is described in a new section in the appendix called “ABM simulation of Ki-67 and BrdU Staining”

With these changes we think we have achieved a more satisfactory agreement between observed and predicted results and updated all figures with Ki67 and BrdU staining simulated results.

However, despite our efforts, we did not reach a perfect agreement between BrdU staining simulated and observed results. The simulation behind our predictions is a complex process describing multiple dynamic processes interacting at different scales. Some of our model assumptions, i.e duration of S-phase, dynamics of BrdU integration into the DNA…, may be too simplistic for these datasets and/or our model maybe be missing some important piece of information on BrdU staining process which is not currently understood. This is difficult to identify at the moment, but we are committed to keep refining our models as new knowledge emerges to enhance their performance and application. On the other hand, BrdU staining involves a considerably large experimental error inherent to both large biological and technical variability (sampling region and sampling time, counting, reagents…). Other authors have also encountered difficulties to obtain a good match between simulated and observed BrdU staining results (see Buske et al, 2011).

– We suggest rearranging the paragraph on lines 175-186: the purpose (why) of BrdU tracking and Ki-67 staining needs to be specified before explaining how it was done. In its current form, the paragraph can be very confusing for readers not familiar with the experimental technique.

We have added a new section in the appendix called “ABM simulation of Ki-67 and BrdU Staining” with further details of these simulations and extended the description of these experiments in the following paragraph of the main manuscript.

“To demonstrate the performance of the model to reproduce the spatiotemporal cell dynamics and composition of a homeostatic crypt, we simulated previous published mouse experiments (Parker et al, 2017; Parker et al, 2019) comprising 5-bromo-29-deoxyuridine (BrdU) tracking (Figure 2E) and Ki-67 staining (Figure 2F). BrdU, is a thymidine analogue often used to track proliferative cells and their descendants along the crypt–villus axis (Gratzner, 1982; Nowakowski et al, 1989). BrdU is incorporated into newly synthesized DNA of dividing cells during the S-phase and transmitted to daughter cells, regardless of whether they proliferate. If the exogenous administration of this molecule is discontinued, the cell label content is diluted by each cell division and is no longer detected after 4–5 generations (Wilson et al, 2008). To simulate the BrdU chase experiment after a single BrdU pulse, we assumed that any cell in S-phase incorporated BrdU permanently into its DNA for the first 120 minutes after injection of BrdU and BrdU cell content was diluted upon cell division such that after five cell divisions, BrdU was not detectable. See Appendix for a complete description. The BrdU chase simulation showed that the observed initial distribution of cells in S-phase as well as division, differentiation and migration of BrdU-positive cells over time were replicated by our model (Figure 2E).

Ki-67 is produced by actively proliferating cells during the S-, G2- and M-phases of the division cycle (Sobecki et al, 2017). Ki-67 is detected after exiting the cycle (Miller et al, 2018) and during G1 in continuously cycling cells (Sobecki et al., 2017), which is related to the time required for the protein to be catabolized (Miller et al., 2018). Our simulations assumed that Ki-67 is detected in continuously cycling cells, cells re-entering the cycle after arrest except during G1, as well as in differentiated cells that were cycling within the past 6 hours and recently drug-arrested cells. See Appendix for a complete description. Similarly, we observed that the ABM-simulated spatial distribution along the crypt of Ki-67 positive cells recapitulated observations in mouse ileum (Figure 2F).”

– Some visualization aid (potentially in Figure 1A) would be helpful to understand retrograde movement defined as "when its velocity is negative in the z direction".

We have added arrows showing the direction of regular cellular migration and retrograde motion to Figure 1A. We have attached a low resolution version of Figure 1 here with this change highlighted in a red box.

– Regarding Figure 3,4: It would be helpful to see side-by-side figures of healthy/homeostatic crypt organization vs the results of the ablation of stem cells or of CDK1 inhibition. Consider including adjacent panels or plotting on top of each other.

We have included an image of a simulated 3D crypt in homeostasis in Figure 4 (CDK1 inhibition plot), because that figure is less crowded. In addition, in this figure, we have extended the simulation duration to achieve full recovery of the crypt after ablation and highlighted in the caption that homeostatic cell counts for comparison can be seen after recovery.

– Regarding Figure 3B: The red line for enterocytes appears in the top right corner.

We have divided the plots to show with clarity all cell lineages in the crypt and villus during stem cell ablation and recovery.

– Regarding Figure 3D: Please use different colors to show directionality. Additionally, red and blue are used to indicate types of cells up to this point. And in general, it is really hard to interpret Figure 3D. Separating the bars (homeostasis vs cell ablation) might help.

We have amended this figure as suggested. The bars corresponding to homeostasis and ablation scenarios are now plotted side-by-side (separate) and represented by different colours. Dark/light colour shows upward and downward motion. We have also clarified the legend to improve understandability.

– Figure 3E is never referenced in the main text, which is in general bad practice. As its goal is to be a comparison point with Figure 3F, we suggest putting both plots in one panel, each with a clear title.

We have corrected this typo, the original sentence referred to Figure 3D instead of 3E, which is now referenced in the main text. As suggested these plots are now combined into one panel and explanatory text has been added into the plots.

– To be able to appreciate Figure 4A, we need a reference point, i.e. the case without CDK1 inhibition. The authors show this in a video, but we suggest including it as a figure here too.

As suggested, we have included an image of a simulated 3D crypt in homeostasis side-by-side with the injured crypt. We have also extended the duration of the simulation to achieve homeostasis and highlighted in the caption that homeostatic cell counts and cell cycle protein dynamics can be seen after recovery for comparison.

– Regarding Figure 4C: Labels describing the cell cycle progress (as in Figure 2B,D) would be helpful.

The labels are now added and colours have been modified to match those in other cell cycle figures.

– Regarding Figure 5: Each panel should have its own caption to increase clarity.

Captions have been addressed as suggested.

– Regarding Figure 5C: Recovery is after 144 hrs (yellow). That is more than two days (main text). Also, data vs experiments in the purple lines looks qualitatively different.

The total number of crypt cells is recovered by day 2 after 5-FU high dose treatment is interrupted, however the crypt cell composition is still affected with increased proliferative cells resulting from the cell proliferation enhancing feedback mechanisms responding to injury. We have clarified this in the following section of the manuscript:

“Figure 5C shows that predicted and observed Ki-67 positive cells declined gradually over time at all positions in the crypt during the 5-FU high dose treatment but the numbers recovered, reaching values above baseline, two days after the interruption of 5-FU administration. The increased rebound of the proliferative crypt compartment after the treatment was captured in our ABM by the implemented BMP-mediated feedback mechanism from mature enterocytes to proliferative cells (see BMP signalling section in the Appendix).”

We made several adjustments to improve model predictions: we refined the counting algorithm that determines cell position and improved the Ki67 and BrdU staining simulations by altering the simulated staining criteria and adding of an estimation of the experimental error to the simulated responses. A description of these changes is described in a new section in the appendix called “ABM simulation of Ki-67 and BrdU Staining”. With these changes we think we have achieved a more satisfactory agreement between observed and predicted results and updated all figures with Ki67 and BrdU staining simulated results.

– Figure 5D: After 6, 24, 96, and 144 hrs, the proportion of proliferative cells is roughly the same. There is no recovery. Please explain why.

Figure 5D describes the proportion of BrdU+ cells during 5-FU low dose treatment. This treatment had a a minimal impact on crypt proliferation and on the number of cells on the epithelium. The plots have been moved to the appendix.

– Figure 5E-H: Not referenced in the main text of the paper. Probably not necessary to include.

These plots including the response to 5-FU low dose have now been moved to the Appendix while new plots describing the response of the cell cycle protein network to the drug have been added to Figure 5 (now Figures 5A and B).

– With exception of Figure 1, the quality of the figures is poor and can be easily improved. For instance, colors for plots Figure 2A-D can be better chosen to facilitate distinction between the lines (ideally considering colorblindness). We appreciate the attempt to keep the color associated with each cell type consistent through the plots, nevertheless in this case that just complicates the readability of the figures: there are too many cell types to keep track of, and some colors are too similar to each other, some colors are hard to see, there are too many lines in individual panels, and some of the designated colors are used in some other contexts. We suggest simplifying Figure 2A-C and similar plots to only show relevant lines, each with its own legend, and always starting with the y-axis in zero (note: the y-axis range in Figure 2A can be misleading, as it doesn't start in zero, but it isn't obvious at a first sight). Colors can then be chosen to help visualization in Figures 2E-F. (If some lines are excluded for readability in the panels, they can also be shown as figure supplements.)

We have simplified the plots in the manuscript main figures to increase clarity and improve line visualization. We have also added new Supplementary Figures S1-S2 with all cell lineages in logarithmic scale. Line colours are changed according to the Seaborn’s colourblind palette.

– Clarify: Line 133 – How do you know what an optimal crypt cell composition is?

All references to optimal cell composition have been changes to ‘homeostatic cell composition’.

– Clarify: Line 149 – Why do you use the Csikasz-Nagy model? Clarify the main text and give some background. Additionally, what is the advantage of this model compared to others models that were described before? Could a similar model be used or was another model tested?

To build the cell cycle in our ABM an exhaustive search of cell cycle models was performed with the following criteria: 1) include the core cell cycle proteins (cyclins, CDKs, Wee1, p21/p27), 2) comprise a variable, such as the mass variable of the Csikasz-Nagy model, that we could link to cell features in our ABM, such as the cell size and 3) include sufficient mechanistic detail to allow the description of drug-cell cycle interactions. The Csikasz-Nagy model fulfils all these criteria and was available in the BioModels repository, with accessible, correct code. It is also a well-established model that builds upon the seminal work of Novak and Tyson on the dynamics of cell cycle protein network. In addition, through experimentation, we found this model particularly amenable and robust to constant, dynamic alterations of the parameters required to dynamically vary the length of the cell cycle according to contact inhibition in our ABM.

Other versions of this model could also be integrated in our ABM and a more sophisticated version could be implemented in the future if required to simulate other mechanisms of cycle protein network disruption. Other models with other hypotheses governing protein dynamics remain to be proven.

We have clarified this in the following paragraph in the appendix:

“We have used the model of Csikasz-Nagy et al. (Csikasz-Nagy et al, 2006), that recreates the mammalian cell cycle and is available in Biomodels (Le Novère & Csikasz-Nagy, 2006). This model is an extension of the seminal work of Novak and Tyson that help reveal the complex nonlinear dynamics of the cell cycle proteins (Novak et al, 2001; Novak & Tyson, 1993, 2004). The Csikasz-Nagy model provides multiple necessary features such as core cell cycle proteins, a mass variable that can be coupled to the volume of the single cells in our ABM and sufficient mechanistic detail to enable a detailed description of drug-cycle interactions.”

– Clarify: Line 176 – BrdU is field-specific knowledge. Please explain what it does to a more general audience.

The following paragraph detailing the use of BrdU has been extended in the main text.

“BrdU, is a thymidine analogue often used to track proliferative cells and their descendants along the crypt–villus axis (Gratzner, 1982; Nowakowski et al., 1989). BrdU is incorporated into newly synthesized DNA of dividing cells during the S-phase and transmitted to daughter cells, regardless of whether they proliferate. If the exogenous administration of this molecule is discontinued, the cell label content is diluted by each cell division and is no longer detected after 4–5 generations (Wilson et al., 2008).”

– Clarify: Line 284 – What is partially published data?

That term has been removed form the text.’ Further clarification has been added in material and methods, which now reads as follows:

“We used BrdU tracking and Ki-67 immunostaining data from previously published experiments in healthy mice (Parker et al., 2017; Parker et al., 2019) and following 5-FU treatment (Jardi et al, 2022). The samples from this later study (Jardi et al., 2022) were analysed again to count Ki-67 positive cells at each position along the longitudinal crypt axis, for 30-50 individual hemi crypt units per tissue section per mouse as previously described (Williams et al, 2016).”

– The full model could not be evaluated, as it is not available yet in BioModels (even if it appears as submitted). The authors should provide access to the reviewers for a full assessment.

Model files are provided with this reply. These are ‘.jl’ files for use with Julia. The model (the files provided with this reply) will be freely publicly available through BioModels upon acceptance of this manuscript for publication.

– We suggest slowing down the cell movements/dynamics in all the videos. It is easy to follow the general dynamic of the crypt but a bit difficult to follow how single cells behave.

All videos have been recreated with slower cell movement and a higher frame rate to allow easier tracking of cell behaviour.

– A reminder in the video of which color represents each cell type in the crypt would help a lot to follow it.

In the new version of the videos we have added a fixed legend.

– Clarify: Line 115 in the Supplementary – Why do the normal distributions have those means and variances? They're very specific.

We have expanded that section to justify the choice of parameters for the normal distributions as follows:

“When a proliferative cell is created, it is assigned a desired final size, rfinal=2 3r, where rN(0.35,0.00875) for stem cells and rN(0.5,0.0125) for all other cells. The mean values, 0.5 and 0.35, of the radius of progenitor and stem cells, respectively, were determined for an average, non-proliferative or proliferative progenitor cell to have, without loss of generality, a diameter of 1 while the diameter of an average stem cell is slightly smaller, 0.7. In this way, the model captures the smaller size described for columnar LGR5+ stem cells (Barker et al, 2008), which additionally helps recapitulate the mechanics and cell composition of the niche. The variance of the radius was determined by our implementation of the cell cycle model in the ABM. In our model, the volume of the cell is equated to the cell’s mass parameter of the Csikasz-Nagy model and, hence, the cell final radius determines the duration of the cell cycle as described above. By simulating the cell cycle model, we observed that large values of the standard deviation resulted in some cells progressing through the cycle too quickly and, therefore, failing to complete the cell cycle correctly. This analysis provided an upper limit to the coefficient of variation (CV) = 0.025 to ensure all cells progress regularly through the cycle during homeostasis. This results in values of the standard deviation of the radius of 0.0125 and 0.00875 for progenitor cells and stem cells, respectively. Of note, a cell radius CV of 0.025 corresponds to a cell volume CV of about 0.075 which is not far from the reported experimental CV for cell volume, about 0.11 (Bell & Anderson, 1967).”

– Clarify: Please clarify the model organism throughout the paper. Both mouse and human are referenced.

We have clarified the model organism throughout the text.

– Discussion: The potential applications of the model should be expanded in the discussion and/or response to reviewers. How exactly can this model be used in different intestinal diseases such as chronic inflammation?

In the text, we discussed three applications of our model for the study of epithelial biology by simulating stem cell ablation and to assess epithelial injury associated with two common mechanisms of action of oncotherapeutics, DNA damage and cell cycle protein network disruption. In principle, the implemented inter- and intra-cellular signalling mechanisms can be modified to simulate the effects of different drugs and diseases. Moreover, the model can be developed further to include additional molecular targets by adding new pathways or expanding molecular networks. In the discussion, we have modified the following section to emphasise this, that reads:

“One of the important applications of our modelling approach lies in the development of safer oncotherapeutics. The model enables the prediction of intestinal injury associated with efficacious dosing schedules in order to minimize toxicity while maintaining efficacy of investigational drugs. We demonstrated the application of our model to predict potential intestinal toxicity phenotypes induced by CDK1 inhibition as well as to describe the disruption of the epithelium at multiple scales triggered by RNA and DNA damage leading to the loss of integrity of the intestinal barrier and diarrhea following 5-FU treatment. The drug-induced perturbation of other cell cycle proteins or signalling pathways, already integrated into the model, is straightforward to simulate with the current version of the model while the resolution of molecular networks can be increased, or new pathways incorporated into the ABM, to describe additional drug mechanisms of action.”

For diseases that involve exclusively the disruption of the epithelium (burns, physical erosion…) our model alone may be sufficient if assuming no extra medical complications. However, for modelling inflammatory diseases, such as IBD, or infections, the role of the mucosal immune system, microbiome and other elements such as the mesenchymal and enteric nervous system needs to be carefully assessed and included in the model when relevant.

– Discussion: In the last years 3D organoid technology has been a major advance in studying intestinal biology. Would it be possible to apply the model to study organoids' responses in vitro? And study the single-cell response using organoids?

Yes organoids and ABMs are mutually supportive of each other. ABMs can be used to test hypotheses behind organoid responses and organoids can provide single cell data suitable to develop ABMs with precision. Moreover, we are actively researching on the application of organoids combined with modelling to predict clinical GI toxicity in patients. If successful, this work will form the basis of a full manuscript of its own.

To capture this point in the manuscript we have added the following section:

“While most of the crypt biology understanding integrated in our model derives from mouse epithelial studies, human-derived intestinal organoids and microphysiological systems, now routinely used in research, can provide highly precise information at the single cell level to inform ABM development. In return, ABMs can help test hypotheses behind organoid responses in health and disease conditions.”

– Discussion: It is well known that small intestinal crypts have a different spatial organization than a colonic crypt, and furthermore, in the small intestinal tissue there are differences among duodenum, ileum, and jejunum. Can this be predicted with the model?

We focussed on the small intestine because of its higher proliferative activity and hence relevance to study oncotherapeutics-induced toxicity but also because of the richness and abundance of datasets.

In the current the model, we assumed that differences between jejunum and ileum of the small intestine, can be implemented using the physical parameters of the crypt (height, circumference) and the length-scale of BMP signals to regulate cell spatial organization. The different parametrization is included in the Parameter Table 1 and was based on data obtained in the different tissues. Likewise, we are actively working on modelling the colonic crypt by modifying the parametrization accordingly to the current understanding of the colonic crypt biology, e.g. Paneth cells replaced by deep crypt secretory cells, larger physical dimensions, slower proliferation and distinctive cell spatial organization resulting from different parametrization of Wnt, BMP and RNF43/ZNRF3 signalling pathways.

– Discussion: Can the authors expand on how different drugs' actions can be predicted using the same model (e.g. 5 Fu compared to oxaliplatin)?

The model is set up to describe perturbation of DNA synthesis and various cell cycle proteins as well as the main signalling pathways. Other molecular targets need to be implemented in the model, (but the model was designed to be modular and easily modifiable).

For 5-FU we have modelled FUTP-induced DNA and RNA damage in proliferative cells leading to cell apoptosis and cycle arrest. In the case of oxaliplatin, the kinetics of platinum-DNA adducts in each proliferative cell leading to cell apoptosis will likely need to be integrated into the model.

We added a section to the manuscript (to address an earlier question) that explains further how to use our model to simulate other drug MoA and we bring here:

“Drug-induced perturbation of other cell cycle proteins or signalling pathways, already integrated into the model, are straightforward to simulate with the current version of the model while the resolution of pathways can be increased, or new molecular networks incorporated into the ABM to describe additional drug mechanisms of action.”

– Discussion: How close is the model to predicting how many cycles of chemotherapy should be given to a patient to eliminate cancer cells and not cause severe toxicity?

Mitigating GI toxicity is one of the main applications to drug development of our model. We are already using this model to predict intestinal injury associated with efficacious dosing schedules in order to minimize toxicity while maintaining efficacy of investigational drugs. Mechanistic models are best placed to generate predictions prior to clinical trials and, hence, maximize patients’ safety.

We have added the following section in the manuscript:

“One of the important applications of our modelling approach lies in the development of safer oncotherapeutics. The model enables the prediction of intestinal injury associated with efficacious dosing schedules in order to minimize toxicity while maintaining efficacy of investigational drugs.”

– Discussion: Could a knock-out mouse model of a specific pathway be used to validate the model?

Yes, in fact knowledge on epithelial biology, such as cell fate and signalling regulation, gathered using knock-out mouse experiments has already been instrumental to develop and validate our model. For instance, the new section on sensitivity analysis in the supplementary material shows the reaction of the model to changes in signalling mechanisms and recreates the response of knockout mouse models of the chosen pathways. Another example is the phenotype associated to ZNRF3 knockout which was revealed in knock out experiment and it is now in the text:

“These premises imply that Paneth cells enhance their own production by generating Wnt signals and inducing prolonged division times, which increases stem and Paneth cell production and could lead to unlimited expansion of the niche, recapitulating the phenotype seen in RNF43/ZNRF3 knockout mice (Koo et al, 2012).”

Knock out mouse models are invaluable tools to study epithelial biology and acquire the deep understanding required to develop ABM models. We will continue to refine our model by adding newly discovered behaviours of the epithelium with this or other research tools.

References

Barker N, van de Wetering M, Clevers H (2008) The intestinal stem cell. Genes Dev 22: 1856-1864

Bell GI, Anderson EC (1967) Cell growth and division. I. A mathematical model with applications to cell volume distributions in mammalian suspension cultures. Biophys J 7: 329-351

Buske P, Galle J, Barker N, Aust G, Clevers H, Loeffler M (2011) A comprehensive model of the spatio-temporal stem cell and tissue organisation in the intestinal crypt. PLoS Comput Biol 7: e1001045

Chelakkot C, Ghim J, Ryu SH (2018) Mechanisms regulating intestinal barrier integrity and its pathological implications. Exp Mol Med 50: 1-9

Csikasz-Nagy A, Battogtokh D, Chen KC, Novak B, Tyson JJ (2006) Analysis of a generic model of eukaryotic cell-cycle regulation. Biophys J 90: 4361-4379

Gratzner HG (1982) Monoclonal antibody to 5-bromo- and 5-iododeoxyuridine: A new reagent for detection of DNA replication. Science 218: 474-475

Helleday T, Petermann E, Lundin C, Hodgson B, Sharma RA (2008) DNA repair pathways as targets for cancer therapy. Nat Rev Cancer 8: 193-204

Jardi F, Kelly C, Teague C, Fowler-Williams H, Rodrigues D, Jo H, ., Ferreria S, Herpers B, Van Heerden M, de Kok TM et al (2022) Mouse organoids as an in vitro tool to study the in vivo intestinal response to cytotoxicants. Archives of Toxicology (Accepted)

Koo B-K, Spit M, Jordens I, Low TY, Stange DE, van de Wetering M, van Es JH, Mohammed S, Heck AJR, Maurice MM et al (2012) Tumour suppressor RNF43 is a stem-cell E3 ligase that induces endocytosis of Wnt receptors. Nature 488: 665-669

Le Novère N, Csikasz-Nagy A, 2006. Cell Cycle Model BioModels.

McQuade RM, Stojanovska V, Donald E, Abalo R, Bornstein JC, Nurgali K (2016) Gastrointestinal dysfunction and enteric neurotoxicity following treatment with anticancer chemotherapeutic agent 5-fluorouracil. Neurogastroenterol Motil 28: 1861-1875

Meineke FA, Potten CS, Loeffler M (2001) Cell migration and organization in the intestinal crypt using a lattice-free model. Cell Prolif 34: 253-266

Miller I, Min M, Yang C, Tian C, Gookin S, Carter D, Spencer SL (2018) Ki67 is a Graded Rather than a Binary Marker of Proliferation versus Quiescence. Cell reports 24: 1105-1112.e1105

Novak B, Pataki Z, Ciliberto A, Tyson JJ (2001) Mathematical model of the cell division cycle of fission yeast. Chaos 11: 277-286

Novak B, Tyson JJ (1993) Numerical analysis of a comprehensive model of M-phase control in Xenopus oocyte extracts and intact embryos. J Cell Sci 106 ( Pt 4): 1153-1168

Novak B, Tyson JJ (2004) A model for restriction point control of the mammalian cell cycle. J Theor Biol 230: 563-579

Nowakowski RS, Lewin SB, Miller MW (1989) Bromodeoxyuridine immunohistochemical determination of the lengths of the cell cycle and the DNA-synthetic phase for an anatomically defined population. J Neurocytol 18: 311-318

Parker A, Maclaren OJ, Fletcher AG, Muraro D, Kreuzaler PA, Byrne HM, Maini PK, Watson AJ, Pin C (2017) Cell proliferation within small intestinal crypts is the principal driving force for cell migration on villi. FASEB J 31: 636-649

Parker A, Vaux L, Patterson AM, Modasia A, Muraro D, Fletcher AG, Byrne HM, Maini PK, Watson AJM, Pin C (2019) Elevated apoptosis impairs epithelial cell turnover and shortens villi in TNF-driven intestinal inflammation. Cell Death Dis 10: 108

Pitt-Francis J, Pathmanathan P, Bernabeu MO, Bordas R, Cooper J, Fletcher AG, Mirams GR, Murray P, Osborne JM, Walter A et al (2009) Chaste: A test-driven approach to software development for biological modelling. Comput Phys Commun 180: 2452-2471

Saltz LB, Cox JV, Blanke C, Rosen LS, Fehrenbacher L, Moore MJ, Maroun JA, Ackland SP, Locker PK, Pirotta N et al (2000) Irinotecan plus fluorouracil and leucovorin for metastatic colorectal cancer. Irinotecan Study Group. N Engl J Med 343: 905-914

Saltz LB, Douillard JY, Pirotta N, Alakl M, Gruia G, Awad L, Elfring GL, Locker PK, Miller LL (2001) Irinotecan plus fluorouracil/leucovorin for metastatic colorectal cancer: a new survival standard. Oncologist 6: 81-91

Sobecki M, Mrouj K, Colinge J, Gerbe F, Jay P, Krasinska L, Dulic V, Fisher D (2017) Cell-Cycle Regulation Accounts for Variability in Ki-67 Expression Levels. Cancer Res 77: 2722-2734

Stein A, Voigt W, Jordan K (2010) Chemotherapy-induced diarrhea: pathophysiology, frequency and guideline-based management. Ther Adv Med Oncol 2: 51-63

Williams JM, Duckworth CA, Vowell K, Burkitt MD, Pritchard DM (2016) Intestinal Preparation Techniques for Histological Analysis in the Mouse. Curr Protoc Mouse Biol 6: 148-168

Wilson A, Laurenti E, Oser G, van der Wath RC, Blanco-Bose W, Jaworski M, Offner S, Dunant CF, Eshkind L, Bockamp E et al (2008) Hematopoietic stem cells reversibly switch from dormancy to self-renewal during homeostasis and repair. Cell 135: 1118-1129

Zhang M, Zhang L, Hei R, Li X, Cai H, Wu X, Zheng Q, Cai C (2021) CDK inhibitors in cancer therapy, an overview of recent development. Am J Cancer Res 11: 1913-1935

https://doi.org/10.7554/eLife.85478.sa2

Article and author information

Author details

  1. Louis Gall

    Clinical Pharmacology and Quantitative Pharmacology, Clinical Pharmacology and Safety Sciences, R&D, AstraZeneca, Cambridge, United Kingdom
    Contribution
    Data curation, Formal analysis, Investigation, Methodology, Writing – original draft, Writing – review and editing
    For correspondence
    louis.gall@astrazeneca.com
    Competing interests
    Employee and shareholder of AstraZeneca Plc
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-1805-2357
  2. Carrie Duckworth

    Institute of Systems, Molecular and Integrative Biology, University of Liverpool, Liverpool, United Kingdom
    Contribution
    Supervision, Writing – review and editing
    Competing interests
    No competing interests declared
  3. Ferran Jardi

    Preclinical Sciences and Translational Safety, Janssen, Beerse, Belgium
    Contribution
    Data curation, Writing – review and editing
    Competing interests
    Employee of Johnson & Johnson
  4. Lieve Lammens

    Preclinical Sciences and Translational Safety, Janssen, Beerse, Belgium
    Contribution
    Data curation
    Competing interests
    Employee and shareholder of Johnson & Johnson
  5. Aimee Parker

    Gut Microbes and Health Programme, Quadram Institute, Norwich, United Kingdom
    Contribution
    Data curation
    Competing interests
    No competing interests declared
  6. Ambra Bianco

    Clinical Pharmacology and Safety Sciences, AstraZeneca, Cambridge, United Kingdom
    Contribution
    Supervision
    Competing interests
    Employee and shareholder of AstraZeneca Plc
  7. Holly Kimko

    Clinical Pharmacology and Quantitative Pharmacology, Clinical Pharmacology and Safety Sciences, R&D, AstraZeneca, Cambridge, United Kingdom
    Contribution
    Supervision
    Competing interests
    Employee and shareholder of AstraZeneca Plc
  8. David Mark Pritchard

    Institute of Systems, Molecular and Integrative Biology, University of Liverpool, Liverpool, United Kingdom
    Contribution
    Supervision, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-7971-3561
  9. Carmen Pin

    Clinical Pharmacology and Quantitative Pharmacology, Clinical Pharmacology and Safety Sciences, R&D, AstraZeneca, Cambridge, United Kingdom
    Contribution
    Conceptualization, Resources, Supervision, Funding acquisition, Investigation, Methodology, Writing – original draft, Writing – review and editing
    For correspondence
    carmen.pin@astrazeneca.com
    Competing interests
    Employee and shareholder of AstraZeneca Plc
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-8734-6167

Funding

European Federation of Pharmaceutical Industries and Associations (Innovative Medicines Initiative 2 No. 116030)

  • Louis Gall
  • Carrie Duckworth
  • Ferran Jardi
  • Lieve Lammens
  • David Mark Pritchard
  • Carmen Pin

Horizon 2020 Framework Programme (Innovative Medicines Initiative 2 No. 116030)

  • Louis Gall
  • Carrie Duckworth
  • Ferran Jardi
  • Lieve Lammens
  • David Mark Pritchard
  • Carmen Pin

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

The authors acknowledge financial support from TransQST consortium. This project has received funding from the Innovative Medicines Initiative 2 Joint Undertaking under grant agreement no. 116030. This Joint Undertaking receives support from the European Union’s Horizon 2020 research and innovation programme and EFPIA.

Ethics

All experiments were performed in an Association for Assessment and Accreditation of Laboratory Animal Care approved rodent facility and in accordance with the applicable animal welfare guidelines and legislation. Experimental procedures were approved by the institutional ethics committee. Ten-week-old male C57/BL6Y mice were obtained from Charles River (France). Mice were housed in polysulfon cages with corncob bedding under standard conditions of room temperature (21°C ± 2), relative humidity (55% ± 15) and a 12-h light cycle. Water and a certified rodent pelleted maintenance diet were supplied ad libitum. Nest material and rodent retreats were provided for environmental enrichment.

Senior Editor

  1. Aleksandra M Walczak, École Normale Supérieure - PSL, France

Reviewing Editor

  1. Mariana Gómez-Schiavon, Universidad Nacional Autónoma de México, Mexico

Reviewer

  1. Mariana Gómez-Schiavon, Universidad Nacional Autónoma de México, Mexico

Version history

  1. Received: December 9, 2022
  2. Preprint posted: December 19, 2022 (view preprint)
  3. Accepted: December 7, 2023
  4. Accepted Manuscript published: December 8, 2023 (version 1)
  5. Version of Record published: January 15, 2024 (version 2)

Copyright

© 2023, Gall et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 675
    Page views
  • 126
    Downloads
  • 1
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Louis Gall
  2. Carrie Duckworth
  3. Ferran Jardi
  4. Lieve Lammens
  5. Aimee Parker
  6. Ambra Bianco
  7. Holly Kimko
  8. David Mark Pritchard
  9. Carmen Pin
(2023)
Homeostasis, injury, and recovery dynamics at multiple scales in a self-organizing mouse intestinal crypt
eLife 12:e85478.
https://doi.org/10.7554/eLife.85478

Share this article

https://doi.org/10.7554/eLife.85478

Further reading

    1. Computational and Systems Biology
    James D Brunner, Nicholas Chia
    Research Article

    The microbial community composition in the human gut has a profound effect on human health. This observation has lead to extensive use of microbiome therapies, including over-the-counter 'probiotic' treatments intended to alter the composition of the microbiome. Despite so much promise and commercial interest, the factors that contribute to the success or failure of microbiome-targeted treatments remain unclear. We investigate the biotic interactions that lead to successful engraftment of a novel bacterial strain introduced to the microbiome as in probiotic treatments. We use pairwise genome-scale metabolic modeling with a generalized resource allocation constraint to build a network of interactions between taxa that appear in an experimental engraftment study. We create induced sub-graphs using the taxa present in individual samples and assess the likelihood of invader engraftment based on network structure. To do so, we use a generalized Lotka-Volterra model, which we show has strong ability to predict if a particular invader or probiotic will successfully engraft into an individual's microbiome. Furthermore, we show that the mechanistic nature of the model is useful for revealing which microbe-microbe interactions potentially drive engraftment.

    1. Computational and Systems Biology
    Tae-Yun Kang, Federico Bocci ... Andre Levchenko
    Research Article

    Angiogenesis is a morphogenic process resulting in the formation of new blood vessels from pre-existing ones, usually in hypoxic micro-environments. The initial steps of angiogenesis depend on robust differentiation of oligopotent endothelial cells into the Tip and Stalk phenotypic cell fates, controlled by NOTCH-dependent cell–cell communication. The dynamics of spatial patterning of this cell fate specification are only partially understood. Here, by combining a controlled experimental angiogenesis model with mathematical and computational analyses, we find that the regular spatial Tip–Stalk cell patterning can undergo an order–disorder transition at a relatively high input level of a pro-angiogenic factor VEGF. The resulting differentiation is robust but temporally unstable for most cells, with only a subset of presumptive Tip cells leading sprout extensions. We further find that sprouts form in a manner maximizing their mutual distance, consistent with a Turing-like model that may depend on local enrichment and depletion of fibronectin. Together, our data suggest that NOTCH signaling mediates a robust way of cell differentiation enabling but not instructing subsequent steps in angiogenic morphogenesis, which may require additional cues and self-organization mechanisms. This analysis can assist in further understanding of cell plasticity underlying angiogenesis and other complex morphogenic processes.