The Role of Fibroblast-Epithelial Cross-Talk on the Distribution of Distinct Fibroblast Phenotypes in the Intestinal Crypt.

Introduction
Colorectal cancer (CRC) caused nearly one million deaths globally in 2020 (Xi and Xu 2021), and mortality is expected to double by 2035, including in high-income countries like the United States, Canada, and Australia (Araghi et al 2019). CRC is the third leading cause of cancer death in adults under 50 in the United States, with incidence rising in those over 30 (Bhandari et al 2017). CRC originates in the intestinal epithelium, which is constantly renewed via cell proliferation, migration, and differentiation (Huels and Sansom 2015). The intestinal epithelium consists of a single monolayer, folded into finger-like structures known as villi, and test-tube-shaped structures known as crypts, with each villus surrounded by multiple crypts. Epithelial stem cells, located in the stem cell niche at the base of the crypt, and transit-amplifying cells located above them, proliferate, generating a mitotic pressure that drives upward movement of epithelial cells along the crypt towards the villi. As these cells migrate towards the villi, they differentiate into specialised cell types, such as nutrient-absorbing enterocytes and mucus-secreting goblet cells (Clevers 2013). Epithelial cells then migrate from the crypt into the villi, where they undergo apoptosis and are shed into the lumen. Under homeostasis, epithelial cell proliferation, differentiation, and cell loss are balanced to maintain the crypt-villi structures. Disruption of this balance, often due to mutations in epithelial stem cells, can lead to excessive proliferation and polyp formation, common features of the early stages of CRC (Huels and Sansom 2015). Understanding the cross-talk mechanisms crucial to crypt homeostasis is essential for understanding how mutations in either epithelial cells or fibroblasts contribute to the emergence of phenotypes characteristic of the early stages of CRC (Nassar and Blanpain 2016).
Epithelial cell fate determination is primarily governed by the balance between two opposing morphogen gradients: WNT and Bone Morphogenetic Protein (BMP). WNT maintains the stem cell niche at the base of the crypt and promotes the proliferation of stem and transit-amplifying cells (Horvay and Abud 2013; Wong et al 2002). Overexpression of the WNT inhibitor Dickkopf-1 has been shown to inhibit epithelial stem cell proliferation and reduce the number of crypts in mouse intestines (Kuhnert et al 2004; Pinto et al 2003). In contrast, BMP promotes epithelial cell differentiation and suppresses proliferation. Multiple studies have shown that BMP is essential for regulating epithelial proliferation (Batts et al 2006; He et al 2004; Haramis et al 2004; Howe et al 2001). WNT and BMP signalling are regulated by distinct phenotypic populations of fibroblasts that surround the crypt. Specifically, fibroblasts at the base of the crypt maintain epithelial stem and proliferative phenotypes by expressing high levels of WNT and BMP inhibitors (BMPi), while fibroblasts at the top of the crypt secrete BMP, promoting epithelial cell differentiation (McCarthy et al 2020). Furthermore, BMPi-expressing smooth muscle cells in the mesenchyme below the crypt base inhibit BMP signalling, thus helping to preserve the stem cell niche (McCarthy et al 2023). The spatial arrangement of these fibroblast populations establishes spatial gradients of WNT, BMP, and BMPi that together regulate epithelial cell fate and are characterised by high concentrations of WNT and BMPi at the base of the crypt, and high concentrations of BMP at the top.
Mutations arising in epithelial cells or fibroblasts can destabilise crypt homeostasis and drive uncontrolled epithelial proliferation, a well-established precursor to polyp formation. Some of the most common mutations affect the Adenomatous Polyposis Coli (APC) gene (Kinzler and Vogelstein 1996). APC is a component of a destruction complex that inhibits WNT target gene transcription. Under normal conditions, WNT binding to epithelial cell receptors prevents the formation of the destruction complex and thereby promotes proliferation. However, mutations in the APC protein also inhibit the formation of the destruction complex, resulting in constitutive activation of WNT signalling and uncontrolled epithelial cell growth (Schneikert and Behrens 2007). Mutations in the BMP pathway also contribute to CRC. For example, SMAD, a downstream target of BMP signalling (Kawai et al 2000 which is frequently mutated in CRC patients, inhibits epithelial cell differentiation (Fleming et al 2013; Koyama et al 1999). Mutations in epithelial cells can alter fibroblast phenotypes. For example, APC mutation-driven changes in WNT concentrations can induce changes in fibroblast phenotypes that inhibit the epithelial-to-mesenchymal transition in CRC tumours (Mosa et al 2020).
It is well established that fibroblasts play a crucial role in regulating epithelial cell behaviour within intestinal crypts, and also influence tumour progression (Alkasalias et al 2018). However, the mechanisms by which fibroblasts acquire distinctive phenotypes along the crypt axis remain unclear. A growing body of evidence supports the view that fibroblasts and epithelial cells have co-evolved a morphogen-mediated cross-talk mechanism that enables mutual regulation of their phenotypes. Studies in chick embryo intestines have shown that fibroblast BMP expression is spatially correlated with epithelial hedgehog (hhog) expression by epithelial cells at the villus tips and the top of the crypt (Shyer et al 2015). Moreover, inhibition of hhog leads to reduced BMP activity and structural defects in the intestinal epithelium. Additional studies have shown that BMP4 expression is regulated by hhog signalling (Ormestad et al 2006; Roberts et al 1995; Walton et al 2012), suggesting that hhog stimulates BMP production by fibroblasts. Since BMP expression is a defining feature of the fibroblast phenotype at the top of the crypt, these findings suggest that hhog plays a central role in regulating subepithelial fibroblasts in this region. Furthermore, Shyer et al (2015) showed that hhog expression itself is modulated by subepithelial fibroblast phenotypes, consistent with the existence of a feedback loop. Together, this evidence supports the idea that epithelial cells and fibroblasts adopt mutually supporting phenotypes through a positive feedback loop mediated by BMP and hhog signalling.
In addition, there is growing evidence that WNT and hhog negatively regulate each other. Studies have shown that blocking WNT increases hhog expression (Degirmenci et al 2018; Ormestad et al 2006; Nik et al 2013). Hhog also indirectly influences BMPi expression. Kraiczy et al (2023) demonstrated that fibroblast phenotypes are influenced by their proximity to BMP-secreting cells. Fibroblasts near these cells adopt a phenotype that supports epithelial cell differentiation by expressing BMP, while fibroblasts located further away express WNT and BMPi, supporting stem-like epithelial phenotypes. Therefore, WNT inhibition of hhog not only reduces BMP expression but also promotes BMPi expression. This relationship supports the idea that fibroblast–epithelial cross-talk, modulated by proximity to BMPi-secreting mesenchymal cells identified by McCarthy et al (2023), is key to the emergence of distinct, mutually supporting fibroblast–epithelial phenotype pairs at the top and base of the crypt. This cross-talk may also explain how mutations in one cell population alter the phenotype of the other. Figure 1 summarises the fibroblast–epithelial interactions that we use to develop our mathematical model of cross-talk in intestinal crypts.
Understanding how mutations in either epithelial cells or fibroblasts influence the phenotype of the opposite cell population is crucial to identify which mutations may contribute to CRC development, and its treatment. Mathematical models offer a valuable tool for understanding how changes in cell phenotypes are regulated by WNT, BMP, and hhog morphogens. Most existing models of epithelial cell fate focus on subcellular signalling pathways, and biochemical reactions involving proteins that regulate gene transcription and, ultimately, cell phenotype and function. Among these models, those describing the WNT pathway have been the most extensively explored. Lee et al (2003) developed the first dynamic model of the WNT pathway, which predicts a cell’s “stemness” based on the concentration of the -catenin/T-cell factor (TCF) complex, a regulator of genes related to cell-cycle progression. Subsequent models have either extended or simplified these equations (Benary et al 2013; Cho et al 2006; Kay et al 2017; Kim et al 2007; MacLean et al 2015; Pedersen et al 2011; Shin et al 2014; Schmitz et al 2011, 2013; van Leeuwen et al 2007; Wawra et al 2007). For example, Kay et al (2017) extended the Lee et al (2003) model to investigate how cross-talk between WNT and Delta-Notch signalling impacts epithelial cell fate specification.
Fewer comparable models exist for the BMP and hhog pathways. Lai et al (2004) proposed the first model of the hhog pathway, Balaskas et al (2012) developed a model of the transcriptional network acting downstream of hhog to study pattern formation in vertebrate neural tube, and subsequent models (Cohen et al 2014, 2015) explored how hhog gradients influence mouse embryo morphogenesis. Several models have examined TGF- regulation of SMAD (Chung et al 2009; Liu et al 2014; Nicklas and Saiz 2013; Tian et al 2013; Warsinske et al 2015; Zhang et al 2014; Zi et al 2012), a key subcellular protein that regulates fibroblast phenotype (Kawai et al 2000). Of these models, only Nicklas and Saiz (2013) accounts for the effect of BMP on subcellular SMAD concentrations.
In this study, we develop a mathematical model of fibroblast–epithelial cross-talk within the intestinal crypt, with the aim of explaining how distinct fibroblast and epithelial cell phenotypes emerge. The model focuses on interactions between the morphogens BMP, WNT, BMPi, and hhog. Morphogen-mediated interactions between epithelial cells and fibroblasts are modelled by modulating extracellular morphogen expression rates as linear functions of the concentrations of bound or unbound morphogen receptors on each cell type. Rather than explicitly representing epithelial cell and fibroblast populations, we capture their influence on the local environment through variables corresponding to morphogen receptor concentrations. Accordingly, we treat epithelial cell and fibroblast phenotypes as existing on a continuous spectrum, determined by the proportion of bound morphogen receptors. For simplicity, we model the concentrations of morphogens and receptors as a well-mixed system of ordinary differential equations (ODEs), based on a set of reactions and the law of mass action. To address the absence of spatial heterogeneity inherent in ODE systems, we introduce a parameter representing the BMPi influx from the BMPi-emitting mesenchyme at the base of the crypt. This BMPi influx acts as a control parameter to implicitly represent spatial variation, with larger values indicating a closer proximity to the base of the crypt.
Changes in fibroblast and epithelial cell phenotypes between the crypt and villus junctions are believed to be abrupt, implying a discrete phenotypic switch rather than a smooth, graded transition (McCarthy et al 2020). One approach to modelling a discrete switch in the epithelial cell phenotype, based on a continuous WNT pathway model, is to identify parameter regimes in which the system exhibits a bistable switch between two steady state solutions. MacLean et al (2015) and Shin et al (2014) demonstrated that WNT pathway models can admit two stable steady state solutions that correspond to proliferative and quiescent phenotypes, and show how epithelial cells may switch between these phenotypes in response to external stimuli. A bistable switch in response to spatial variation along the crypt axis provides a potential mechanistic explanation for the stem-supporting and BMP-emitting fibroblast phenotypes observed by McCarthy et al (2020). We use our model to: Demonstrate how fibroblast–epithelial cross-talk interactions and implicit spatial variation establish at least two pairs of mutually supporting fibroblast and epithelial cell phenotypes;

Identify cross-talk interactions that are necessary for the formation of these phenotypic pairs, and time-dependent switching between the phenotypic pairs;

Investigate how mutations in epithelial cells or fibroblasts dysregulate cross-talk dynamics, driving shifts in the phenotypes of adjacent cells, a common feature of the early stages of CRC.

The remainder of the paper is structured as follows. In Section 2, we introduce our model of fibroblast–epithelial cross-talk. Section 3 presents numerical results from the full model, as well as from a series of submodels in which specific cross-talk interactions are altered to investigate their roles in maintaining crypt homeostasis. In Section 4, we examine the implications of our findings, and in Section 5, we outline potential directions for future research.

Model development
Our mathematical model describes the evolution of four key morphogens: BMP (), hhog (), WNT () and BMPi (); the BMP/BMPi complex (); and the concentrations of bound/unbound fibroblast and epithelial cell morphogen receptors. We assume that WNT receptors are only expressed by epithelial cells, hhog receptors by fibroblasts, while BMP receptors are expressed by both. We denote the concentrations of bound and unbound epithelial cell WNT receptors by and , the concentrations of bound and unbound fibroblast hhog receptors by and , and the concentrations of bound and unbound fibroblast and epithelial cell BMP receptors by (, ). All dependent variables are measured in units of . We use the principle of mass balance and the law of mass action to construct a system of ODEs describing the time evolution of the dependent variables. All model parameters are non-negative and represent reaction rate constants, unless otherwise specified.
We assume that morphogens bind reversibly to their respective receptors, and that the total concentration of receptors on each cell type are maintained at a constant level. Under these assumptions, the evolution of the receptor concentrations can be described by the following ODEs:where and denote the forward and backward binding rates respectively. As a consequence of our assumption of mass balance for the morphogen receptors, the pairs of equations for bound and unbound receptors for each type sum to 0. Therefore, the differential equations for the receptor concentrations can be simplified using the conservation of the total receptor concentrations:where and represent the total concentrations of fibroblast BMP and hhog receptors, respectively, while and represent the total concentrations of epithelial cell BMP and WNT receptors. These parameters represent concentrations with units of M and are not rate constants. Their values correspond to the sum of the initial concentrations of bound and unbound receptors.
In addition to binding to fibroblast and epithelial cell receptors, BMP also binds reversibly with extracellular BMPi, preventing BMP from binding to cell receptors. We also assume that the BMP/BMPi complex decays exponentially, as it resides in the extracellular environment and is not subject to cellular regulation. Combining these assumptions, we deduce that the ODE governing the evolution of the BMP/BMPi complex is as follows:where and are the forward binding rate, backward binding rate and decay rate respectively.
We assume that three processes dominate/regulate the dynamics of the morphogens: exponential decay, reversible binding, and source terms. The source terms represent the rate of morphogen expression by epithelial cells and fibroblasts at rates which depend on the concentration of bound receptors on each cell type. Cross-talk between the cell types is mediated by the morphogens. We classify source terms as either stimulatory or inhibitory interactions, depending on whether a morphogen’s expression is positively or negatively correlated with another morphogen’s expression. Fibroblasts are the primary source of BMP, WNT and BMPi, while epithelial cells are the primary source of hhog. BMP expression is positively correlated with BMP and hhog; hhog is positively correlated with BMP and negatively correlated with WNT; and WNT and BMPi are negatively correlated with BMP and hhog. Figure 2 summarises these eight interactions: blue arrows represent reversible binding of morphogens, black arrows denote stimulatory interactions, and blunted arrows indicate inhibitory interactions. By combining the above assumptions and the interactions in Figure 2, we derive the following ODEs for the time evolution of the four morphogens:where and , as defined by Equations (9)–(12), represent the contributions to morphogen concentrations arising from morphogens binding to, and unbinding from, their target receptors. - are the rate constants for the cross-talk interactions, and - are the decay rates for the morphogens. In Equation (17), represents the influx of BMPi into the system from the mesenchyme beneath the crypt. acts as a surrogate for spatial position, with larger values of corresponding to regions closer to the base of the crypt. Therefore, we use as a control parameter for distance from the mesenchyme layer at the base of the crypt.
Before analysing the behaviour of our model, it is convenient to non-dimensionalise the system of ODEs. Using the variable and parameter reductions detailed in Tables 1 and 2, the dimensionless form of Equations (9)–(17) are given by:
In Section 3, we identify the steady-state solutions of the dimensionless model and show how their existence and linear stability change as varies. The steady-state solutions solve the following system of algebraic equations:where:

Results
The primary objective of this paper is to explore how epithelial-fibroblast cross-talk, coupled with implicit spatial variation, gives rise to the experimentally observed morphogen gradients and multiple pairs of mutually supporting epithelial cell and fibroblast phenotypes. We begin by identifying parameter values for which the model admits at least two physically realistic steady-state solutions, each representing a mutually supporting pair of epithelial-fibroblast phenotypes, using the Advanced Deficiency Algorithm (Ellison 1998) (see Appendix A.1 for details). These parameter values and their corresponding solutions are then used to generate bifurcation diagrams showing how the steady-state solutions change as varies, using numerical continuation (see Appendix A.2). Recall that represents the local influx of BMPi from the mesenchyme at the base of the crypt, and thus acts as a proxy for spatial position along the crypt axis. We then use linear stability analysis to determine the local stability of these steady-states (see Appendix A.3). Finally, we present time-dependent solutions in which is a prescribed function of time. These simulations allow us to examine how epithelial cell phenotypes evolve as epithelial cells migrate up the crypt.
In Section 3.1, we characterise the bifurcation structure of the steady-state solutions to equations (27)–(35) as varies, and identify parameter values that give rise to multiple pairs of mutually supporting fibroblast and epithelial cell phenotypes. Particular attention is given to parameter regimes that produce fibroblast phenotypes consistent with those observed by McCarthy et al (2020). These results describe the behaviour of a healthy crypt under homeostatic conditions. In Section 3.2, we analyse a series of submodels in which one or more of the cross-talk interactions between fibroblasts and epithelial cells have been modified. This analysis reveals interactions which may be essential for homeostasis of epithelial cell and fibroblast phenotypes in the crypt. In Section 3.3, we explore how the bifurcation structure of the model changes as key cross-talk parameters vary, in order to understand how mutations may disrupt homeostasis. Since we study the behaviour of the model for different values of the cross-talk parameters, each figure represents the model behaviour for a different set of model parameters. Therefore, we have included a list of the model parameter values used for each figure is included in Appendix A.

Model Behaviour
The bifurcation diagrams in Figures 3 and 4 show how the steady-state solutions for the proportion of bound receptors and morphogens change as varies. For large values of (), the model admits a unique physically realistic steady-state solution that is linearly stable. This solution is characterised by fibroblasts with low proportions of bound hhog and BMP receptors, and epithelial cells with a low proportion of bound BMP receptors and a high proportion of bound WNT receptors. Furthermore, extracellular concentrations of both BMP and hhog are low, while extracellular concentrations of WNT and BMPi are high. This pattern of extracellular morphogen levels and receptor binding is consistent with the epithelial stem cell and stem cell supporting fibroblast phenotypes found at the base of the crypt.
As decreases the epithelial stem cell and stem supporting fibroblast solution branch persists. Two new solution branches emerge from a fold bifurcation at a critical value of (where ), one stable and one unstable. Steady-state solutions on the new stable branch exhibit a reversed pattern of bound receptor proportions to the stem cell solution branch: the proportions of bound fibroblast BMP and hhog receptors are high, and the proportion of bound epithelial cell BMP receptors is high while the proportion of bound epithelial cell WNT receptors is low. The morphogen concentrations are also reversed, with high concentrations of extracellular BMP and hhog, and lower concentrations of WNT and BMPi. This behaviour is consistent with differentiated epithelial cells and their supporting fibroblasts found near the top of a healthy crypt.
All three solution branches persist as decreases until passes through a critical value (), at which the unstable solution branch and the solution branch that corresponds to the epithelial stem cells and their supporting fibroblasts meet at a fold bifurcation. For smaller values of (), the only stable solution branch that persists corresponds to differentiated epithelial cells and their supporting fibroblasts.
Recall that BMPi is secreted by the cells in the sub-epithelial mesenchyme at the base of the crypt, and that the value of is a surrogate for proximity to the base of the crypt. Therefore, the results in Figures 3 and 4 show how epithelial-fibroblast cross-talk partitions the crypt into three spatially distinct zones, each characterised by a different combination of cell phenotypes:Stem epithelial cell zone: near the base of the crypt, where and the model admits a unique stable steady-state characteristic of epithelial stem cells and stem supporting fibroblasts;

Differentiated epithelial cell zone: near the top of the crypt, where , the model exhibits a unique stable steady-state characteristic of differentiated epithelial cells and their supporting fibroblasts;

Transitional epithelial cell zone: in the middle of the crypt, when , the model admits two stable and one unstable steady-state solutions.

Figure 5 illustrates how the steady-state solutions in Figures 3 and 4 organise the cell phenotypes into three distinct zones. The Differentiated epithelial cell zone contains only differentiated epithelial cells and their mutually supporting fibroblasts; the Stem epithelial cell zone contains only epithelial stem cells and their supporting fibroblasts; and the Transitional epithelial cell zone contains a mixture of both pairs of phenotypes, because both steady-states are viable within this interval. Figure 5 also shows that the morphogen gradients exhibit a discrete step change between the Transitional and Differentiated epithelial cell zones, as predicted by the fold bifurcations of the steady-states. Although diffusion in vivo is likely to smooth this transition, the steady-states of our model suggest that the underlying morphogen gradients are unlikely to evolve gradually at the transitions between phenotypic zones, in contrast to the assumptions made in many existing models.
We now use our model to investigate how epithelial cell phenotypes change as varies according to the function shown in Figure 6a. We choose this function of because this function demonstrates the effect of the value of smoothly transitioning through all three phenotype zones. As decreases from its initial value, the system remains on the solution branch associated with stem epithelial cells until . As decreases beyond this fold point, the system evolves towards the solution branch characterised by differentiated epithelial cells and their supporting fibroblasts. The system remains on this solution branch until increases through , when it returns to the stable solution branch associated with epithelial stem cells. The hysteresis loop demonstrates two key behaviours of the epithelial cell and fibroblast phenotypes. The sharp switch of the time-dependent solution as decreases mimics the process of epithelial cell differentiation as they migrate up the crypt. The second is the robustness of the differentiated epithelial cell phenotypes. The time-dependent solution only switches from the differentiated to the stem branch once the value of belongs to the Stem epithelial cell zone, indicating that naturally occurring fluctuations in the value of are unlikely to trigger epithelial cell de-differentiation.

Identifying Key Cross-Talk Interactions
The steady-state solutions of the full model show how distinct phenotypic zones may form at different positions along the axis of a healthy crypt. We now investigate how the cross-talk interactions included in the model contribute to the formation of the phenotypic zones by considering submodels in which one or more of the cross-talk interactions are altered.
The rates of morphogen expression by fibroblasts and epithelial cells are linear functions of the proportion of bound morphogen receptors, with the -parameters (), listed in Table 3 as the coefficients. These parameters represent the maximal rate of morphogen expression for each interaction in Figure 2. We consider two modifications to these interactions: setting parameters to zero or setting the contribution of the cross-talk interaction to a constant value determined by the value of the parameter. For instance, changing interaction such that results in no BMP-mediated BMP stimulus, even when all fibroblast BMP receptors are bound. Conversely, another possible alteration is to set the rate of interaction to a constant rate of . In effect, setting the rate of BMP-mediate BMP stimulus to the maximal possible rate regardless of the proportion of bound fibroblast BMP receptors.
These modifications have two biological interpretations depending on whether the interaction is a stimulatory or inhibitory interaction. As an example of a stimulus interaction, fixing effectively removes BMP-mediated BMP stimulus. In contrast, setting results in maximal inhibition of hhog expression regardless of the proportion of bound epithelial cell WNT receptors, effectively modelling the scenario in which the WNT pathway is permanently activated. Setting the contribution of interaction to a constant value of , is interpreted as inactivating WNT-mediated inhibition of hhog.

Inactivation of hhog Mediated BMP Stimulus
The No BMP stimulus model is derived by fixing . So that hhog-mediated BMP stimulus is inactive. Figure 7 shows how the proportions of bound BMP and WNT receptors on the epithelial cells vary with when . As for the full model (see Section 3.1), when is large, the submodel admits a single steady-state solution characteristic of epithelial stem cells and their supporting fibroblasts: the proportion of bound epithelial cell WNT receptors is high and the proportion of bound epithelial cell BMP receptors is low, while the proportions of bound BMP and hhog fibroblast receptors are low (results not shown). As for the full model, a pair of steady-state solutions emerge at a fold bifurcation at a critical value of : one unstable, and one stable steady-state solution characteristic of differentiated epithelial cells and their supporting fibroblasts. In contrast to the full model, the unstable solution branch and epithelial stem cell branch do not meet at a fold bifurcation for some value of . Instead, the unstable solution branch intersects with the epithelial stem cell branch at a transcritical bifurcation. The existence of a transcritical bifurcation follows naturally from the steady-state equations when . The steady-state equations for the No BMP stimulus model are the same as the full model (Equations (27)-(35)), except that Equation (32) is replaced by:Substituting the expressions for and from Equations (27) and (31) into Equation (36), we derive the following equation:Therefore, when one steady-state solution exists for all values of , such that:Substituting and from Equations (28) and (30) with and into Equations (33) and (34) yields the following equations for and :Elimination of between Equations (39) and (40) yields the following quadratic for the corresponding steady-state solutions for :It is straightforward to show that Equation (41) only admits one positive root, . Negative variable values are considered physically unrealistic; therefore, we only consider the positive solution of Equation (41). The value of the other variables can be defined in terms of as follows:Since this solution branch exists for all values of , it cannot terminate at a fold bifurcation. is a repeated route of Equation (37) when the values of and are:Therefore, the two solutions branches must intersect when . The Jacobian of the equations at the intersection point admits a zero eigenvalue, and corresponds in a change in the stability of the epithelial stem cell solution branch (See Appendix C.1). Thus the intersection of the two curves is a transcritical bifurcation, and exists for any set of physically positive parameter values, but may not exist for and .
Motivated by our analysis of the hysteresis loop exhibited by the full model (see Figure 6), we now examine phenotype switching in the No BMP stimulus model as decreases monotonically from to 0. We omit the investigation of time-dependent behaviour as increases because the boundary between the Stem and Transitional epithelial cell zones is formed by a fold bifurcation regardless of the value of . As a result, the qualitative behaviour of the switch between the differentiated epithelial cell solution branch and the epithelial stem cell solution branch is unchanged. To quantify the effect of hhog-mediated BMP stimulus on switching dynamics, we compare time-dependent solutions of the model with against the time-dependent solution when . It is possible to show that a small increase in the value of always exists for which the transcritical bifurcation splits into two fold bifurcations, provided the system admits bistability when (see Appendix C.1 for details). Figures 8a and 8bshows the bifurcation diagrams of the model when and . Figure 8dshows the corresponding time-dependent solutions for different values of , for the same prescribed function of shown in Figure 8c. These results demonstrate that the No BMP stimulus model does admit phenotype switching, but the switching is slow relative to the full model with a positive rate of hhog-mediate BMP stimulus. One interpretation of the slower switching time is that epithelial cells migrating up the crypt retain their stem phenotype until they are eventually shed into the lumen. This would imply that differentiated epithelial cells may not be present in the crypt, therefore, we conclude that cross-talk is necessary crypt homeostasis.

Maximal WNT Inhibition
The Maximal WNT inhibition model, is obtained by setting , and results in complete WNT-mediated hhog inhibition regardless of the proportion of bound epithelial cell WNT receptors. When Equation (33) reduces to give:As in Section 3.2.1, it is possible to show that, at steady-state, and solves the following equation:We note that the steady-state with exists for all positive parameter values, and the other model variables values are given by:This solution branch is characteristic of epithelial stem cells and their supporting fibroblasts. As with the No BMP stimulus model, the unstable solution branch and the epithelial stem cell branch intersect at a transcritical bifurcation when:(see Appendix C.2 for the details)
The time-dependent behaviour of the model when is similar to the behaviour of the model when . Similarly a small increase in the value of , causes the transcritical bifurcation to split into two fold bifurcations recovering faster switching time between the solution branches (See Appendix C.2 for details).

Inactive WNT Receptors
Next we consider the Inactive WNT receptor model, describes a scenario where WNT no longer inhibits hhog expression by epithelial cells. In this model, the source term corresponding to interaction in Figure 2 is now a constant value of . This effectively mimics the condition by which the epithelial cell WNT pathway is completely inactive. As a result, downstream processes typically triggered in WNT binding, including the inhibition of hhog, are entirely absent. Detailed modifications to the equations accompanying this submodel are described in Appendix C.3.
Similarly to the full model discussed in Section 3.1, the bifurcation structure of the model partitions the crypt into three distinct zones separated by two fold bifurcations (Figure 18). The presence of the two fold bifurcations demonstrate that the inhibition of hhog by WNT is not a critical requirement for the formation of different phenotypic zones within the crypt.

Minimal Viable Model
Comparison of Figures 3, 7, and 9 reveal that two different bifurcation structures can generate the three phenotypic zones. However, when either or phenotype switching is much slower than when . We conclude that and are required for the bifurcation structure that facilitates the existence of three phenotypic zones and rapid switching between epithelial cell phenotypes for the time-dependent solutions.
These results pose a natural question with respect to the minimum set of parameters necessary to maintain the three phenotypic zones and fast phenotype switching. It is possible to infer from the steady-state equations that if all parameters other than and are equal to 0, the system only admits one solution with non-negative values for all variables, indicating that additional parameters must be nonzero for the model to admit multiple steady-state solutions. To understand how many additional alpha parameters are required to be positive, we constructed the Minimum viable model by setting all , and to zero. The only non-zero parameters and , are the rate constants for interactions , , and in Figure 2. The minimal viable model also admits the three phenotypic zones, separated by the value of at two fold bifurcations (Figure 19). This demonstrates that only one additional cross-talk dynamic is necessary to recover the key behaviour of the full crypt model.
The inclusion of any single additional parameter alongside and does not necessarily reproduce the behaviour observed in the full crypt model. For instance, a submodel for which only , , and are positive retains the three phenotypic zones separated by two fold bifurcations. However, the WNT concentration remains zero across all values of , rendering the model phenomenologically unrealistic. Table 3 highlights that interactions , , and possess a unique structure not replicated by any other choice of three interactions. These three interactions correspond to one interaction for BMP, hhog, and WNT as the bound and product morphogens, and must include interactions and . Furthermore, these three interactions form a closed cycle of regulatory interactions: for interaction , hhog is the bound morphogen and BMP is the product; for interaction , BMP is the bound morphogen and WNT is the product; and for interaction , WNT is the bound morphogen and hhog is the product. No other combination of three interactions completes such a cycle with BMP, hhog and WNT.

Variation in Cross-Talk Dynamics
In Section 3.2.4, we showed that the existence of the three phenotypic zones and fast phenotype switching are robust to removing most cross-talk interactions. While simpler submodels share the same characteristics as the full model in Section 3.1, the parameter values for the simpler models are significantly different from the parameter values for the full model. This suggests that the values of and must change to compensate for the loss of the other cross-talk interactions. Although the submodels are useful for determining which interactions are important for maintaining the bifurcation structure and time-dependent behaviour of the full crypt model, they do not realistically capture the gradual impact of mutations in either fibroblasts or epithelial cells on cross-talk. Initially, mutations typically arise in a small subset of cells. Consequently, the crypt behaviour remains initially unchanged and evolves progressively as mutated cells accumulate. Therefore, gradual changes in the bifurcation structure as the parameters continuously vary are a more realistic framework to understand how mutations change the organisation of the fibroblast and epithelial cell phenotypes. In this section, we investigate how the bifurcation diagrams of the full model as varies change for different values of . We will focus on how the three phenotypic zones that characterise the the model change as varies for different values of .
In Figure 10 we plot the steady-state solutions of the full model for , and 40 as varies. When , the model admits three distinct phenotypic zones (see Figure 11a). Two types of disruption to this zonal structure occur as varies. First, when , the two fold bifurcations disappear, and all solutions belong to a unique stable steady-state solution branch for all values of . For small values of , the steady-state solutions are characterised by small proportions of bound WNT receptors and large proportions of bound BMP and hhog receptors, resembling the steady-state solution associated with the Differentiated epithelial cell zone at the top of the crypt. For larger values of , the receptor-binding profile resembles the steady-state solution associated with the Stem epithelial cell zone. However, unlike the behaviour of the model when , the steady-state solutions have a graded transition between these extremal solutions. Although the transition remains sharp, it becomes more gradual with increasing (results not shown), indicating that as increases intermediate phenotypes occur for a larger range of values of .
The second type of disruption to the phenotypic zones are shifts in the positions of the fold bifurcations, specifically, the BMP high/WNT low and BMP low/WNT high folds at and , respectively. Since represents the influx rate of BMPi from the mesenchyme at the base of the crypt, only non-negative values of are physically realistic. The values of and define the boundaries between the Stem, Transitional, and Differentiated epithelial cell zones. If either bifurcation occurs at , then the Differentiated or the Differentiated and Transitional epithelial cell zones ceases to exist within the physiologically relevant range. When , , indicating a loss of the Differentiated epithelial cell zone (Figure 11b). When , , and the Differentiated and Transitional zones disappear (Figure 11c). In both cases, the model predicts that the epithelial stem cells fail to differentiate as decreases (i.e as they migrate up the crypt). This is because the time-dependent solution no longer crosses a fold bifurcation as decreases, and therefore, does not trigger phenotype switching. and are also both less than 0 when , but as described in Section 3.2.1, the bifurcation forming the boundary between the Differentiated and Transitional zones is a transcritical bifurcation rather than a fold bifurcation.
Figure 12 shows the plots of the values of and as varies, and the corresponding changes to the phenotypic zones. The vertical coloured lines correspond to the solution plots in Figure 10. As increases, the two fold bifurcations emerge at values of and respectively. Both fold bifurcations exist for larger values of , but for negative values of which are biologically unrealistic. Therefore, and represent the loss of the Differentiated and Transitional epithelial cell zones respectively. The values of and both increase until they meet at a cusp point at . As a consequence, the Transitional epithelial and the start of the Differentiated epithelial cell zone migrate down the crypt as increases, until two fold bifurcations meet at a cusp, creating a region of monostability.
We performed a similar analysis for all other parameters. For each -parameter, except , varying the -parameters continuously resulted in the loss of the phenotypic zones and the merging of the fold bifurcations at a cusp observed in Figure 12. For some of the parameters, only the Differentiated epithelial cell zone was lost. However, the loss of just the Differentiated epithelial cell zone is sufficient for the disruption of the behaviour of a healthy crypt because the time-dependent solution does not switch between the epithelial stem cell solution branch to the differentiated epithelial solution branch without the existence of Differentiated epithelial cell zone. The results of varying the different parameters are provided in Appendix D.
In addition to studying the disruption of the phenotypic zones by varying each of the parameters individually, we also examined whether phenotypic zones lost due to changes in one parameter could be recovered by adjusting other parameters. We show that this is, in principle, possible by demonstrating in Appendix D that the loss of the phenotypic zones for small values of can be reversed by changing the value of .
The two types of disruption to the phenotypic zones also occur when the bifurcation marking the boundary between the Differentiated and Transitional epithelial cell zones is a transcritical bifurcation instead of a fold bifurcation. For example, in Appendix D, we show the loss of the Differentiated and Transitional epithelial cell zones for the No BMP stimulus model as increases, and the formation of a monostable zone when the value of decreases to a critical value. This is the same behaviour observed for the full model when varying any of the parameters.

Discussion
Our mathematical model of fibroblast–epithelial cross-talk exhibits bistability for certain parameter values. The cross-talk model is formulated as a system of time-dependent ODEs without explicit spatial heterogeneity. Instead, a control parameter representing the influx of BMPi, , serves as a surrogate for distance from the base of the crypt. The different steady-state solutions admitted by the model as varies, represent mutually supporting fibroblast and epithelial cell phenotypes. We identify distinct zones of the crypt, characterised by the number and nature of the steady-state solutions and their local stability. Based on our analysis of the steady-state and time-dependent solutions of the model, we have made the following principal findings:Spatially distinct phenotypic zones arise as a consequence of bistability: We assume that the steady states of the model represent pairs of mutually supporting fibroblast and epithelial cell phenotypes; thus, these phenotype pairs exist along a continuum as varies. However, the bifurcation analysis of the model enables a broad classification of all pairs as belonging to two possible groups of phenotypes defined by the two stable solution branches: one corresponding to epithelial stem cells and their supporting fibroblasts, and the other to differentiated epithelial cells and their supporting fibroblast phenotype. These two phenotype pairs are distinct, as no stable solution branch connects them. Consequently, the system does not permit intermediate phenotypes. This suggests that, although a range of fibroblast–epithelial phenotype pairs exist within the crypt, all fall into one of these two functional groups. The bifurcation diagrams of the steady-state solutions as varies delineate three distinct intervals, corresponding to phenotypic zones within the crypt: Stem, Transitional, and Differentiated. The Stem and Differentiated zones only contain the stem and differentiated epithelial cell solution branches, respectively. The Transitional zone lies between these two, defined by the interval of in which both stable branches coexist with an additional unstable branch. The boundaries between these zones are marked by fold bifurcations, which occur at critical values of and are characterised by higher or lower BMP and WNT concentrations. The BMP low/WNT high fold bifurcation defines the boundary between the Differentiated and Transitional zones, while the BMP high/WNT low fold bifurcation marks the boundary between the Transitional and Stem epithelial cell zones. Because serves as a surrogate for distance from the base of the crypt, these zones are also spatially distinct: the Stem zone resides at the crypt base, the Differentiated zone at the top, and the Transitional zone lies between them. This prediction is consistent with experimental observations that fibroblasts at the base and top of the crypt adopt distinct phenotypes that support stem and differentiated epithelial cells, respectively (McCarthy et al 2020).

Distinct phenotypes emerge during epithelial cell migration: Time-dependent solutions of the governing equations reveal a hysteresis loop as is varied across the interval of BMPi defining the Transitional epithelial cell zone (bistable region), according to a prescribed function of time. If the time-dependent solution starts on the epithelial stem cell steady-state solution branch, as decreases, the time-dependent solution switches to the solution branch associated with the differentiated epithelial cells and their supporting fibroblasts. The transition from the stem to the differentiated phenotype is triggered when decreases past the critical value marking the boundary between the Differentiated and Transitional zones. Conversely, the transition from differentiated to stem phenotypes occurs when increases past the larger critical value defining the boundary between the Transitional and Stem epithelial cell zones. A natural interpretation of a time-varying is the changing distance of an epithelial cell from the BMPi-emitting mesenchyme layer as it migrates up the crypt. Thus, the switch of the time-dependent solution between steady-state solution branches models epithelial stem cells differentiating once they reach a certain height up the crypt. Importantly, the model also predicts that differentiated epithelial cells may, in principle, de-differentiate into stem cells. This behaviour aligns with experimental observations of ectopic crypt formation in response to increased BMPi concentrations (Haramis et al 2004). However, the hysteresis loop of the model predicts that de-differentiation is not triggered by natural fluctuations of the expression of BMPi by the cells in the mesenchyme below the crypt.

Hhog is essential for rapid phenotype switching during epithelial cell migration: A central component of our model analysis involves identifying the most critical cross-talk interactions using informative submodels. In our model, cross-talk is represented by modulation of a morphogen’s expression rate proportionally to the proportion of a bound morphogen receptor. One key submodel studied in this article is the No BMP stimulus model, in which the hhog-mediated stimulus of BMP is removed by setting its associated production rate to zero (see Section 3.2.1). While this submodel still admits all three phenotypic zones, the fold bifurcation separating the Differentiated and Transitional epithelial cell zones is replaced by a transcritical bifurcation. As decreases past this bifurcation, the time-dependent solution transitions from the epithelial stem cell solution branch to the differentiated epithelial cell solution branch, but at a slower rate than the full model. This slow transition likely reflects prolonged retention of the epithelial stem cell phenotype during migration. Although this result shows that fibroblasts can, in principle, regulate epithelial fate without cross-talk, such cross-talk is necessary for the rapid phenotype switching observed in vivo. Therefore, hhog signalling ensures the zones of the fibroblasts and epithelial cell phenotypes predicted by the steady-state equations for spatial static populations of fibroblasts and epithelial cells are present for a static fibroblast population with a dynamic epithelial cell population which migrates up the crypt.

Overactive WNT inhibition prevents fast switching of epithelial cell phenotypes: Our analysis indicates that overactive WNT-mediated hhog inhibition can also prevent rapid phenotype switching as epithelial cells migrate up the crypt (Section 3.2.2). As for the No BMP stimulus model, the boundary between the Differentiated and Transitional epithelial cell zones is defined by a transcritical bifurcation. Therefore, the transition of the time-dependent solution between the stem and differentiated epithelial cell solution branches is much slower than for the full model.

Most cross-talk interactions are redundant: After identifying the key cross-talk interactions required to generate the three phenotypic zones and enable fast switching between stable steady-states, we assessed how many interactions could be removed without disrupting normal crypt behaviour. Out of the eight cross-talk interactions in the full model, our results indicate that only three are essential to preserve the core behaviours. Specifically, the model retains its qualitative behaviour provided; there is at least one interaction in which each of BMP, hhog, and WNT acts as either a stimulus or an inhibitor; there is at least one interaction in which each of BMP, hhog, and WNT is a product morphogen for a cross-talk interaction; and both hhog-mediated BMP stimulus and WNT-mediated hhog inhibition are present. Only the cross-talk interactions of the minimum viable model meet these criteria. The presence of the three phenotypic zones separated by two fold bifurcations suggests that the remaining five interactions can be omitted, with their functional roles compensated for by adjusting the rate constants of the remaining three interactions. Our model assumes that experimentally observed correlations in morphogen expression are the result of direct receptor-mediated interactions, specifically, morphogen binding or unbinding events. While this assumption may not hold in all cases (some correlations may arise from series of interactions as opposed to a single interaction), our analysis shows that such simplifications do not compromise the qualitative predictions of the model. Therefore, the model’s behaviour remains robust even when the assumption of direct regulation for all correlated morphogens is relaxed.

The existence of three phenotypic zones depends on a delicate balance of cross-talk interactions: In addition to analysing the steady-states of informative submodels, we investigated how the bifurcation structure of the full model changes as the values of the parameters governing cross-talk interactions vary. In Section 3.3, we selected a parameter set for which the model exhibits three phenotypic zones, each defined by fold bifurcations as varies. We then examined how this structure changes as the rate constant for hhog-mediated BMP stimulus is varied. Two distinct types of disruption to the bifurcation structure were observed. The first is the meeting of the two fold bifurcations at a cusp point, resulting in a single monostable zone with a unique steady-state solution for all values of . The second is a shift in the positions of the fold bifurcations, which alters the boundaries between phenotypic zones. If the value of at the BMP low/WNT high bifurcation is negative, the Differentiated epithelial cell zone is lost. If the value of at the BMP high/WNT low bifurcation is negative, both the Differentiated and Transitional epithelial cell zones are lost. We extended this analysis to all other cross-talk parameters and found that both types of disruption, cusp formation and boundary migration, can be induced by varying any of the cross-talk rates, with one exception. In some cases, altering the cross-talk rate resulted only in the loss of the Differentiated epithelial cell zone, but this alone is sufficient to prevent phenotype switching.

Mutations in either fibroblasts or epithelial cells alter the phenotype of the opposite cell population: A possible interpretation of changes in the rate constants explored in Section 3.3 is that they reflect the effects of cancerous mutations in either epithelial cells or fibroblasts. Variations in these rate constants shift the boundaries between phenotypic zones, thereby altering the spatial regions within the crypt where viable epithelial cell and fibroblast phenotype pairs can be maintained. Because these disruptions can be triggered by changes in either fibroblast- or epithelial-specific parameters, the model predicts that a mutation in one cell population is sufficient to alter the phenotype of the other. In particular, if the Differentiated epithelial cell zone is lost due to a mutation in either cell type, epithelial cells no longer undergo phenotype switching as they migrate up the crypt. Although proliferation is not included in our model, preventing epithelial cells from differentiating would shift the crypt toward sustained stem-like states. This loss of differentiation is known to drive hyperproliferation and the emergence of adenomatous polyps, a hallmark of early CRC lesions.

Our results demonstrate a possible mechanistic explanation for the formation of distinct fibroblast populations, that epithelial cell fate, driven by implicit spatial variation and epithelial-fibroblast cell cross-talk. Furthermore, our model shows that alterations in cross-talk interactions, potentially resulting from mutations in either epithelial cells or fibroblasts, can influence the phenotype of the opposing cell population. A natural extension to our findings would be to test if our model predicts the existence of distinct phenotypic populations with explicit spatial variation using Agent-Based Models (ABMs). ABMs are widely used to study the dynamic behaviours of epithelial cells in intestinal crypts (Bravo and Axelrod 2013; Buske et al 2011; Du et al 2015; Dunn et al 2012, 2016; Ingham-Dempster et al 2017; Van Leeuwen et al 2009; Ward et al 2020). In these models, epithelial cell fate is typically governed by local WNT concentrations; however, they often neglect the contributions of fibroblasts or fibroblast–epithelial cross-talk. These models often use fixed WNT gradients (Buske et al 2011; Dunn et al 2016; Van Leeuwen et al 2009), although some models, such as those by Du et al (2015) and Ward et al (2020), incorporate WNT gradients generated by epithelial cells. Of these, only Du et al (2015) considers the role of a fixed BMP gradient. As such, existing models are limited in that they omit the dynamic role of plastic fibroblast phenotypes and their critical involvement in the establishment of morphogen gradients that regulate epithelial cell fate. Incorporating our cross-talk model into an ABM framework could address these limitations. By eliminating the use of an abstract control parameter, and instead modelling spatial dynamics explicitly, we could evaluate whether our model still predicts the emergence of spatially distinct fibroblast populations. Furthermore, this approach would also enable investigation into whether variations in cross-talk interaction parameters can realistically simulate the crypt’s response to cancerous mutations, offering a powerful tool for understanding disease progression and therapeutic targeting.
Another avenue for future research is coupling the cross-talk model with models of subcellular signalling pathways associated with the BMP, WNT and hhog morphogens. In the current framework, morphogen expression rates are modelled as being directly proportional to the proportion of bound receptors on fibroblasts and epithelial cells. However, in reality, these rates are likely to depend on subcellular dynamics. For example, Lee et al (2003) proposed a model of the WNT pathway, in which the transcription of target genes depends on the concentrations of several subcellular proteins. MacLean et al (2015) extended this model to show how the WNT pathway can admit bistability for some combinations of parameter values, and, by extension, separate phenotypes. Extending the crypt model by integrating equations for the morphogen pathways such as those for WNT proposed by Lee et al (2003) or MacLean et al (2015) could result in the existence of more than two stable steady-states, and by extension more than two phenotypic pairs. Our model classifies all fibroblasts and epithelial cells as belonging to two broad classes of phenotype pairs corresponding to each stable steady-state solution branch. More than two stable steady-state solution branches would allow for the discrimination of different phenotypes within the two groups predicted by our model. For example, differentiating between epithelial stem cells in the stem cell niche and transit-amplifying cells.

Conclusions
To our knowledge, the model presented in this article is the first to describe the impact of cross-talk between fibroblasts and epithelial on the formation of multiple populations of mutually supporting epithelial cell and fibroblast phenotypes. The model proposes a mechanism by which the morphogen gradients may form within the crypts. The model predicts that different phenotypic fibroblast populations responsible for epithelial cell fate observed by McCarthy et al (2020), are an emergent property of epithelial-fibroblast cross-talk and spatial heterogeneity. The biological significance of this result is that fibroblasts do not attain permanent phenotypes within the crypt, but instead modulate their phenotype depending on their proximity to the mesenchyme beneath the crypt. These results are consistent with the view that fibroblasts and epithelial cells have co-evolved to create the crypt structure and regulatory mechanisms, as suggested by Shyer et al (2015). Furthermore, our analysis shows how mathematical modelling can be used to test the impact of mutations on the pathways of morphogens responsible for the fate of the epithelial cells and fibroblasts within the crypt. In future work, we will embed the model within an ABM to examine how mutations in epithelial cells or fibroblasts reshape the local microenvironment. By capturing these early dysregulations, the extended model will enable us to simulate the microenvironmental conditions that precede polyp formation and contribute to CRC initiation.