A Hybrid Model for Cellular Dynamics in Colorectal Crypts

Sousa, I. J. L.; Romanazzi, G.

DOI 10.5433/1679-0375.2025.v46.53590

Citation Semin., Ciênc. Exatas Tecnol. 2025, v. 46: e53590

Received: August 17, 2025 Received in revised for: October 24, 2025 Accepted: November 12, 2025 Available online: December 12, 2025

Abstract:

Colorectal cancer is believed to originate from abnormal cell proliferation in small cavities of its epithelium called crypts. We present a hybrid framework that couples a Voronoi tessellation-based cell model with a continuous differential model to describe cellular dynamics, proliferation, and differentiation in a colorectal crypt. The framework is implemented numerically using finite difference and finite element methods. This allows us to simulate how cells move, proliferate, and differentiate in normal and abnormal scenarios. We consider in particular the normal homeostatic case for validating the framework, and simulate two abnormal scenarios where perturbations in cell distribution or proliferation are imposed. The benefit of such a hybrid model is that we can freely impose perturbations and show results for cellular dynamics, proliferation, and differentiation at two different scales: a macro continuous scale and a micro cell scale. This framework provides insights into how cell behavior changes under abnormal conditions and how crypt lesions and adenomas may originate from cellular dynamics perturbations.

Keywords: colon crypt, partial differential equations, tessellation, cell proliferation and differentiation, finite element method

Introduction

A tightly regulated process of continuous cell renewal and differentiation in colonic crypts plays a crucial role in maintaining healthy digestion. Disruptions in this process, such as abnormal cell proliferation within a crypt, are associated with the onset of colorectal cancer (Osborne et al., 2015; Vanleeuwen et al., 2006). To deepen our understanding of how such abnormalities evolve into carcinomas, it is essential to model and simulate crypt cellular dynamics under both physiological and pathological conditions.

In the field of mathematical modeling, two main approaches are widely used to describe cellular behavior in biological tissues: discrete based-cell and continuous models. Discrete models, such as the cell-centered Voronoi Tessellation (VT) model, focus on the spatial location and mechanical interactions of individual cells (Meineke et al., 2001). In contrast, continuous models describe the dynamics in space and time of macroscopic quantities, such as cell density, using systems of differential equations (Olivrom et al., 2022).

Motivated by the need for a more comprehensive representation of crypt dynamics, we propose a hybrid cellular model (HCM) which integrates both modeling strategies to describe cell dynamics in a single crypt. The HCM couples a continuous partial differential equation (PDE) system with a discrete VT model. While the continuous model describes the pressure as well as the dynamics, proliferation, and differentiation of crypt cell families, the VT model accounts for the spatial location, mechanical forces and adhesion of each cell. In particular the discrete framework is used to update each cell according to the cell densities provided by the continuous model, enabling the simulation of both local mechanical interactions and global population-level behaviors.

The continuous model is solved numerically using a piecewise linear finite element method, which provides an approximation of cell family densities at each mesh node on the crypt epithelium surface. Then, an integration step updates the distribution for each VT cell allowing the VT model implementation.

The proposed HCM allows us to simulate both normal and abnormal proliferation scenarios. It offers valuable insights into how disordered patterns of cell growth can potentially lead to adenoma or lesion formation that may progress to cancer. In fact, it is known that abnormal crypt cell proliferation patterns can be responsible for the early stages of colorectal cancer (Vanleeuwen et al., 2006; Romset et al., 2022).

Materials and methods

Colorectal crypts are formed by different cell types whose differentiation and proliferation follow a regular process that is important to model since, as mentioned before, any process disruptions may lead to early stages of colorectal cancer.

Here we classify only three cell types (also called families) by listing their characteristics. However, more cell families exist, such as the ‘Paneth’ cells, that will be then immersed in our proposed three family model. The three cell families modeled are: stem cells, semi-differentiated cells and fully differentiated cells (Murray et al., 2011).

Stem cells are characterized by dividing and generating descendant cells that differentiate into other cell types. Semi-differentiated cells also undergo a differentiation process (called specialization) to perform specific and vital functions in the intestine. Moreover, similar to stem cells, they can proliferate (divide) but faster than stem cells.

The last cell family is formed by fully differentiated cells that are mature cells completely specialized. They perform vital functions in the epithelium tissue despite lacking the ability to generate new cells of the same lineage through cell division, thus fully-differentiated cells do not proliferate.

The behavior of all these cells is crucial to maintain a healthy intestinal homeostasis (Asarian et al., 2012), in fact stem cells ensure a continuous supply of differentiated cells to maintain the intestinal lining. Instead, semi-differentiated and fully differentiated cells perform vital functions necessary for digestion and nutrient absorption, as well as for protecting the intestinal mucosa from external aggression.

It is known that crypts can be represented geometrically as a cylinder (Kershaw et al., 2013), with an open top orifice and a closed rounded bottom see Figure 1.

(a)

(b)

Crypt epithelium cells are located on the crypt surface, as seen in the endoscopic pictures in Figure 1(a). Note that stem cells are located bottom of the crypt, whereas semi-differentiated and fully differentiated cells are located in the middle and top of the crypt, respectively. The crypt geometry is simplified by removing its closed bottom thus without considering stem cells located at the very crypt bottom. Then the cylinder is unfolded to obtain a rectangle as shown in Figure 1(b). This simplified representation is commonly used in works involving Voronoi tessellation models (Meineke et al., 2001; Vanleeuwen et al., 2009).

