Modeling and Computer Simulation of Viscoelastic Crypt Deformation

ABSTRACT

Colorectal cancer morphogenesis begins at the cellular level from cell mutations in the intestinal epithelium cavities called crypts. These mutations lead to a pressure difference in the epithelium crypt walls, which can cause deformation and generate visible abnormalities in the epithelium. The geometrical modeling of these crypts and the mathematical modeling of the biomechanical process that leads to deformations can be simulated by using a Finite Element Method. The method solves numerically the system of partial differential equations (PDEs) that governs this phenomenon and permits to estimate the deformations of the crypt walls. In this work we simulate the crypt deformation when the cell mutations appear in several regions of the crypt epithelium.

Keywords:
colorectal cancer; crypts; finite element method; computer simulation

1 INTRODUCTION

In the process of the development of colorectal cancer, called carcinogenesis, cell mutations can arise in colon cavities, called crypts, containing transit and stem cells ³⁰. Mutant cells usually show higher proliferation rate than normal cells, which causes a larger pressure on the crypts walls than that generated by normal cells. This leads to a deformation of the crypt walls and its top orifice. If mutant cells fill adjacent crypts, a set of abnormal and deformed crypts appear. Thereafter there is a growth of the epithelium tissue caused by these abnormal crypts inside the colon that evolves to an adenoma. Adenomas in the colon are characterized in early stages by the formation of larger and deformed orifices of the epithelium crypts, which have a displacement in the intestinal lumen ^{14), (30}.

Regarding the cell and tissue growth modeling, we present in Section 2 a spatial continuous model based on PDEs, similar to that used in ¹⁷, that determines the density and adhesion cell-cell pressure of fully differentiated and proliferative transit cells in a single crypt.

The exact morphology resulting from an abnormal proliferation of cells at some locations along the crypt axis is not known. However this leads to a deformation that can proceed to a budding, a orifice enlargement or an abnormal crypt fission ^{1), (15}. We follow a model similar of that introduced in ¹⁷, but instead of using a biomechanical model considering only the crypt orifice deformation, we consider the entire crypt walls deformation which permits to simulate the colon epithelium tissue growth generated by abnormal crypt cells located in any position in relation to the crypt base. This allows to be more effective in understanding how orifice enlargement and deformation observed in usual colonoscopy images can be originated by abnormal cell proliferation along the usual colonoscopy camera hidden crypt axis. Understanding such origin of crypt morphology permits to associate different orifice deformation shapes to particular locations of the abnormal cells along the crypt axis. It is then feasible to apply a specific target therapy to these locations to block the adenoma formation since cells at different crypt axis quotes have different level of Wnt/-catenin/TCF4 signaling pathway and APC protein. Wnt pathway is for instance more active in the lower part of the crypt and decreases moving upwards up to be null at the crypt top, whereas APC protein has an opposite gradient along the crypt axis ^{17), (23}. Such a model based on longitudinal crypt section can allow a model calibration based on real data of a three dimensional crypt as that obtained by the scanning electron microscope SEM as that used in ^{4), (25}.