Let \(\Omega=\left\{ (x,y)\in\mathbb{R}^2:\; 0<x<c,\; 0<y<h \right\}\) be the rectangular crypt domain, where \(c\) and \(h\) are the orifice circumference and the crypt height, respectively. We denote its boundary by \(\Gamma\), see Figure 2, where \(\Gamma_1\) is its bottom edge, \(\Gamma_3\) and \(\Gamma_4\) are the left and right edges, respectively, and \(\Gamma_2\) is its crypt top orifice edge.

(a)		(b)		(c)

Continuous model

The continuous model in HCM describes the crypt cellular dynamics in \(\Omega\) over the time interval \([0,T]\). It is based on a system of partial differential equations (PDEs) with convection, diffusion and reaction terms that has been modified with respect those used in (Figueiredo et al., 2019; Olivrom et al., 2022) to add an equation for stem cell dynamics and a more well-balanced semi-differentiated cell proliferation.

The unknown system solutions, defined at the points \((x, y) \in \Omega\) and at time \(t\in (0, T]\), are the densities of stem, semi-differentiated, fully differentiated cells denoted respectively by \(C_{\ell}(x, y, t)\), with \(\ell=A,B,C\), and the cell pressure \(p(x, y, t)\).

This pressure is generated by the cell mitotic activity which induces a passive cell displacement in the crypt with velocity \(\mathbf{v}=-\xi \nabla p\) usually directed upwards to the crypt top orifice. Here \(\xi\) is a tissue parameter depending on the permeability and viscosity of the epithelium which may depend on the cell (Greenspan et al., 1976). Notwithstanding each cell type may have different diffusion \(D\) and tissue parameter function \(\xi\), we assume for model simplification that they are independent of cell type.

Based on these assumptions, the space-time continuous model for cell densities \(C_A\), \(C_B\), \(C_C\) is governed by the following convective-diffusive-reaction PDE system

The rates \(\alpha_B\), \(\beta_{AB}\), and \(\beta_{BC}\) represent, respectively, the proliferation rate of semi-differentiated cells (\(B\)), the differentiation rate from stem (\(A\)) to semi-differentiated cells, and the differentiation rate from semi-differentiated to fully differentiated cells (\(C\)). All these rates are functions depending on the space coordinate \(y\) along the vertical crypt axis. Proliferation rate \(\alpha_B\) increases from the crypt bottom (\(y=0\)) reaching a maximum at an intermediate height, and then decreases to zero at the crypt top \((y=h)\) (Almet et al., 2021).

Experimental observations (Nicolas et al., 2007) indicates that stem cell differentiation occurs after asymmetric stem cell divisions producing one stem cell and one semi-differentiated cell. Thus stem cell population is constant over time under physiological conditions. This observation justifies that differentiation term \(\beta_{AB} C_A\) is present only in the second equation of the system of equations (1).

Both rates \(\beta_{AB}\) and \(\beta_{BC}\) progressively increase towards the crypt top and are negligible at the crypt bottom. This spatially dependent behavior is consistent with the literature (Drasdo et al., 2001; Di et al., 2010), which report a predominance of differentiation far from the crypt bottom where only stem cells reside. We assume that no void exists in the epithelium crypt and then, since \(C_A\), \(C_B\), and \(C_C\) represent volumetric fractions of the respective cell types in the point \((x,y)\) at time \(t\), the following constraint is valid

This implies that domain \(\Omega\) is fulfilled only by three cell families. Then, we can eliminate the third equation in the system of equations (1) and rewrite it in terms of \(C_A\) and \(C_B\) only. Summing the resulting two equations yields, after simplifications, the following system

The following boundary and initial conditions complete the system:

The system of equations (3) with boundary and initial conditions in equations (4) and (5) is the final version of our continuous model accounting for cell diffusion, pressure-driven convection, proliferation, and differentiation. This is solved to get \(C_A\), \(C_B\) and \(p\) in the crypt \(\Omega\) at each time \(t\in[0,T]\), and then, by using the no-void relation in equation (2) which ensures the conservative nature of the model, we can simply get the stem cell density \(C_C=1-C_A-C_B\).

Due to its construction, see Figure 2(a), domain \(\Omega\) has periodic boundary conditions on \(\Gamma_3\) and \(\Gamma_4\) and then we have continuity for density and pressure on these edges. This is imposed by the homogeneous Neumann conditions on pressure and cell densities in equations (4) and (5). At the crypt top \(\Gamma_2\), the Neumann pressure condition \(\xi \frac{\partial p}{\partial y} = - g_y\) uses a known experimental velocity \(g_y\) obtained by histological observations (Vanleeuwen et al., 2006).

Boundary conditions on the crypt bottom \(\Gamma_1\) are obtained instead based on the known predominance of stem cells therein. This means that the crypt bottom is an auto-regenerative tissue with a relevant internal renewal. While in the crypt middle most of cells are semi-differentiated, fully differentiated cells with density \(C_C\) predominate at the crypt top. These mature cells are then shed to the lumen and replaced by fully differentiated cells coming from the crypt middle in a continuous homeostatic renewal process. Analytical expressions for the initial densities \(C_A^0\), \(C_B^0\), \(C_C^0\) in the Normal and Abnormal Cases can be found in the Results and discussion section.

As it will be shown therein, the continuous model given in equations (3)–(5), will allow faithful reproduction of physiological patterns that are consistent with experimental and theoretical assumptions of crypt cell renewal processes in normal and abnormal scenarios.

Voronoi Tessellation model