In literature, there are many spatial models describing the cell spatial location, also called cellbased models, including in-lattice or grid-models ^23),(30),(31 and off-lattice or lattice-free models ^{2), (9), (19}. The choice to preferring a spatial continuum model with respect to a cell-based model is in the possibility of extend its application to simulate cell dynamics in millions of crypts that can be done using a multiscale strategy, as done in ¹⁶.

For the tissue growth modeling, as highlighted in reviews ^3),(11),(24, many possible models based on cell dynamics in the intestinal crypts have been considered, among them: continuum mechanics models ^{15), (27}; compartmental models describing the transition between cell types, regardless of their position ^{7), (8}; and nonspatial stochastic models based on the use of epidemiological data to predict cancer formation or some genetic instability ^{12), (13}. Our model determines the displacement and deformation of the epithelium tissue along the crypt based on a continuum equilibrium relation, this allows us to use directly the pressure gradient between cells as a force acting on the tissue.

2 MODELING CRYPT CELL DYNAMICS AND DEFORMATION

In this Section we present the geometrical domain used to represent a crypt and two half crypts at its left and right in a cross section of the colon epithelium. Then we introduce the PDE system that describes the cell dynamics in this domain, see Subsection 2.2. Later in the Subsection 2.3 we define a viscoelastic model responsible for the tissue deformation along the crypt, which coupled with the cell dynamics model allows to deform the geometrical domain of the crypt and to observe, in some abnormal cases, tissue growth as it is shown in Section 4.

2.1 Geometrical Domain

The colon epithelium is formed by a periodic distribution of crypts that in three dimension can be locally represented as in Figure 1. We study and analyze the deformation occurring in a small region Ω of the colon epithelium during a time interval [0, T]. Its external crypt subdomain Ω_E is formed by the cross section of the external epithelium tissue, which has one crypt in its center and two walls of two adjacent crypts at its left and right, see Figure 2.

Figure 1:
Small portion of the epithelium tissue formed by 9 crypts.

Figure 2:
Domain Ω_E formed by a cross section of a single crypt in the center with two adjacent crypt walls at left and right.

The domain Ω has also a tissue region Ω_C that is formed by the connective tissue of the colon epithelium between the crypts. Thus we model the deformations and cell dynamics in $\Omega ={\Omega }_E\cup {\Omega }_C$, see Figure 3. In Ω_E , we use a diffusive-convective model for the dynamics of proliferative (also called transit) cells and fully differentiated cells. Both Ω_E and Ω_C have elastic structures that can move during the time interval [0, T]. Figure 4 illustrates the boundary domain subdivisions and their notations. Region Ω_E has an external boundary formed by Γ_r that contains its upper circular borders, Γ₃ is the vertical boundary of the Ω_E bottom and Γ_top is the horizontal boundary of the Ω_E top. We note that normal mature (fully differentiated) cells are released in the colon lumen through this last boundary Γ_top³⁰. The remaining external Ω_E boundary is denoted by Γ₁. In the interface between Ω_E and Ω_C we localize Γ_bot that represents the circular bottom boundary, and Γ₂ that is the remaining interface boundary. Furthermore we denoted by Γ₄ the lower Ω_C boundary where the tissue is supposed to be fixed.

Figure 3:
Full Problem domain Ω with the external crypt tissue Ω_E and the internal connectivity tissue Ω_C .

Figure 4:
Subdivisions of the boundary.

Since we deform such domains in time, we denote by $\Omega \left(t\right)$ the geometrical domain of the crypt at time $t\in [0,T]$ that is deformed in relation to the original domain Ω represented in the figures above. Similarly we denote by ${\Omega }_E\left(t\right)$ and ${\Omega }_C\left(t\right)$ the deformed external and connective tissue domain at time t.

2.2 PDE model for cell dynamics

We suppose, as done in ^{17), (18}, that a crypt contains two families of cells: proliferative transit cells with density $N_1=N_1(\mathbf{x},t)$ and fully differentiated living cells with density $N_2=N_2(\mathbf{x},t)$ where $\mathbf{x}=(x,y)\in {\Omega }_E\left(t\right)$ and $t\in [0,T]$. It is also supposed that N ₁ and N ₂ satisfy the overall density condition $N_1+N_2=1$ and that cell families have the same diffusion coefficient D as it is asserted in ²⁸. According to the model used in ^17),(18, N ₁ and N ₂ satisfy in ${\Omega }_E\left(t\right)$ for $t\in (0,T]$ the following PDE system

$\left\{\begin{array}{c}\begin{array}{l}\frac{{\displaystyle \partial N_1}}{{\displaystyle \partial t}}-\nabla ·(\nabla p{N}_1)=\nabla ·(D\nabla N_1)+\alpha N_1-\beta N_1\\ \frac{{\displaystyle \partial N_2}}{{\displaystyle \partial t}}-\nabla ·(\nabla p{N}_2)=\nabla ·(D\nabla N_2)+\beta N_1,\\ N_1+N_2=1\end{array}\end{array}\right.,$ (2.1)

where α is the proliferation rate of N ₁, β is the proliferation rate of N ₂, p is the cell-cell adhesion pressure responsible for the upwards cell transition along the crypt axis. Since $N_1+N_2=1$, if we sum the first two equations of (2.1) and take $N=N_1$, the temporal and spatial derivatives cancel. Thus we obtain that $N=N(\mathbf{x},t)$ and $p=p(\mathbf{x},t)$ satisfy in ${\Omega }_E\left(t\right)$ for $t\in (0,T]$:

$\left\{\begin{array}{l}\frac{\partial N}{\partial t}-\nabla .(\nabla pN)=\nabla .(D\nabla N)+\alpha N-\beta N\\ -\Delta p=\alpha N.\end{array}\right.$ (2.2)

The boundary conditions for the pressure are

$\left\{\begin{array}{ll}p& =0,\text{if}\mathbf{x}\in {\Gamma }_{top}\\ \eta p+(1-\eta )\frac{\partial p}{\partial n}& =0,\text{if}\mathbf{x}\in {\Gamma }_r\\ \frac{\partial p}{\partial n}& =0,\text{if}\mathbf{x}\in {\Gamma }_3\cup {\Gamma }_2\cup {\Gamma }_1\cup {\Gamma }_{bot}\end{array}\right.$ (2.3)

In the second equation of (2.3), η is the modulus of y-coordinate of outward normal vector in Γ_r . This is a Robin boundary condition that allows the transition between the Dirichlet and Neumann conditions in Γ₁ and Γ_top .

For the density, the boundary conditions are

$\left\{\begin{array}{ll}N& =1,\text{if}\mathbf{x}\in {\Gamma }_3\\ D\frac{\partial N}{\partial n}+N\frac{\partial p}{\partial n}& =0,\text{if}\mathbf{x}\in {\Gamma }_{top}\cup {\Gamma }_r\cup {\Gamma }_{bot}\\ \frac{\partial N}{\partial n}& =0,\text{if}\mathbf{x}\in {\Gamma }_1\cup {\Gamma }_2\end{array}\right.$ (2.4)

The proliferation rate α is higher at the crypt bottom and becomes null at the top, normally in the upper third of the crypt height where most of the fully differentiated cells, that are not proliferative, reside ³⁰. On the other hand, the transformation rate β behaves in the opposite way. It means that it is higher at the crypt top and null at the bottom, thus we use as in ^{17), (18} the rates

$\alpha \left(y\right)=\left\{\begin{array}{ll}\tau (y-\frac{2}{3}h{)}^2,& 0\le y\le \frac{2}{3}h\\ 0,& \mathrm{elsewhere}\end{array}\right.,\beta \left(y\right)=\left\{\begin{array}{ll}0,& 0\le y\le \frac{2}{3}h\\ \tau (y-\frac{2}{3}h{)}^2,& \mathrm{elsewhere}\end{array}\right.$ (2.5)

where h denotes the crypt height and $\tau >0$ is a constant that is determined by using some biological information of the analyzed crypt such as the cell size and the cell convective velocity along the crypt axis.

The full model (2.2) with the boundary conditions (2.3)-(2.4) guarantees that transit cells diffuse in the crypt, proliferate with rate α and differentiate with rate β . Moreover these cells have a convective upward velocity $v=-\nabla p$. This velocity increases along the vertical crypt axis up to reach the crypt top, where fully differentiated cells are released in the lumen, as it has been observed experimentally ³⁰. We note in fact that the boundary conditions for the pressure and density in Γ permit cells to move upwards, and to be released in the colon lumen across the Γ_top boundary.

2.3 The Viscoelastic Model

As discussed in ¹⁷, the crypt domain Ω can be considered as a viscoelastic structure. Then it can have a displacement $\mathbf{u}(\mathbf{x},t)=\left(u_1\right(\mathbf{x},t),u_2(\mathbf{x},t\left)\right)$ in relation to its original position, where u ₁ and u ₂ are the displacement along the x−axis and y−axis of point x at time t, respectively. Since connective tissue domain base is normally considered as fixed we allow only Γ₁, Γ₂, Γ_bot , Γ_top , Γ_r to be free boundaries whereas Γ₃, Γ₄ remain fixed. The gradient of pressure differences acts as a force on the external crypt boundary ${\Gamma }_1\cup {\Gamma }_r\cup {\Gamma }_{top}$, which causes a deformation of Ω, with a displacement u. The stress tensor $\sigma \left(\mathbf{u}\right)=\left({\sigma }_{ij}\right(\mathbf{u}){)}_{i,j=1,2}$ has elastic ${\sigma }_{ij}^{el}$ and viscous ${\sigma }_{ij}^{vs}$ components, so that ${\sigma }_{ij}={\sigma }_{ij}^{el}+{\sigma }_{ij}^{vs}$. Since Ω_E contains only cells it relaxes more linearly than Ω_C and so it can be represented as follows with an elastic structure, see ¹⁷:

${\sigma }_{ij}^E\left(\mathbf{u}\right)={\lambda }_E\nabla ·\mathbf{u}{\delta }_{ij}+{\mu }_E\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right).$ (2.6)

In (2.6) λ _E and µ _E are the Lamé coefficients, which describe mechanical properties of elastic materials ²⁹. On the other hand, in Ω_C , we have a viscoelastic structure that has elastic and viscous components. We use a Standard Linear Solid (SLS) model for modeling the viscoelastic tissue domain Ω_C . This model is largely used in the representation of biological tissues ²⁶. Then, using the Einstein summation convention, we have

${\sigma }_{ij}^C\left(\mathbf{u}\right)=D_{ijkl}(\mathbf{x},0){ϵ}_{kl}\left(\mathbf{u}\right)(\mathbf{x},t)-{\int }_0^t\frac{\partial D_{ijkl}}{\partial s}(\mathbf{x},t-s){ϵ}_{kl}\left(\mathbf{u}\right)(\mathbf{x},s)ds,$ (2.7)

where

${ϵ}_{ij}\left(\mathbf{u}\right)=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right),$ (2.8)

$D(\mathbf{x},t)=(D_{ijkl}{)}_{ijkl=1,2}$ with

$$\begin{gathered} D_{11kl}(\mathbf{x},t)=\left(\begin{array}{cc}\lambda (\mathbf{x},t)+2\mu (\mathbf{x},t)& 0\\ 0& \lambda (\mathbf{x},t)\end{array}\right),D_{12kl}(\mathbf{x},t)=\left(\begin{array}{cc}\mu (\mathbf{x},t)& 0\\ 0& \mu (\mathbf{x},t)\end{array}\right), \\ {D}_{21kl}(\mathbf{x},t)=\left(\begin{array}{cc}0& \mu (\mathbf{x},t)\\ \mu (\mathbf{x},t)& 0\end{array}\right),D_{22kl}(\mathbf{x},t)=\left(\begin{array}{cc}0& \lambda (\mathbf{x},t)+2\mu (\mathbf{x},t)\\ \lambda (\mathbf{x},t)& 0\end{array}\right) \end{gathered}$$ (2.9)

where $\lambda (x,t)={\lambda }_E\left(1+e^{-\frac{t}{{\tau }_0}}\right),\mu (x,t)={\mu }_E\left(1+e^{-\frac{t}{{\tau }_0}}\right)$ with ${\tau }_0=5.43s$, see ^{10), (17}.

At the equilibrium, the crypt structure displacements satisfy for $i=1,2$ for all $t\in (0,T]$

$\left\{\begin{array}{ccc}-\nabla ·{\sigma }_i^E\left(\mathbf{u}\right)& =& f_i,\text{in}{\Omega }_E\left(t\right)\\ -\nabla ·{\sigma }_i^C\left(\mathbf{u}\right)& =& 0,\text{in}{\Omega }_C\left(t\right)\\ \mathbf{u}& =& 0,\text{on}{\Gamma }_3\cup {\Gamma }_4\\ {\sigma }_i\left(\mathbf{u}\right)·\mathbf{n}& =& 0,\text{on}\partial {\Gamma }_1\\ u_i(\mathbf{x},0)& =& \frac{\partial u_i}{\partial t}(\mathbf{x},0)=0,\text{in}\Omega \left(t\right)\end{array}\right.$ (2.10)

where ${\sigma }_i\left(\mathbf{u}\right)=\left({\sigma }_{i1}\right(\mathbf{u}),{\sigma }_{i2}(\mathbf{u}\left)\right)$ and $\mathbf{n}=(n_1,n_2)$ is the outward normal vector in Γ₁. Furthermore, we have $\mathbf{f}=-\gamma \nabla (p-p^{*})$, with $\mathbf{f}=(f_1,f_2),p^{*}$ being the pressure exerted between the cells in the normal case and p the exerted pressure in the abnormal case, that causes the crypt deformation, and γ is an adimensional term related with the current geometry of the tissue domain and cell dimensions.

3 METHODOLOGY

We simulate the crypt deformations and the associated tissue growth during the first stages of colorectal cancer in which abnormalities occur in a single crypt as shown in Section 4. We use a classic Finite Element Method with piecewise linear basis to solve the differential problem in different time instants and we advance in time with a fixed time interval ∆t by using a Backward Euler’s method. Once the crypt deforms as the process goes forward, the domain mesh changes in each time iteration. The computational code is implemented using the FreeFem++ software ²², that solves numerically the density N, pressure p, and displacement u by using a piecewise linear finite element basis of the domain $\Omega \left(t\right)$ at each $t\in (0,T]$.

4 RESULTS AND DISCUSSIONS

In order to simulate different crypt deformations and the associated tissue growth due to abnormal crypt cell proliferation, we need to assert the normal case and how it differs with respect the abnormal cases.

4.1 Stable Solution for N

In order to have a stable solution of system (2.2) with boundary conditions we find an expression for the transformation rate β, for a given α and a stable solution N. We suppose that

$\alpha \left(y\right)=\left\{\begin{array}{cc}\tau {\left(y-\frac{2h}{3}\right)}^2& \text{if}0\le y\le \frac{2h}{3}\\ 0& \text{if}y>\frac{2h}{3}\end{array}\right.$ (4.1)

in which τ is the relaxation constant, and the stable solution is

$N=N(x,y,t)=\left\{\begin{array}{cc}1& 0\le y<y_b\\ N_s\left(y\right)& y_b\le y\le y_t\\ 0& \begin{array}{l}{y}_t<y\le h\end{array}\end{array}\right.$ (4.2)

Where $N_s\left(y\right)=\frac{y-y_t}{y_b-y_t}$, in which yt is the minimum (vertical) coordinate y of the top region and y _b is the maximum (vertical) coordinate of the bottom region. Since the solution N is stable, we have by hypothesis that $\frac{\partial N}{\partial t}=0$. Then, from the first equation of (2.2), we have

$\Delta pN+\nabla p·\nabla N+\nabla ·(D\nabla N)+(\alpha -\beta )N=0.$ (4.3)

But, from the second equation, $-\Delta p=\alpha N$, with D constant, we obtain in $[y_b,y_t]$

$\beta N_s\left(y\right)=\alpha \left(y\right)[1-N_s(y\left)\right]N_s\left(y\right)+\nabla p·(0,N{\text{'}}_s(y\left)\right)+DN\text{'}{\text{'}}_s\left(y\right).$ (4.4)

Moreover since $N{\text{'}}_s\left(y\right)=\frac{1}{y_b-y_t}$ and $N\text{'}{\text{'}}_s\left(y\right)=0$, we obtain

$\beta N_s\left(y\right)=\alpha \left(y\right)[1-N_s(y\left)\right]N_s\left(y\right)+\frac{\partial p}{\partial y}·\frac{1}{y_b-y_t}.$ (4.5)

Since $\frac{\partial p}{\partial x}=0$ on ${\Gamma }_1\cup {\Gamma }_2$ and $p=0$ on Γ_top , we assume that $p=p\left(y\right)$ in order to get a simpler solution. Then, recalling that $\Delta p=-\alpha N_s\left(y\right)$, we have

$p\text{'}\text{'}\left(y\right)=-\alpha \left(y\right)N_s\left(y\right).$ (4.6)

Integrating both sides of the equation we obtain

$p\text{'}\left(y\right)-p\text{'}\left(y_t\right)={\int }_y^{y_t}\alpha \left(\xi \right)N_s\left(\xi \right)d\xi ,$ (4.7)

for all $y\in [y,y_t]$.

We know that, at the crypt top, $v_{top}=-\frac{\partial p}{\partial y}\left(y_t\right)=0.85$ position/hour = 0.85 h _cell µm/h³⁰. Moreover, since a crypt has a height 433µm on average, and cells have height 5.9µm²⁰, then for a crypt in the model with height h _crypt , we have

$h_{cell}=\frac{5.9h_{crypt}}{433}.$ (4.8)

Thus $v_{top}=-\frac{\partial p}{\partial y}\left(y_t\right)=0.85·\frac{5.9h_{cell}}{433}=0.01158h_{crypt}$. Then

$\frac{\partial p}{\partial y}=-v_{top}+{\int }_y^{y_t}\alpha \left(\xi \right)\left[\frac{\xi -y_t}{y_b-y_t}\right]d\xi .$ (4.9)

If $0<y\le \frac{2h}{3}$, we have

${\int }_y^{y_t}\alpha \left(\xi \right)\left[\frac{\xi -y_t}{y_b-y_t}\right]d\xi ={\int }_y^{\frac{2h}{3}}\tau {\left(\xi -\frac{2h}{3}\right)}^2\left(\frac{\xi -y_t}{y_b-y_t}\right)d\xi .$ (4.10)

On the other hand, if $y>\frac{2h}{3}$, we have

${\int }_y^{y_{top}}\alpha \left(\xi \right)\left[\frac{\xi -y_y}{y_b-y_t}\right]d\xi =0.$ (4.11)

Thus, we get

$\frac{\partial p}{\partial y}=-v_{top}+{\int }_y^{\frac{2h}{3}}\tau {\left(\xi -\frac{2h}{3}\right)}^2\left(\frac{\xi -y_t}{y_b-y_t}\right)d\xi ,$ (4.12)

if $0<y\le \frac{2h}{3}$ and

$\frac{\partial p}{\partial y}=-v_{top},$ (4.13)

if $y>\frac{2h}{3}$. At last, if we define β only for $y_b\le y\le y_{lim}$ with $y_{lim}<y_t$ and recall that $N_s\left(y_t\right)=0$, we obtain

$\beta \left(y\right)=\alpha \left(y\right)[1-N_s(y\left)\right]+\frac{\partial p}{\partial y}·\frac{1}{y-y_t},$ (4.14)

where $N_s\left(y\right)$ is a solution of the problem

$\frac{\partial N}{\partial t}-\nabla ·(\nabla pN)=\nabla ·(D\nabla N)+\alpha \left(y\right)N-\beta \left(y\right)N.$

Figures 5 and 6 show the stable solution for density and the pressure solution of the equation $-\Delta p=\alpha N_s$.

Figure 5:
Stable solution N _s .

Figure 6:
Pressure p associated with the stable density solution N _s .

4.2 Numerical Simulations

We tested many abnormal examples each leading to a different cell dynamics and crypt deformation. In each test the abnormality, characterized by a region with an unexpected high proliferation rate α, is located in a different region. We use the following parameter values in each simulation: crypt height $h=433\mu m,D=1$ (µm)²/hour, $\gamma =1·{10}^{-6},\Delta t=1·{10}^{-6}$ hour, $\tau =1·{10}^{-5}$ hour⁻¹(µm)⁻² and $H_{max}=6.67369\mu m$µm, with H _max being the maximum diameter of the mesh triangles. Moreover we use the boundary conditions for the pressure and density (2.3)-(2.4) and the initial condition $N(x,y,0)$ given in (4.2) with $y_t=h=433\mu m$ and $y_b=37\mu m$.

In all simulations, meshes are automatically generated by the software FreeFem++. The finite element mesh used at the beginning of the simulations, at time t=0, have 17174 triangles and 9307 vertices. After each time step, the mesh is updated based on the calculated displacements using the “movemesh” routine. If this update causes overlapping of the elements, FreeFem++ adapts the mesh using the “adaptmesh” routine based on the BAMG mesh adaptive method ²¹.

4.2.1 Test I

In the first test the abnormality is centered in the lower region of the crypt, as shown by the black circle in Figure 7.

Figure 7:
Initial mesh and the abnormality location marked by a black circle located in the crypt bottom.

After 15 time steps, the solutions u ₁ and u ₂ calculated on the updated mesh are shown in 8 and 9. Moreover, solutions N and $p-p^{*}$ are plotted respectively in Figures 10 and 11.

Figure 8:
Displacement u ₁ after 15 time steps.

Figure 9:
Displacement u₂ after 15 time steps.

Figure 10:
Cell density N after 15 time steps.

Figure 11:
Pressure difference $p-p^{*}$ after 15 time steps.

We observe that the displacements along the x-axis are symmetric since they are the same in the right and left crypt walls, according to the expectations. It results then in a symmetric deformation problem. Moreover, we note a relevant opening of the crypt orifice along the time, that characterizes microadenomas such as the Aberrant Crypt Foci¹, see ^{16), (17), (30}. Besides that, we have a large pressure difference in the region of the abnormality. This pressure difference decreases as we ascend the crypt. Regarding cell density N, we note that it is equal to 1 at the crypt base and decreases its value moving upwards in the direction of the crypt top where it has a null value.

4.2.2 Test II

In the second simulation the abnormality is located in the left branch (wall) of the crypt, in a centered position in relation to the y-axis, as it can be seen by the black circle in Figure 12.

Figure 12:
Initial mesh and the abnormality location in the middle of the crypt axis.

After 20 time steps, the displacements u ₁ and u ₂ computed on the updated mesh are shown in Figures 13 and 14. N and $p-p^{*}$ are respectively shown in Figures 15 and 16.

Figure 13:
Displacement u ₁ after 20 time steps.

Figure 14:
Displacement u ₂ after 20 time steps.

Figure 15:
Density N after 20 time steps.

Figure 16:
Pressure difference $p-p^{*}$ after 20 time steps.

In this simulation we observe that there is a larger displacement in the left crypt wall in relation to the portion on the right. This behavior agrees with what is expected because the abnormality is on the left side and, therefore, exerts larger influence in this region than on the right. Once again the pressure difference is large close to the abnormality and this difference decreases as we move away from this region. In this case, the density behavior is similar to that observed in test 1: it varies from one to zero, from the bottom to the top.

4.2.3 Test III

In the third situation, the abnormality is on the left branch of the crypt, in a higher position, as shown by Figure 17.

Figure 17:
Initial mesh with abnormality location in the top crypt left corner.

After 20 time steps, the solutions to u ₁ and u ₂ calculated on the updated mesh are shown in Figures 18 and 19. Besides that, solutions to N and $p-p^{*}$ are respectively shown in Figures 20 and 21.

Figure 18:
Solution to u ₁ after 20 time steps.

Figure 19:
Solution to u ₂ after 20 time steps.

Figure 20:
Solution to N after 20 time steps.

Figure 21:
Solution to $p-p^{*}$ after 20 time steps.

In this simulation, the region on the left of the crypt moves, forming an opening in this direction. The right-hand portion of the crypt has a very small displacement. This occurs because the difference in pressure in the right portion is very close to zero, since the abnormality is found in the left portion. Again, cell density decreases along the crypt vertical axis, from one at the bottom to zero at the top.

4.2.4 Test IV

In the fourth simulation, the abnormality is centered in the upper right portion of the crypt, as shown by the black circle in Figure 22.

Figure 22:
Initial mesh with abnormality location.

After 10 time steps, the solutions to u ₁ and u ₂ calculated on the updated mesh are shown in Figures 23 e 24. Moreover, the solutions N and $p-p^{*}$ are respectively shown in Figures 25 and 26.

Figure 23:
Horizontal displacement u ₁ obtained after 10 time steps.

Figure 24:
Vertical displacement u ₂ obtained after 10 time steps.

Figure 25:
Solution to N after 10 time steps.

Figure 26:
Solution to $p-p^{*}$ after 10 time steps.

In this test we observe an elongation of the left domain region of the epithelium due to the presence of upward-oriented forces in this region. Such forces come from the pressure difference in the abnormality region. The displacements in y-axis, shown by u ₂, are larger than those in x-axis. Moreover, the pressure difference is larger in the upper left portion of the crypt, being practically null in the entire right branch. Cell density presents a behavior similar to that of the other tests, decreasing from the crypt bottom to the crypt top.

5 CONCLUSIONS

The process of colorectal cancer development arises from the presence of abnormal cells in the intestinal crypts. This leads to a pressure difference, which generates tissue displacements. Using a PDE model, it was possible to numerically simulate the tissue deformation over time by using the Finite Element Method. This is a problem in which the domain is updated in each time iteration.

We analyzed how different positions of the abnormality in the tissue affects the process evolution. The displacements occurs in different ways depending on the inital location of abnormal proliferative cells. This is in line with what is biologically expected since abnormal cell proliferation can lead to forces acting in the tissue epithelium ^{1), (15}.

The crypt orifice deformation is an early sign of the formation of adenomas in the colon ^{17), (30}. Therefore, with this research, we were able to model how abnormal cells characterized by excessive proliferation can lead to crypt deformation that can occur before cancer appears ³⁰. This work suggests that we can understand where the abnormal cells are acting, observing the deformed orifice and crypt walls. This can be useful to force target drugs to act in specific colon epithelium regions in order to block the colorectal cancer formation or reducing its growth.

REFERENCES

• ¹
  A.A. Almet, H.M. Byrne, P.K. Maini & D.E. Moulton. The role of mechanics in the growth and homeostasis of the intestinal crypt. Journal of Theoretical Biology, 20 (2021), 585-608.
• ²
  A.A. Almet , B.D. Hughes, K.A. Landman, I.S. Näthke & J.M. Osborne. A Multicellular Model of Intestinal Crypt Buckling and Fission. Bulletin of Mathematical Biology, 80(2) (2018), 335-359.
• ³
  A.A. Almet , P.K. Maini , D.E. Moulton & H.M. Byrne . Modeling perspectives on the intestinal crypt, a canonical system for growth, mechanics, and remodeling. Current Opinion in Biomedical Engineering, 15 (2020), 32-39.
• ⁴
  K. Araki, T. Ogata, M. Kobayashi & R. Yatani. A Morphological Study on the Histogenesis of Human Colorectal Hyperplastic Polyps. Gastroenterology, 109 (1995), 1468-1474.
• ⁵
  R.P. Bird. Role of aberrant crypt foci in understanding the pathogenesis of colon cancer. Cancer Letters, 93 (1995), 55-71.
• ⁶
  R.P. Bird & C.K. Good. The significance of aberrant crypt foci in understanding the pathogenesis of colon cancer. Toxicology Letters, 112-113 (2000), 395-402.
• ⁷
  M. Bjerknes. Expansion of Mutant Stem Cell Populations in the Human Colon. Journal of Theoretical Biology , 178(4) (1996), 381-385.
• ⁸
  B.M. Boman, J.Z. Fields, O. Bonham-Carter & O.A. Runquist. Computer Modeling Implicates Stem Cell Overproduction in Colon Cancer Initiation. Cancer Research, 61(23) (2001), 8408-8411.
• ⁹
  P. Buske, J. Galle, N. Barker, G. Aust, H. Clevers & M. Loeffler. A Comprehensive Model of the Spatio-Temporal Stem Cell and Tissue Organisation in the Intestinal Crypt. PLoS Computational Biology, 7(1) (2011), e1001045.
• ¹⁰
  E. Carniel, M. Mencattelli & G.e.a. Bonsignori. Analysis of the structural behaviour of colonic segments by inflation tests: Experimental activity and physio-mechanical model. Proc Inst Mech Eng H., 229(11) (2015), 794-803.
• ¹¹
  G. De Matteis, A. Graudenzi & M. Antoniotti. A review of spatial computational models for multicellular systems, with regard to intestinal crypts and colorectal cancer development. Journal of Mathematical Biology, 66(7) (2013), 1409-1462.
• ¹²
  A. Di Garbo, M.D. Johnston, S.J. Chapman & P.K. Maini . Variable renewal rate and growth properties of cell populations in colon crypts. Physical Review E, 81(6) (2010), 061909.
• ¹³
  A. d’Onofrio & I.P.M. Tomlinson. A nonlinear mathematical model of cell turnover, differentiation and tumorigenesis in the intestinal crypt. Journal of Theoretical Biology , 244(3) (2007), 367-374.
• ¹⁴
  K. Drasdo & M. Loeffler . Individual-based models to growth and folding in one-layered tissues: intestinal crypts and early development. Nonlinear Analysis: Theory, Methods & Applications, 47(1) (2001), 245-256.
• ¹⁵
  C.M. Edwards & S.J. Chapman . Biomechanical Modelling of Colorectal Crypt Budding and Fission. Bulletin of Mathematical Biology , 69(6) (2007), 1927-1942.
• ¹⁶
  I.N. Figueiredo, C. Leal, G. Romanazzi & B. Engquist. Homogenization model for aberrant crypt foci. SIAM Journal on Applied Mathematics, 76(3) (2016), 1152-1177.
• ¹⁷
  I.N. Figueiredo , C. Leal , G. Romanazzi & B. Engquist . Biomathematical model for simulating abnormal orifice patterns in colonic crypts. Mathematical biosciences, 315 (2019), 108221.
• ¹⁸
  I.N. Figueiredo , C. Leal , G. Romanazzi , B. Engquist & P.N. Figueiredo. A convection-diffusionshape model for aberrant colonic crypt morphogenesis. Computing and Visualization in Science, 14(4) (2011), 157-166.
• ¹⁹
  J. Galle , M. Hoffmann & G. Aust . From single cells to tissue architecture -a bottom-up approach to modelling the spatio-temporal organisation of complex multi-cellular systems. Journal of Mathematical Biology , 58(1) (2009), 261-283.
• ²⁰
  D.V. Guebel & N.V. Torres. A computer model of oxygen dynamics in human colon mucosa: implications in normal physiology and early tumor development. Journal of Theoretical Biology , 250(3) (2008), 389-409.
• ²¹
  F. Hecht. BAMG: bidimensional anisotropic mesh generator. User Guide. Technical report, INRIA, Rocquencourt (1998).
• ²²
  F. Hecht . New development in FreeFem++. J. Numer. Math., 20(3-4) (2012), 251-265. URL https://freefem.org/
  » https://freefem.org/
• ²³
  P. Hogeweg. Evolving Mechanisms of Morphogenesis: on the Interplay between Differential Adhesion and Cell Differentiation. Journal of Theoretical Biology , 203(4) (2000), 317-333.
• ²⁴
  S.K. Kershaw, H.M. Byrne , D.J. Gavaghan & J.M. Osborne . Colorectal cancer through simulation and experiment. IET Systems Biology, 7(3) (2013), 57-73.
• ²⁵
  Y. Kuratani, S. Tamura, Y. Furuya & S. Onishi. Morphogenesis of a colorectal neoplasm with a type III S pit pattern inferred from isolated crypts. J Gastroenterol, 43 (2008), 597-602.
• ²⁶
  N.K. Kyslstad. “Simulating the viscoelastic response of the spinal cord”. Ph.D. thesis, Faculty of Mathematics and Natural Sciences, University of Oslo, Oslo (2014).
• ²⁷
  M.R. Nelson, J.R. King & O.E. Jensen. Buckling of a growing tissue and the emergence of twodimensional patterns. Mathematical Biosciences, 246(2) (2013), 229-241.
• ²⁸
  T. Roose, S.J. Chapman & P.K. Maini . Mathematical Models of Avascular Tumor Growth. SIAM REVIEW, 49(2) (2007), 179-208.
• ²⁹
  J. Salençon. “Handbook of Continuum Mechanics: General Concepts, Thermoelasticy”. Springer Science & Business Media (2012).
• ³⁰
  I.M.M. Van Leeuwen, H.M. Byrne , O.E. Jensen & J.R. King . Crypt dynamics and colorectal cancer: advances in mathematical modelling. Cell proliferation, 39(3) (2006), 157-181.
• ³¹
  S.Y. Wong, K.H. Chiam, C.T. Lim & P. Matsudaira. Computational model of cell positioning: directed and collective migration in the intestinal crypt epithelium. Journal of The Royal Society Interface, 7(3) (2010), S351-S363.

• 1
  Aberrant crypt foci (ACF) ^{5), (6} are clusters of crypts in the colon epithelium, containing cells with a deviant behavior with respect to the normal ones such as those hyper proliferating that due to their action can deform the crypt structure giving rise to aberrant crypt orifices. It is believed that ACF are precursors of colorectal cancer morphogenesis.

Publication Dates

• Publication in this collection
  04 Apr 2022
• Date of issue
  Jan-Mar 2022

History

• Received
  21 Mar 2021
• Accepted
  24 Aug 2021