The discrete component of HCM uses the cell-centered Voronoi Tessellation (VT) model to obtain the locations and mechanical interactions of each crypt cell. For each cell \(k\), its geometric center, denoted by \(\mathbf{q}_k\), is also referred as the VT generator point.

Mechanical interactions, also mentioned as forces, between neighboring cells are originated by two main mechanisms: cell adhesion and migration towards preferred regions.

In particular, a force between two neighboring cells acts along their connecting line and depends on the distance between the cell centers. Cell adhesion is modeled as follows: if neighboring cells are very close, the force is repulsive, whereas at intermediate distances an attractive adhesion predominates. In this model, the force between two cells is considered to be linearly proportional to the overlap of their volumes, defined as the difference between the sum of their natural radii and the distance between their centers. This simplification allows capturing the dominant mechanical effects of cell adhesion without introducing excessive modeling complexity. Thus the force acting on cells \(k\) due to cell \(\ell\) is

where \(\mathbf{q}_k\) and \(\mathbf{q}_\ell\) are respectively the position vectors of the \(k\) and \(\ell\) cell centers, \(\mu\) is a tissue stiffness and

is the overlap between the cells \(k\) and \(\ell\) that have radius \(R_k\) and \(R_\ell\) respectively.

The total force \(\mathbf{F}_k\) acting on cell \(k\) is obtained by summing forces \(\mathbf{F}_{k\ell}\) in equation (6) generated by cells \(\ell\) adjacent or in contact with cell \(k\), that is

where \(\mathcal{N}_k\) is the set of cells adjacent or in contact with cell \(k\).

Assuming that cell dynamics occurs in an over-damped regime, which is typical in viscous domains such as the colon epithelium, inertial effects can be neglected, and the equation governing the cell \(k\) dynamics is:

where \(\eta\) is the viscosity coefficient associated with the medium’s resistance to cell displacement.

Equation (9) is solved numerically using the explicit Euler method, which updates the cell position at time \(t + \Delta t\), for \(\Delta_t>0\), as follows

This is done for a fixed number \(N_t\) of time steps up to covering the chosen simulation time interval.

In Voronoi Tessellation models, cells are represented by polygonal regions determined by the cell center positions (Masteritalo et al., 2024), and connectivity is inferred from the Delaunay triangulation, the dual of the Voronoi mesh. However, when a purely linear force is used as done above, unrealistic connections may exist between distant cells, compromising biological realism. To avoid this, we adopt a cutoff linear modification, in which the force is neglected if the overlap \(d_{k\ell}\) defined in equation (7) is below a negative threshold \(d_{\text{min}}\). This modification makes the model more robust, preventing spurious interactions and ensuring a more accurate representation of tissue organization and mechanics throughout the simulation. Thus, we have \[F_{k\ell} = \begin{cases} \mu d_{k\ell}+(0,-\dfrac{d}{dy_k} (y_k-y_{\tau_k}^2)), & \text{if } d_{k\ell} \geq d_{\text{min}} \\ (0,-\dfrac{d}{dy_k} (y_k-y_{\tau_k}^2)) & \text{if } d_{k\ell} < d_{\text{min}} \end{cases}, \quad \text{ with } y_{\tau_k}=\left\{\begin{array}{cc} \frac{1}{6}h & \text{ if }\tau_k=A \\ \frac{1}{2}h & \text{ if }\tau_k=B \\ \frac{5}{6}h & \text{ if }\tau_k=C \end{array}\right.,\] where the added gradient term \((0,-\dfrac{d}{dy_k} (y_k-y_{\tau_k})^2)\) mathematically models the experimental observation that stem cells (\(A\)) prefer to be located at the lower third crypt level near the crypt bottom, semi-differentiated cells (\(B\)) in the middle and mature fully differentiated cells (\(C\)) at the upper third crypt level.

Numerical method for the continuous model

The initial boundary value problem given in equations (3)-(5) is numerically solved by using a piecewise finite element method (FEM), implemented in Python using FEniCS library (Fenics et al., 2015) with the spatial discretization carried out using a triangular mesh, see Figure 2(b), and time integration is performed using a first-order implicit-explicit scheme. We give below some overview of the piecewise linear finite element method implemented.

Let \(\Delta_T=\frac{T}{N_T}>0\) the time step of the \(N_T\) intervals \([t_{n-1},t_n]\) with \(t_n=n\Delta_T\), for each \(n=1,\ldots,N_T\), our iterative method obtain a numerical pressure approximation \(p^{n-1}_h\) for the pressure \(p(\cdot,\cdot, t_{n-1})\) at time \(t_{n-1}\) and a after compute the numerical densities \(C^{n}_{A,h}\) and \(C^{n}_{B,h}\) that approximate \(C_A(\cdot,\cdot,t_n)\) and \(C_B(\cdot,\cdot,t_n)\) which are solutions of equations (3)-(5) at time \(t_n\). Let also consider the space \(V=\left\{v \in H^{1}(\Omega): ~v|_{\Gamma_{1}}=0\right\}\), we obtain the variational formulation at time \(t=t_{n-1}\) of (3)\(_1\), that is of the first equation of the system of equations (3), as follows: find \(p\in V\) such that

where \[\begin{aligned} \small a_p(u,v)= \int_\Omega \xi \nabla u\nabla v dxdy,\quad \ell^{n-1}_p(v)=\int_\Omega \left(\alpha_BC^{n-1}_B(1-C^{n-1}_B)+\beta_{AB}C^{n-1}_A \right)v dxdy -\int_{\Gamma_2}\xi g_yvds \end{aligned}\] for all \(u,v\in V\), where \(C_A^{n-1}\), \(C_B^{n-1}\) are approximations of \(C_A(\cdot,t_{n-1})\), \(C_B(\cdot,t_{n-1})\) obtained in the previous iteration \(n-1\). Let \(V_h\subset V\) the subspace of piecewise linear functions on the triangulation \(\Omega_h\subset\Omega\), see Figure 2(b), with basis \(\{\Phi_k\}_{k=1,\ldots,N_h}\) where \(N_h\) is the number of \(\Omega_h\)-nodes except that in \(\Gamma_1\). The piecewise linear finite element method consists in finding \(p^{n-1}_h=\displaystyle\sum_{k=1}^{N_h}p^{n-1}_k\Phi_k\in V_h\), an approximation of the solution \(p\) of equation (11), such that This reduces to finding \(\overline{p}^{n-1}=\left(p^{n-1}_i\right)\in\mathbb{R}^{N_h}\) solution of \(A\overline{p}^{n-1}= b\) where \(A=(A_{i,j})\) is the \(N_h\times N_h\) and \(b=(b_i)\) is the \(N_h\times 1\) vector, defined as \(A_{i,j}= a_p(\Phi_i,\Phi_j)\) and \(b_i=\ell^{n-1}_p(\Phi_i)\) for \(i,j=1,\ldots,N_h\).

We now construct a time-discrete Backward Euler and space-discrete piecewise linear finite element method for approximating \(C_A\) at time \(t_n\) by building the variational formulation of the second equation of the system of equations (3).

Let \(W=\{w\in H^{1}(\Omega): ~w|_{\Gamma_1\cup\Gamma_2}=0\}\) and \(W_h\subset W\) the \(W-\)subspace of linear piecewise functions in \(\Omega_h\) with basis \(\{\Psi_k\}_{k=1,\ldots,M_h}\), where \(M_h\) is the number of \(\Omega_h-\)nodes except those in \(\Gamma_1\cup\Gamma_2\). Using Backward Euler method in time, the piecewise linear finite element method applied to the second equation of the system of equations (3) at each time step \(t_n\) seeks \(\widetilde{C}^{n}_{A,h}\in W_h\) such that

where \[\begin{aligned} a^n_A(u_h,w_h)= \int_{\Omega} D\nabla u_h \nabla w_h +\xi u_h\nabla p^{n-1}_h \nabla w_h dxdy, \quad \forall u_h,w_h \in W_h \end{aligned}\] with \(p^{n-1}_h\) obtained in equation (12).

The problem in equation (13) reduces to finding \(\overline{C}^n_{A}=\left( \widetilde{C}^n_{A,i}\right)\in\mathbb{R}^{M_h}\) such that \((M+\Delta_t S_A^n)\overline{C}^n_A = M\overline{C}^{n-1}_A\), where \(M=(M_{ij})\) and \(S_A^n=(S^n_{A,i,j})\) are \(M_h\times M_h\) matrices whose coefficients are \(M_{ij}=\int_{\Omega}\Psi_i\Psi_j dxdy\), and \(S^n_{A,i,j}=a^n_A(\Psi_j,\Psi_i)\). To obtain an approximation \(C_{A,h}^n\) of \(C_A\) at time \(t_n\) that satisfies the conditions \(C_A=1\) if \(y=0\) and \(C_A=0\) if \(y=h\), we need to add the function \[\Psi^{n-1}_p(x,y)=\left\{\begin{array}{ll} \displaystyle e^{-\frac{\xi}{D} p^{n-1}_h(x,y)} & \text{ if } y\ne 0\\ 0 & \text{elsewhere} \end{array}\right.\] to the solution \(\widetilde{C}^{n}_{A,h}\) of equation (13). Thus, we get

as numerical approximation of \(C_A\) at time \(t_n\). Finally we need to approximate numerically the density \(C_B\) at time \(t_n\) using the nonlinear third equation of the system given in equation (3).

This is achieved using a time-discrete implicit-explicit scheme combined with a space-discrete piecewise linear finite element method: find \(C_{B,h}^n\in W_h\) such that

where \(\forall u_h,v_h,w_h \in W_h\): \(a^n_B(u_h,v_h,w_h)= \int_{\Omega} D\nabla u_h \nabla w_h + \xi u_h\nabla p^{n-1}_h \nabla w_h -(\alpha_B-\beta_{BC})u_hw_h + \alpha_B u_hv_h w_h ~dxdy\) and \(l^n_B(w_h)=\int_{\Omega}\beta_{AB}C^n_{A,h} w_h dxdy,\) with \(C_{A,h}^n\) defined as the numerical approximation obtained in equations (13) and (14).

The problem in equation (15) reduces to finding \(\overline{C}^n_{B,h}=\left( C^n_{B,i}\right)\in\mathbb{R}^{M_h}\) such that \((M+\Delta_t S_{B}^n)\overline{C}^n_{B,h} = M\overline{C}^{n-1}_B\), where \(M=(M_{ij})\) and \(S_B^n=(S^n_{B,i,j})\) are \(M_h\times M_h\) matrices whose coefficients are\(M_{ij}=\int_{\Omega}\Psi_i\Psi_jdxdy\) and \(S^n_{B,i,j}=a^n_B(\Psi_j,C^{n-1}_{B,h},\Psi_i)\), respectively.

Hybrid cellular model

The hybrid cellular model (HCM) couples the Voronoi Tessellation model with the continuous model which have been described in the previous sections. This is done by using the following iterative HCM algorithm:

• Step 1. Solve numerically the continuous problem in equations (3)-(5), to obtain pressure \(p\) at time \(t=t_n\) and cell densities \(C_A\), \(C_B\), \(C_C\) at time \(t_{n+1}=t_n+\Delta_T\), as described in the previous section.

• Step 2. Project the continuous densities to the VT generator points \(\mathbf{q}_k\) and select the cell type for each cell.

• Step 3. Compute the interaction forces between cells as defined in section Voronoi Tessellation Model and update their locations using equation (10) for \(N_t\) steps with time step \(\Delta_t\).

• Step 4. Get the continuous density interpolating the VT cell positions of previous steps on the triangulation \(\Omega_h\)

• Step 5. Update \(n\): \(n+1\to n\).

Let \(T>0\) the chosen simulation time, the algorithm Steps are repeated \(N_T\) times until reaching time \(t=T\), where \(N_T=T/\Delta_T\). Since Step 1 has been fully explained in the section Numerical Method for the Continuous Model, we describe here with further details the remaining HCM algorithm steps.

In Step 2, we project the cell density from the continuous model to the VT model built on the Voronoi generator points \(\mathbf{q}_k\) that coincide with the \(\Omega_h\) triangulation nodes. In order to determine the correct \(k\) cell type we choose the type \(A\), \(B\), or \(C\) that has the maximum continuous density in the generator point \(\mathbf{q}_k\).

In Step 3, a Forward Euler method with time step \(\Delta_t\) is applied to update \(\mathbf{q}_k\) according to the adhesion forces described in the section Voronoi Tessellation Model. Note that time step \(\Delta_t\) in the VT model is smaller than the time step \(\Delta_T\) used in the continuum model, since cellular adhesion occurs on a shorter timescale (hours) compared to differentiation or proliferation (days). Note also that the VT step number \(N_t\) should be sufficiently large to guarantee that a stable/homeostatic solution of equation (9) is reached. Since cells \(k\) have changed their location in order to update the continuous densities \(C_{\tau}(x,y)\) with \(\tau=A,B,C\) in the triangulated mesh \(\Omega_h\) we perform in Step 4 an interpolation that averages VT cell density in small square regions \(s_{ij}\) centered in the \(\Omega_h\) nodes \((x_i,y_j)\) as follows:

\[C_{\tau}(x_i,y_j)= \dfrac{\text{Sum of cell areas of type } \tau \text{ that intersect the square } s_{ij}}{|s_{ij}|} ,\quad \tau=A,B,C\]

where \(|s_{ij}|\) denotes the area of \(s_{ij}\).

A graphical representation of these squares \(s_{ij}\), highlighted in yellow, on a VT cell distribution with regular hexagonal cells is illustrated in Figure 2(c). Herein the \(s_{ij}\) centers are marked by black dots and coincide with the \(\Omega_h\) nodes located on the edges of the triangulation \(\Omega_h\) represented in Figure 2(b).

Finally in Step 5 we update the iteration time index \(n\), thereafter we obtain at time \(t_n\) the discrete VT cell distribution and the continuous cell densities in all the crypt domain \(\Omega\).

Results and discussion

Here, we present results from the implementation of HCM described in the previous section. In order to simulate cellular dynamics based on experimental data we begin to describe the cell and crypt geometrical dimensions. Crypt cells have an average spherical shape with diameter \(5.9 \mu m\) (Guebeltorres et al., 2008).

The average normal crypt has a crypt orifice surrounded by 23 or 24 cells (Baker et al., 2014) and height formed by a number of cells between 75 and 110 (Bernstein et al., 2010). Assuming an average of 23.5 cells along the crypt orifice circumference \(c\), its length is \(c=23.5\times 5.9= 138.65\mu m\). Whereas the crypt height \(h\) is formed on average by \(92.5=\frac{75+110}{2}\) cells, then the crypt height measures \(92.5\times 5.9= 545.75 \mu m\).

Thus, the relation between height and circumference is \(\frac{h}{c}=\frac{545.75}{135.7}\approx 4.02\approx 4\) and the crypt can be discretized using \(c=23\) along the circumference and \(h=92\) cells along the height. However, this would yield to a total of \(23\times 92=2116\) crypt cells, which is computationally expensive.

To reduce computational cost, each VT cell is represented by a compartment containing \(3\times 3\) actual crypt cells, leading to a discretized crypt with \(c=\frac{24}{3}=8\) and \(h=\frac{90}{3}=30\) cells.

We simulate and track crypt cellular dynamics in \(\Omega\) under three scenarios:

• Normal Case: Parameters, initial conditions and rates are adjusted to reflect typical cell homeostasis in a crypt.

• Abnormal Case 1: The semi-differentiated cell proliferation rate \(\alpha_B\) and top velocity are increased compared to the Normal Case.

• Abnormal Case 2: Perturbed initial conditions for \(C_A\), \(C_B\), and \(C_C\) are considered, where stem cells are also present in middle of the crypt in addition to their expected location at the crypt bottom.

We first tested our model under the Normal Case to validate its performance. This task is important since the most common situation in a crypt is that the cellular dynamics is a homeostatic process. Any perturbation of this case will be considered as an Abnormal Case and then this will be used to simulate and understand how adenomas or lesions can be originated and evolve in the crypt. Note that our test-cases have been built for understanding adenoma formation since the exact mechanism by which a disruption in a single crypt leads to adenoma formation, and subsequently to carcinoma in the colorectal tract, remains unclear (Vanleeuwen et al., 2006). The results of the HCM implementation for these three cases are presented and analyzed in the following sections.

Normal Case

In the Normal Nase a crypt cell homeostasis occurs, that is characterized by almost constant cell densities, pressure and cell distribution along the simulated time interval. In this scenario, parameters, as well as proliferation and differentiation rates, are set to match experimental measurements in a crypt exhibiting continuous, well-regulated cell renewal and differentiation, without disruptions that could lead to adenoma formation. These parameters are derived from medical observations and measurements reported in the literature (Vanleeuwen et al., 2006, 2009; Romset et al., 2022) and references therein.

We first observe that the theoretical cell densities \(C_A\), \(C_B\), and \(C_C\) in a homeostatic crypt, should be stationary solutions of the continuous model in equations (1), (4) and (5).

Specifically, as described in the section continuous model, stem cells (\(C_A\)) predominate at the crypt bottom, while semi-differentiated (\(C_B\)) and fully differentiated cells (\(C_C\)) are mainly located in the crypt middle and top, respectively. Moreover it is known that stem cells are not found in Normal Case over one third of crypt height, and fully differentiated cells are absent from the crypt bottom. Therefore, the following simple linear expressions for the initial cell densities \(C_A^0=C_A^0(x,y)\), \(C_B^0=C_B^0(x,y)\), \(C_C^0=C_C^0(x,y)\) are used in the Normal Case:

The corresponding graphs, equations (16)–(18), are given in Figure 3.

(a)

(b)

(c)

The associated discrete cell distribution of the VT model at time \(t=0\), where each Voronoi region is centered at an \(\Omega_h\) node, is represented in Figure 4(a). For completeness, Figure 4(b) shows the corresponding distribution at the final time T for the Abnormal Case 1, whereas Figures 4(c)–(d) present the initial and final distributions for Abnormal Case 2, respectively. These abnormal cases are described in the next sections.

In order to

(a)

(b)

(c)

(d)

In order to obtain the numerical pressure approximation \(p_h\) using equation (12), it is necessary to specify the tissue parameter \(\xi\), the cell velocity \(g_y\) at the crypt top, the proliferation rate \(\alpha_B\) and differentiation rate \(\beta_{AB}\).

Since stem cells at crypt bottom are nearly immobile (Vanleeuwen et al., 2009), we assume that the tissue parameter \(\xi\) is null at the very bottom of the crypt and increase upward, where cells are less resistant to pressure gradients due to higher permeability and lower viscosity. Then cells have an increased velocity going upwards along the crypt vertical axis as expected. Assuming that the tissue parameter remains constant at the crypt top, we adopt the following piecewise linear \(\xi\), which depends solely on the vertical coordinate \(y\)

\[\xi(y)= \begin{cases} \dfrac{3}{2h}y & \text{if } 0<y\leq \dfrac{2h}{3},\\[4pt] 1 & \text{if } y> \dfrac{2h}{3}. \end{cases} \]

It is known from (Potten et al., 1992) that cells move at most one cell position per hour near the crypt top orifice, then since our VT model each cell corresponds to \(3\times 3\) real cells, we have that the maximum velocity located on the crypt top is \(g_y=\frac{1}{3}\) of a cell position per hour.

The semi-differentiated proliferation rate \(\alpha_B\) in the Normal Case exhibits instead a bimodal behavior along the vertical coordinate \(y\), as observed experimentally (Almet et al., 2021).

Based on the experimental proliferation rate reported by (Potten et al., 1992), we derived the \(\alpha_B\) analytical expression in units of cells per hour, presented in equation (20)

\[\alpha_B(y)= \begin{cases} -a_1y^2 + a_2y + a_0 & \text{if } 0<y<\dfrac{h}{5},\\[4pt] -a_3\left(y-\dfrac{h}{5}\right)^2 + a_4 & \text{if } \dfrac{h}{5}\leq y \leq \dfrac{2h}{3},\\[4pt] a_5(y-h)^2 & \text{if } \dfrac{h}{5}\leq y < \dfrac{2h}{3}. \end{cases} \tag{19}\]

where \(a_0=\dfrac{1}{86}\), \(a_4=\dfrac{1}{29.9}\), \(a_1=\dfrac{25}{h^2}(a_4-a_0)\), \(a_2=\dfrac{10}{h}(a_4-a_0)\), \(a_5= \dfrac{9}{2h^2}a_0\), \(a_3=\dfrac{225}{49h^2}\left(a_4-a_5\dfrac{h^2}{9}\right)\). Details on how these coefficients \(a_i\) were determined are provided in (Masteritalo et al., 2024).

The remaining differentiation rates, \(\beta_{AB}\) and \(\beta_{BC}\), are also measured in units of cell number per hour. They are provided in Equations (21)–(22):.

The graph of \(\alpha_B\) is shown in Figure 5(a), while the differentiation rates \(\beta_{AB}\) and \(\beta_{BC}\) are shown in Figures 5(b) and 5(c), respectively.

(a)

(b)

(c)

Their expression arise from the experimental observation cited in (Nicolas et al., 2007) that semi-differentiated cells at the crypt top differentiate at most once per week, whereas stem and semi-differentiated cells do not differentiate at the crypt bottom and in the lower two-thirds of the crypt height, respectively.

The remaining parameters and time steps used in HCM are: simulation time \(T=1\), number of time steps \(N_T=100\) , continuous time step \(\Delta_T= 0.01\), diffusion coefficient \(D=1\), tissue stiffness coefficient \(\mu=50\), viscosity coefficient \(\eta=1\), negative threshold \(d_{min}=-1.5\).

Using the above initial densities in equations (16)–(18), proliferation, differentiation rates in equations (20)–(22), and problem parameters, we observe in Figures 6(a)–6(c), that cell densities do not change significantly compared to the initial ones, in fact only a regularization of the densities can be noted. Moreover, pressure solutions are also almost invariant in the time period \([0,T]\), see Figures 7(a)–7(b).

(a)

(b)

This is validating that we are modeling a homeostatic case with almost time constant velocity and cell densities.

Abnormal Case 1

In this scenario we induce abnormal proliferation by using a proliferation rate that generates five cells (per hour) more than the normal rate \(\alpha_B\) defined in equation (20). We also perturb the velocity \(g_y\) at the crypt top, specifically, it is increased fivefold compared to the Normal Case.

Thus, five cells per hour are now expelled from the crypt top orifice to the colon lumen instead of one as in the Normal Case. Starting by the normal initial conditions for densities given in equations (16)–(18), and considering the mentioned perturbations we have a large difference in the cell distribution at final time \(T\) as can be seen in Figures 8(a)–8(c).

(a)

(b)

(c)

Moreover, as it can be seen in Figure 8(b) and Figure 4(b) semi-differentiated cells fill almost completely the crypt.

Abnormal Case 2

In this second abnormal scenario different initial conditions are used: stem cells are initially located in the middle of the crypt (\(C_A=1\) in \([\frac{h}{2}-1,\frac{h}{2}+1]\)), see Figures 9(a)–9(c).

(a)

(b)

(c)

This is a significant modification with respect to the cell distribution in a healthy normal crypt where only semi-differentiated cells can reside at the crypt middle. In this scenario we want to simulate what is observed in lesions of the colon that present abnormal stem cells located in the middle of the crypts, because it is known that this can lead to or be responsible for adenoma formation. In fact as described in (Romanazzi & Settanni, 2022) and references therein there is a correlation between adenoma formation and stem cell location far from the bottom that are characterized by high Wnt.

Applying our computational method at time \(t=T\) we get that stem cells continue to be present in the crypt middle, see densities plots in Figures 10(a)–10(c). This demonstrates that a such abnormality is not removed along the time and then it can persist and potentially lead to lesions and adenoma formation.

(a)

(b)

(c)

The associated VT cell distribution in this Abnormal Case at initial time and final time are given in Figures 4 (c)-(d).

Computational Tools

All simulations were implemented in Python (Python et al., 2024). The continuous system was solved using the FEniCS (Fenics et al., 2015) package, which enables the formulation and solution of partial differential equations using the finite element method. The discrete model was also implemented in Python, using libraries for geometric processing and simulation of cell interactions. Results were visualized using the Matplotlib Python library.

Conclusions

In this work we have presented a hybrid cellular model for simulating cell dynamics in colonic crypts. This framework integrates a cell-centered Voronoi Tessellation (VT) model with a continuous PDE system to track cell densities and position in a cylindrical crypt. It reproduces the normal homeostasis state and enables the testing of different abnormal scenarios. In particular we have shown in an Abnormal Case that stem cells initially located at the crypt middle, that is uncommon in homeostasis regime, continue to reside there along the simulated time interval. Thus such an initial abnormal condition perturbs the cellular dynamics permanently. This result is important because it validates the biological observations and theoretical hypotheses that modifying cell proliferation or stem cell locations in a crypt can lead to abnormal crypts with lesions.

Our hybrid framework allows to describe simultaneously cell proliferation, convective cell velocity due to mitotic activity directed to the crypt top and cell differentiation which occur over longer timescales (days), while cell-scale dynamics occur over shorter timescales (hours). Numerical simulations have been implemented by a finite element method which has been coupled with the VT cell model solved in time by a Forward Euler method. This methodology is computationally efficient with respect a fully cellular discrete model, such as that presented in (Vanleeuwen et al., 2009), which solves many differential ordinary equations presenting parameters linked with a provided Wnt gradient.

As future work, we plan to enhance our hybrid model by introducing a Wnt gradient which can perturb the model parameters and rates at cellular scale rather than just on the continuous scale as done in this work.

Acknowledgments

Authors thanks the financial support received by "Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Código de Financiamento 001",  "Fundação de Amparo à Pesquisa e ao Desenvolvimento Científico e Tecnológico do Maranhão (FAPEMA) - Código de Financiamento BM-06193/22", and "Programa de Pos-Graduação de Matemática Aplicada do IMECC/UNICAMP".

Author Contributions

I. J. L. Sousa participated in conceptualization, programming, data curation, formal analysis, investigation, methodology, writing – original draft, and visualization. G. Romanazzi participated in conceptualization, data curation, formal analysis, investigation, methodology, writing – original draft, and visualization.

Conflicts of interest

The authors declare no conflict of interest.

References

Almet, A. A., Byrne, H. M., Maini, P.K. & Moulton, D. E. (2021). The role of mechanics in the growth and homeostasis of the intestinal crypt. Biomechanics and Modeling in Mechanobiology, 585-608. https://doi.org/10.1007/s10237-020-01402-8

Alns, Martin S., Blechta, Jan, Hake, Johan, Johansson, August, Kehlet, Benjamin, Logg, Anders, Richardson, Chris N., Ring, Johannes, Rognes, Marie E. & Wells, Garth N. (2015). The FEniCS Project Version 1.5. Archive of Numerical Software 3(100), 9-23. https://doi.org/10.11588/ans.2015.100.20553

Asarian, L., Gloy, V., & Geary, N. (2012). Homeostasis. Academic Press, 324-333. https://doi.org/10.1016/B978-0-12-375000-6.00191-9

Baker, A. M., Cereser, B., Melton, S., Fletcher, A. G., Rodriguez-Justo, M., Tadrous, P. J., Humphries, A., Elia, G., McDonald, S. A., Wright, N. A., Simons, B. D., Jansen, M. & Graham, T. A. (2014). Quantification of crypt and stem cell evolution in the normal and neoplastic human colon. Cell reports 8(4), 940-947. https://doi.org/10.1016/j.celrep.2014.07.019

Bernstein, C., Facista, A., Nguyen, H., Zaitlin, B., Hassounah, N., Loustaunau, C., Payne, C. M., Banerjee, B., Goldschmid, S., Tsikitis, V. L., Krouse, R. & Bernstein, H. (2010). Cancer and age related colonic crypt deficiencies in cytochrome c oxidase I. World Journal of Gastrointestinal Oncology 2(12), 429-442. https://doi.org/10.4251/wjgo.v2.i12.429

Chandrasinghe, P. C. (2018). Basics in molecular evolution of colorectal cancer and their implications for the surgeon: is it a 'big-bang' or a 'survival of the toughest'?. The Sri Lanka Journal of Surgery 36(2), 18-21. https://doi.org/10.4038/sljs.v36i2.8511

Di Garbo, A., Johnston, M. D., Chapman, S. J. & Maini, P. K. (2010). Variable renewal rate and growth properties of cell populations in colon crypts. Physical Review E 81(6), 061909. https://doi.org/10.1103/PhysRevE.81.061909

Drasdo, Dirk & Loeffler, Markus (2001). Individual-based models to growth and folding in one-layered tissues: intestinal crypts and early development. Nonlinear Analysis: Theory, Methods & Applications 47(1), 245-256. https://doi.org/10.1016/S0362-546X(01)00173-0

Figueiredo, I. N., Leal, C., Romanazzi, G. & Engquist, B. (2019). Biomathematical model for simulating abnormal orifice patterns in colonic crypts. Mathematical Biosciences 315, 108221. https://doi.org/10.1016/j.mbs.2019.108221

Greenspan, H. P. (1976). On the growth and stability of cell cultures and solid tumors. Journal of Theoretical Biology 56(1), 229-242. https://doi.org/10.1016/S0022-5193(76)80054-9

Guebel, D. V. & Torres, N. V. (2008). A computer model of oxygen dynamics in human colon mucosa: Implications in normal physiology and early tumor development. Journal of Theoretical Biology 250(3), 389-409. https://doi.org/10.1016/j.jtbi.2007.09.035

Kershaw, S. K., Byrne, H. M., Gavaghan, D. J. & Osborne, J. M. (2013). Colorectal cancer through simulation and experiment. IET Systems Biology 7(3), 57-73. https://doi.org/10.1049/iet-syb.2012.0019

Meineke, Frank A., Potten, Christopher S. & Loeffler, Martin (2001). Cell migration and organization in the intestinal crypt using a lattice-free model. Cell Proliferation https://doi.org/10.1046/j.0960-7722.2001.00216.x

Murray, P. J., Walter, A., Fletcher, A. G., Edwards, C. M., Tindall, M. J. & Maini, P. K. (2011). Comparing a discrete and continuum model of the intestinal crypt. Physical biology 8(2), 026011.

Nicolas, P, Kim, K. M,, Shibata, D. & Tavaré, S. (2007). The stem cell population of the human colon crypt: Analysis via methylation patterns. PLoS Computational Biology 3(3), e28. https://doi.org/10.1371/journal.pcbi.0030028

Oliveira, E. P. & Romanazzi, G. (2022). Modeling and computer simulation of viscoelastic crypt deformation. Trends in Computational and Applied Mathematics 23(1), 193-211. https://doi.org/10.5540/tcam.2022.023.01.00193

Osborne, J. M. (2015). A multiscale model of colorectal cancer using the cellular Potts framework. Cancer Informatics 14, 83-93. https://doi.org/10.4137/CIN.S193

Potten, C. S., Kellett, M., Roberts, S. A., Rew, D. A. & Wilson, G. D. (1992). Measurement of in vivo proliferation in human colorectal mucosa using bromodeoxyuridine. Gut 33(1), 71-78. https://doi.org/10.1136/gut.33.1.71

Python Software Foundation (2024). Python Language Reference, Version 3.2. https://www.python.org/download/releases/3.2/

Romanazzi, G. & Settanni, G. (2022). Mathematical Model for Simulation of Morphological Changes associated to Crypt Fission in the Colon. Discrete and Continuous Dynamical Systems - Series S 15(12), 3781-3805. https://doi.org/10.3934/dcdss.2022055

Sousa, I. J. L. (2024). Modelos Celulares nas Criptas do Cólon: Análise e Implementação Numérica [Dissertação de Mestrado, Universidade Estadual de Campinas]. Repositório Institucional. https://doi.org/10.47749/T/UNICAMP.2024.1408687

van Leeuwen, I. M. M., Byrne, H. M., Jensen, O. E. & King, J. R. (2006). Crypt dynamics and colorectal cancer: Advances in mathematical modelling. Cell Proliferation 39(3), 157-181. https://doi.org/10.1111/j.1365-2184.2006.00378.x

van Leeuwen, I. M. M., Mirams, G.R., Walter, A., Fletcher, A., Murray, P., Osborne, J. M., Varma, S., Young, S. J., Cooper, J., Doyle, B.,et al. (2009). An integrative computational model for intestinal tissue renewal. Cell Proliferation, 617-636. https://doi.org/10.1111/j.1365-2184.2009.00627.x