Biomathematical model for simulating abnormal orifice patterns in colonic crypts

Mathematical Biosciences 315 (2019) 108221

Contents lists available at ScienceDirect

Mathematical Biosciences journal homepage: www.elsevier.com/locate/mbs

Biomathematical model for simulating abnormal orifice patterns in colonic crypts

T

Isabel N. Figueiredo ,a, Carlos Leala, Giuseppe Romanazzib, Björn Engquistc ⁎

a

CMUC, Department of Mathematics, Faculty of Sciences and Technology, University of Coimbra, Coimbra, Portugal Department of Applied Mathematics, Institute of Mathematics, Statistics and Scientific Computing (IMECC), State University of Campinas (UNICAMP), Brazil c Department of Mathematics and the Institute for Computational Engineering and Sciences, The University of Texas at Austin, Austin, USA b

ARTICLE INFO

ABSTRACT

Keywords: Colonic crypt Convection-diffusion Viscoelasticity

Colonic polyps, which are abnormal growths in the colon, are a major concern in colon cancer diagnosis and prevention. Medical studies evidence that there is a correlation between histopathology and the shapes of the orifices in colonic crypts. We propose a biomathematical model for simulating the appearance of anomalous shapes for the orifices of colonic crypts, associated to an abnormal cell proliferation. It couples a mechanical model that is a mixed elastic/viscoelastic quasi-static model describing the deformation of the crypt orifice, with a convection-diffusion model that simulates the crypt cell dynamics in space and time. The coupling resides in the variation of pressure generated by abnormal proliferative cells that induce a mechanical force and originate the change in shape of the crypt orifice. Furthermore the model is formulated in a two-dimensional setting, for emulating the top view of the colonic mucosa, observed in vivo in colonoscopy images. The primary focus of this study is on the modeling of this complex biological phenomenon, by defining an appropriate reduced biomathematical model. Additionally, a numerical procedure to determine its solution is also addressed. The overall numerical simulations indicate that an excess of cell proliferation, in different crypt locations, creates some of the anomalous patterns of the colonic crypt orifices, observed in vivo in medical images.

1. Introduction

the identification and early removal of adenomatous polyps can stop or slow its evolution [13]. To distinguish, among the adenomatous polyps, between those with a large potential to progress into cancer and those without any risk, a classification by Kudo [14,15] is normally used. It relies on the pattern of the pits of the crypts and consists of six different types: Type I, roundish (normal); Type II, star like or papillary; Type IIIL, tubular or round pits, larger than normal; Type IIIS, tubular or roundish pits smaller than normal; Type IV, sulcus-, branch- or gyrus-like; Type V, irregular or non-structural. Types I and II are considered normal, hyperplastic, inflammatory regions, whereas pit patterns of Types IIIL, IIIS, IV and V are considered neoplastic [16]. In Fig. 1 there is an illustration, taken from [1], that shows Types I, II and IIIL. In this work, which is a continuation of our previous works [17,18], we propose a biomathematical model to simulate the appearance of abnormal orifice patterns in colonic crypts. In order to derive this model, we assume that the deformation of the orifices is due to the large proliferative capacity acquired by epithelial cells inside and outside crypts. These abnormal cells can be originated by different biological perturbations (such as APC gene, mutations in KRAS, see for instance

The colon epithelium is densely filled by small cavities, called crypts, that are glands that descend from the intestinal lumen into the underlying connective tissue. The crypts have a cylinder tube shape, closed at the bottom and with an open orifice at the top and directed to the intestinal lumen. In the normal case, crypt orifices are circular, however when some abnormal biological alterations occur, such as genetic alterations (see for example [2]) the crypt geometry and its orifice can deform and change its shape. Essentially, these abnormal biological alterations can cause an uncontrolled cell proliferation and/ or block cell differentiation. This can lead to premalignant intestine lesions, such as aberrant crypt foci and small polyps [3–9]. In particular, aberrant crypt foci, which are clusters of colonic crypts with abnormal orifices and shapes, are believed to be the first stage of the formation of adenomatous polyps [10–12]. There is also medical evidence showing that most colorectal cancers have, as precursors, benign adenomatous polyps (abnormal protrusions in the colonic mucosa directed to the intestinal lumen) that can be detected by colonoscopy. Since the adenoma-carcinoma process is slow and can last many years,

⁎

Corresponding author. E-mail address: [email protected] (I.N. Figueiredo).

https://doi.org/10.1016/j.mbs.2019.108221 Received 9 April 2019; Received in revised form 20 June 2019; Accepted 26 June 2019 Available online 02 July 2019 0025-5564/ © 2019 Elsevier Inc. All rights reserved.

Mathematical Biosciences 315 (2019) 108221

I.N. Figueiredo, et al.

Fig. 1. Top row: Magnification chromoendoscopic images (in vivo medical images) showing different shapes of the colonic crypt orifices. Bottom row: Corresponding schematic illustrations of the different shapes. (a) and (d) Type I, (b) and (e) Type II, (c) and (f) Type IIIL. (Images from [1] with permission of Elsevier).

[19,20]), that alter the balance of the colonic cell renewal process. This is considered to be one of the possible causes of carcinogenesis [10]. The proposed biomathematical model is formulated in an appropriate two-dimensional domain, to imitate the top view of the colonic mucosa, observed in vivo in colonoscopy images, and couples a colonic cell evolution model with a mechanical model used to simulate the deformation of the epithelium tissue. The associated numerical simulations show a relation between abnormal cell proliferation, placed in different crypt locations, with different and abnormal pit patterns of the crypt. The colonic cell evolution model, is a mixed parabolic-elliptic system of partial differential equations, and simulates the dynamics of proliferative colonic cells. The unknowns of this model, in the spatial point x and at time t, are the density N(x, t) of proliferative cells (by proliferative cells we mean transit and/or semi-differentiated cells that also include the smaller subgroup of stem cells) and the pressure p(x, t) exerted by the cells due to their mitotic activity that induces a passive movement with velocity v. This convective cell velocity is modeled by a p (with ξ a parameter related to the perDarcy’s law, that is v = meability and the viscosity of the epithelium tissue, see [21]). It is known that the normal cell renewal process is associated to a rate of proliferation of semi-differentiated transit cells that is high at the bottom of the crypt and decreases along the crypt’s axis, being null in the last top third of the crypt [22,23]. Therefore, in this colonic cell evolution model, we also assume that the cell proliferation rate is dominated by the proliferation rate of semi-differentiated cells and depends only on the cell vertical position along the crypt’s axis and is constant over time. Since the pressure p is a consequence of the cell division process, then in the normal case, p decreases along the crypt axis towards the top orifice of the crypt, whereas the velocity increases along the same direction [11]. Genetic alterations may change the normal proliferative capacity of cells, and this can lead to an excess of proliferation. In these cases, the stability of the cell renewal process is lost giving rise to a different pressure. In our colonic cell evolution model the abnormal case is simulated by assuming that there are abnormal cells, located inside the crypt or in the inter-cryptal region, surrounding the crypt orifice, with a

high proliferation rate. This anomalous situation, according to our model, originates a pressure variation, with respect to the normal case, that generates a mechanical force responsible for the deformation of the normal round shape of the crypt orifice. In fact, it is this force, that we adopt as the coupling element, linking the colonic cell evolution model with the mechanical model that simulates the deformation of the crypt orifice. This mechanical model is a mixed elastic/viscoelastic quasi-static model, that gives the deformation of crypt top, that is considered to be an elastic medium, together with its inter-cryptal region, surrounding the crypt, considered to be viscoelastic (see for instance [24–27] for the justification of media used for biological tissues). The abnormal model of cellular evolution is obtained by the hypothesis that there are abnormal cells in a region of the crypt or in the inter-cryptal zone, where they should not be. These abnormal cells are characterized by having a higher proliferation rate with respect to the (normal) expected proliferation rate associated to their position. The biomathematical model, proposed in this paper, shows that the epithelium tissue containing a crypt deforms, in particular the shape of the crypt orifice deforms, when there are abnormal cells. After this introduction this paper has five more sections. In Section 2 we describe the colonic cell evolution model, then in Section 3 the mechanical model and the proposed biomathematical model in Section 4. Afterwards, the numerical procedure for the solution of the proposed model is explained in Section 5, and finally, in Section 6, the results are shown, including as well the discussion and conclusion subsection. 2. Colonic cell evolution model 2.1. General model It is known that the colon epithelium has a fast cell renewal, that lasts approximately only six-seven days in humans. The epithelium is characterized by a thin layer, whose thickness is equal to the diameter of a single cell, that fills the colon surface and the inner crypt walls. It is precisely along the crypt walls that most of the epithelium cell 2

Mathematical Biosciences 315 (2019) 108221

I.N. Figueiredo, et al.

proliferation occurs: it is high in the lower part of the crypt where semidifferentiated cells reside, and small at the top of crypt. Stem cells give rise to semi-differentiated cells, that are characterized to be proliferative and able to transit towards the top of the crypt up to two thirds of the crypt height, and along the crypt wall. Fully differentiated cells, that are not proliferative, are found in the third upper part of the crypt; these cells move further upwards in the direction of the top orifice, undergo apoptosis near the orifice, and then they are released in the lumen, ending the cell renewal process [11,22,28]. For defining the colonic cell evolution model, we suppose that the cells are grouped into two broad families, proliferative cells, that include semi-differentiated cells and also the smaller group of stem cells at the bottom of the crypt, and fully differentiated cells, with densities N1 and N2, respectively. Moreover, we suppose that both broad types of cells, with densities N1 and N2, follow a diffusive-convective flux with a p, constant diffusion coefficient D and with the same velocity v = where p is the pressure exerted by the cells due to their mitotic activity in the crypt walls (hereafter we assume = 1). Then, we adopt the following continuous model, similar to that used in [17,18], to characterize the colonic cell dynamics

N1 . ( pN1) = t N2 . ( pN2) = t N1 + N2 = 1.

. (D N1) + N1

into a malignant tumor): the top-down [29] and the bottom-up [30] - see also [31] for a related an interesting paper that includes a synthetic explanation of these two theories. For the top-down theory, dysplastic cells that reside in the inter-cryptal zone can spread laterally and invade healthy crypts. For the bottom-up theory, abnormal cells initially at the bottom of the crypt can migrate towards the top of the crypt and colonize/invade the inter-cryptal zone. Hereafter we rely on these biological hypothesis, and consequently, the abnormal case in our proposed model, means that the rate of proliferation α in (2), deviates from its normal definition, that is, it might be non zero in the inter-cryptal region or it might be anomalously high inside the crypt (in the normal case, as mentioned before, α is a decreasing function along the crypt height and is null in the third upper part of the crypt height). The following boundary conditions complement the equations (2)

p = 0, p N n

. ( pN ) =

N1 ,

. (D N2) + N1 , (1)

p =

. (D N ) + N N.

S,

N = N , in n

S,

(3)

where ∂S is the boundary of the hexagon in S1. Note that the first condition is a Dirichlet homogeneous condition for the pressure, imposing the zero value only in ∂S, and not in the crypt orifice or in the inter-cryptal region S1. In fact, let ∂S1 be the total boundary of S1, then 3: x 2 + x 2 = R } S = S1 {(x1, x2 , L) . The second boundary 1 2 condition is a Robin-type free-flux condition for the density N, that permits the proliferative cells to move beyond the boundary of the crypt domain, and to dislocate then into the neighbour epithelium tissue regions.

Here ∇. denotes the divergence operator, is the partial derivative t operator with respect to time, α is the proliferative rate of the cells with density N1 (in the model we assume that the proliferation rate is dominated by the proliferation rate of semi-differentiated cells), β is the transformation rate of proliferative cells (with density N1) into the fully differentiated cells (with density N2), and N1 + N2 = 1 is the no void condition [17,18]. Both functions α and β depend only on the cell position relative to the crypt vertical axis. In a normal cell renewal process, α is a decreasing function along the crypt height and it is null in the third upper part of the crypt height, whereas β has the reverse behavior, is null at the bottom of the crypt and large at the top. These rates allow then to have most of the proliferative cells at the bottom of the crypt, and fully differentiated cells at the top, as it is reported in the literature [22,28]. We remark that equations (1) incorporate the stem cells that generate the transit amplifying cells and are located always at the bottom part of the crypt. By summing the two first equations in (1), and defining N = N1, we have that the pressure p and the density N of proliferative cells verify

N t

D

in

2.2. Two-dimensional representation of a colonic crypt In order to reformulate the colonic cell evolution model (2)-(3) in a two-dimensional framework, to obtaining a model that simulates the top view of the colonic mucosa (see Fig. 1), we adopt a 2D Riemann manifold representation of the above domain S = S1 S2 S3 . To this end we firstly recall some definitions and notations. 3 be a 2D (two-dimensional) Riemann manifold and (U, φ) Let S 2 and 2 3. :U S its associated local chart, with U Any point x in U is denoted by x = (x1, x2) and to any scalar function 2 3 f: S we associate the function f¯ : U defined by f¯ = f . The Riemannian metric g connected to (U, φ) is defined by

g = (gij )ij = 1,2 ,

with

g ij = <

xi

,

xj

>

where < , > is the scalar product in

and the partial xj xi derivatives with respect to xi and xj, respectively. The determinant of g is denoted by det g and its inverse by g 1 = (g ij )ij = 1,2 . Then, the gradient of a scalar function f, defined on the manifold S, is

N, (2)

2

where Δ is the Laplace operator. Since our goal is to understand how a deregulation of the cell cycle can lead to a deformation of the crypt orifice, we model the colonic cell dynamics not only inside the crypt, where the cell renewal process occurs, but also in the inter-cryptal region at the epithelium surface. For this we admit that system (2) is valid in all the three-dimensional extended crypt domain depicted on the right of Fig. 2. This three-dimensional crypt domain contained in 3 (where a general point is x = (x1, x2 , x3) ), herein denoted by S, represents a three-dimensional crypt of height L, and is formed by three regions S = S1 S2 S3 in 3, where S1 is the hexagon-like inter-cryptal region at the height L with a circular hole (the orifice of the crypt) with radius R. The regions 3: x 2 + x 2 = R2 , 0 S3 = S2 = {(x1, x2 , x3) x3 L} and 1 2 3: x 2 + x 2 {(x1, x2 , 0) R2} are respectively the inner crypt wall and 1 2 the closed bottom of the crypt. At this point it is worth recalling the two existing theories for explaining the formation of an adenoma (that can subsequently develop

g ij

f= i, j = 1

2,

and

f¯ , xj x i

the divergence operator, div V, of a vector field V, defined on S, is

div V

·V =

1 | det g|

2 i=1

xi

( | det g| Vi ),

and finally the Laplacian operator Δf (or more exactly the Laplace–Beltrami operator) of a scalar function f, defined on the manifold S, is

f=

· f=

1 | det g|

2 i, j = 1

xi

| det g| g ij

f¯ . xj

We consider now the manifolds S1, S2 and S3, exhibited in Fig. 2 (right), with the associated charts (U1, U ), (U2, U ) and (U3, U ), 3 1 2 where the sets U1, U2 and U3 are represented in Fig. 2 (left), and 3

Mathematical Biosciences 315 (2019) 108221

I.N. Figueiredo, et al.

Fig. 2. Left: Two-dimensional representation of the crypt shown in the right. Right: three-dimensional representation of the geometry of a colonic crypt and its inter-cryptal region at the top. S1, shown in dark gray, is the top hexagon-like inter-cryptal region with a centered circular hole (the orifice of the crypt) of radius R, S2, shown in light gray, is a cylinder that represents the inner crypt wall and S3, shown in blue, is a circle that simbolizes the closed bottom of the crypt. S1, S2 and S3 correspond, respectively to U1, U2 and U3, in two-dimensions. The correspondence between left and right pictures is emphasized by the identical colors chosen for depicting each pair (Si, Ui), with i = 1, 2, 3 . Finally, r is the radius of the inner circle (associated to the base of the crypt), such that 0 < r < R. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

U1

U2

(x1, x2 ) = (x1, x2 , L),

in U1,

Rx1 , (x1, x2)

(x1, x2 ) =

1

Rx2 L , ( (x1, x2) (x1, x2) R r

r) ,

A11 (x1, x2) =

in U2, R (x x ) = (x1, x2 , 0), U3 1 2 r

in U3.

r R

(4)

ij,

in U1,

gij = ( 1)i + j 2

R r

gij =

r R

R2xi x j [x12

ij,

+

x 22]2

+

(R

+

x 22]2

, in U2,

det g =

2

in U2, in U3, in U1,

2

x 22 R (x1 x2)

2

2

in U2, in U3,

0

in U3,

A21 (x1, x2 ) =

x1 x2

in U1, R r hR (x1, x2)

0

with δij the Kronecker delta, and the determinants are defined by

det g = 1,

x12 R (x1, x2 )

A12 (x1, x2) = A21 (x1, x2 ),

L2xi x j r )2 [x12

2

2

1 x12 R r + 2 R h A22 (x1, x2 ) =

Here L > 0 represents the crypt height, r is the radius of the inner circle (associated to the base of the crypt), such that 0 < r < R, and (x1, x2 ) = x12 + x 22 . Thus we have the following expressions for the Riemannian metric g = (gij )ij = 1,2 in each region Uk (with k = 1, 2, 3):

gij =

in U1,

x 22 R r + R2 h

2

1 R2

in U2, in U3,

and

in U1, R2L2 r ) 2 (x12 + x 22),

in U2, b1 (x1, x2) =

0 in U1, x1 in U2, R2 0 in U3,

2.3. Reformulation of the colonic cell evolution model

b2 (x1, x2 ) =

Therefore, using the formulas for the gradient, divergence and Lapce–Beltrami operators indicated in the previous Section 2.2, the colonic cell evolution model (2) can be reformulated in two-dimensions in the set U = U1 U2 U3, as follows:

0 in U1, x2 in U2, R2 0 in U3.

Since in U1, the chart U1 is the identity function in (x1, x2), then the boundary conditions (3) defined in the boundary of S, ∂S, are also valid in the boundary of U, ∂U:

(R R r

det g =

4

,

in U3.

2

N t

Akl k, l = 1

2

Akl k,l= 1

xk

2

Akl k,l= 1

D

xk N xl

2p

x k xl

2

p xl

N

N k=1 2

+D

bk k=1

2

bk k=1

bk N +( xk

p = N, xk

p = xk

p = 0, p N n )N,

in U , N D = N , in U . n

(6)

Finally, the initial condition associated to (5) is the prescribed cell density at time t = 0, that is

N (x1, x2 , 0) = N0 (x1, x2),

(5)

where N0(x1, x2) is a fixed pre-defined function in the domain U.

where we removed the bar notation, that is, we used N instead of N¯ and p instead of p¯ . The coefficients Aij = Aij (x1, x2) and bi = bi (x1, x2) are respectively 4

(7)

Mathematical Biosciences 315 (2019) 108221

I.N. Figueiredo, et al.

Fig. 3. Ephitelial region filled by regular crypts. (a) Histological image of the human colon [32]. (b) Schematic picture representation of the colonic epithelium.

Fig. 4. (a) Two-dimensional representation U of the crypt geometry, used for defining the colonic cell evolution model (see Fig. 2). (b) The same two-dimensional representation U of the crypt geometry, where a slim ring, the region between Γ2 and Γ4 (shown in yellow), and a region (shown in pink) are depicted to represent, respectively, the region occupied by the last layer of cells (in the vertical direction) but still inside the crypt and the region of the crypt around the crypt orifice. (c) Domain Ω for the mixed elastic/viscoelatic model. It represents a portion of a top view of the colonic mucosa, and is composed of two different subdomains: Ωve made of a viscoelastic, and Ωe made of an elastic material. In particular, Ωve is in the region U1 defined in (a), and Ωe is the union of the pink region external to Γ2 in (b) with a close-up and a top view of the yellow slim ring defined in (b), that is, the region between Γ2 and Γ4. The light gray circle is a hole and represents the normal round crypt orifice. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

3. Crypt orifice deformation model

the last layer of cells (in the vertical direction), located at the extremity of the crypt at the top, but still inside the crypt. In Fig. 4(c) the global domain = ve e , represents a top view of the colonic mucosa, that involves the crypt orifice, shown in light gray, and includes the inter-cryptal region, for the mixed elastic/viscoelatic mode. It is composed of two different subdomains, Ωve (shown in brown), representing the inter-cryptal region and made of a viscoelastic material, and the other Ωe (shown in pink and yellow), around the crypt orifice at the top, with total thickness equal to the size of a single cell, and made of an elastic material. We remark that Ωe, in subfigure (c), corresponds to the union of the pink region external to Γ2 in (b) with a close-up and a top view of the yellow slim ring defined in (b), that is, the region between Γ2 and Γ4. We also note that Ωe is a thin ring that has thickness equal to the diameter of a cell and is filled by a layer of cells: the thickness of the yellow and pink region measures, each, half of a diameter of a cell. Finally, Σ1 is the boundary of the hexagon, Σ2 separates Ωe from Ωve, Σ3 represents the circle of radius R shown in subfigure (a), and the light gray region bounded by Σ4 is not included in Ω and represents the normal pattern of a crypt orifice. Finally we remark the analogy between the hexagonal domain depicted in Fig. 3 (real image in the left and corresponding synthetic picture in the right), showing the top view of the colonic mucosa, and the domain Ω described and shown in Fig. 4(c).

3.1. Schematic domains Fig. 3 (left) is an histological image showing a top view of the of colon epithelium, where one can visualize the cells located at the crypt orifice. This image indicates that the top of the crypts are distributed as in a honeycomb, with an hexagonal periodic domain. In Fig. 3 (right) we represent schematically this periodic and hexagonal crypt distribution. Each hexagon has a centred ring that defines the region where the crypt cells are located at the crypt orifice; outside this ring there is the inter-cryptal region that in normal cases does not have any proliferative colonic cell. In order to define the proposed biomathematical model (see Section 4), which is a coupled model connecting the colonic cell evolution model (5) and a mixed elastic/viscoelastic deformation model (see Section 3.2), we need now to describe in detail some parts of the domain U, already introduced in Section 2.2, and which is the domain for the colonic cell evolution model, and, in addition, to introduce a new domain for defining the mixed elastic/viscoelastic deformation model. This second domain emulates the hexagon displayed in Fig. 3, and it is by using this second domain that the proposed biomathematical model is able to capture different patterns for the top orifices of the crypts. In Fig. 4(a) it is shown the domain U of the colonic cell evolution model. It consists of three regions U1, U2 and U3 that represent, respectively, the inter-cryptal region at the top, the crypt wall and the base of the crypt. The curves Γ1, Γ2 and Γ3 represent the boundary of the hexagon, the boundary between U1 and U2, and the boundary between U2 and U3, respectively. In Fig. 4(b) the pink region (external with respect Γ2) represents the region occupied by the cells at the top of the crypt that are at the surface of the colon, directed to the lumen, while the yellow slim rim, between Γ2 and Γ4, represents the cells at the top of the crypt, that is,

3.2. Mixed elastic/viscoelastic deformation model The model used in this work for the deformation of the crypt orifice is based on the deregulation of the normal colonic cell renewal process, that, as already mentioned in Section 1, is connected to some genetic mutations that lead colonic cells to increase their proliferative capacity [3–5]. In this model, we assume that in an abnormal case, proliferative cells located in the crypt or at the inter-cryptal region, induce an abnormal high pressure, responsible for the orifice crypt deformation (we 5

Mathematical Biosciences 315 (2019) 108221

I.N. Figueiredo, et al.

note that a normal cell pressure cannot change the crypt shape). Thus, the deformation of the crypt orifice is modeled then, by a mixed elastic/ viscoelastic quasi-static equilibrium problem associated to some forces, which are the bond with the colonic cell evolution model defined in (5). These forces, represented by f, are defined by the variation of the difference between abnormal and normal pressure, p and p*, respectively, (p p*), with γ a pre-defined constant. that is f = We emphasize that the abnormal pressure p is obtained by solving the second equation of the colonic cell evolution problem (5) with an abnormal rate of proliferation α. The proposed model is defined in the space-time domain Ω × [0, T], where Ω is the flat space domain = e ve , represented in Fig. 4(c), and [0, T] the interval representing the time domain. By solving this model it is determined the displacement u (x , t ) = (u1 (x , t ), u2 (x , t )) in Ω, at each spatial point x ∈ Ω and at time t ∈ ]0, T[, by solving the following system e j ij (u )

=

fi

in

e

ve j ij (u )

=

fi

in

ve

u =

0

in

1

× ]0, T [,

0

in

2

× ]0, T [,

[[u]] = e ij (u ) nj

=

e ij (u ) nj

=

ve ij (u ) nj

fi

orifice Ωve, as an elastic material and a viscoelastic material respectively, is in good agreement with what is described in the literature for biological soft tissues, as well as for the colon tissue and colonic cells, see for instance [36–42]. 4. Biomathematical model The proposed biomathematical model is thus the coupled system (5)–(8), where (5) is associated to the boundary and initial conditions (6) and (7), respectively. The two problems are coupled by the cell pressure, obtained by solving the colonic cell evolution problem (5), and used in the definition of the force introduced in the mixed elastic/ viscoelastic problem (8). In order to solve this coupled system, we need to solve the problems (5) and (8) that are formulated in different domains, U and Ω, represented respectively in Fig. 4(a) and (c). In fact this discrepancy between the domains, where Ω changes with time, is one of the major difficulties for finding the solution of system (5)–(8). The introduction of the intermediate domain, shown in Fig. 4(b), is a strategy to overcome this difficulty by using an appropriate numerical approach, explained in Section 5. Problem (5) is non-trivial, but it can be solved analytically, using a fixed-point strategy similar to that presented in [18]. However, the coupled problem (5)–(8) is very complex to be solved in an analytical way, and to the best of our knowledge it does not exist in the literature an appropriate technique to solve it. Therefore we decide to present here only a numerical approach summarized in the next Section 5: there it is included a pseudocode that describes the steps for solving the coupled problem (5)–(8). Regarding the choice of the parameters in (5)–(8), we remark that the qualitative definition, reported in the literature, for the normal rates α and β, in (5), allows semi-differentiated cells to be located from the bottom to up to two thirds of the height of the crypt, and to have living cells with no reproductive capacity in the last third upper part of the crypt. Therefore, in this proposed biomathematical model (5)–(8), the distinction between normal and abnormal cases depends on the definition of the rate of proliferation α. In the normal case, α is a radial and decreasing function that is null in U1 (see Fig. 4(a) and also (20)). The so-called abnormal case is instead characterized by the appearance of cells with reproductive capacities greater than they should have, in irregular locations (deviant from the normal locations) in the domain U. This will be done using a Gaussian function that is added to the normal proliferation rate α* (see Section 6.2).

× ]0, T [, × ]0, T [,

on

2

× ]0, T [,

on

4

× ]0, T [.

(8)

In (8), we use the Einstein summation convention for repeated indices, n = (n1, n2 ) is the notation for the unit outward normal to the boundaries Σ2 or Σ4, the notation [[.]] symbolizes the jump, ∂i is the partial derivative with respect to xi, fi are the components of the force f, that is (f1 , f2 ) = ( 1 (p p*), 2 (p p*)), and σe and σve are the stress tensors, defining the constituent laws, for the elastic and viscoelastic materials, respectively (see for example [33–35]). These are defined by e ij (u)

= aijkl

ve ij (u)(t )

kl (u),

= bijkl (0)

t

kl (u )

0

bijkl s

(t

s)

kl (u ) ds ,

where the subscripts i, j, k, l vary in the set {1, 2}, (u) = (

(9) 2 ij (u ))ij = 1

is

the linearized strain tensor defined by ij (u) = 2 ( j ui + i uj ),aijkl are the coefficients of elasticity (a linear Hooke’s law is assumed) and bijkl(t) the coefficients of the so-called relaxation function, that defines the specific viscoelastic material (in this case a biological tissue) representing the effects of memory of the material. These are defined by 1

aijkl =

e ij kl

bijkl (t ) =

+ µe (

ve (t ) ij kl

ik jl

+

5. Approximate solution of the biomathematical model

il jk ),

+ µ ve (t ) (

ik jl

+

il jk )

(10)

The solution of the biomathematical model (5)–(8) is obtained by defining an appropriate numerical discretization of the system (5)–(8) (based on finite elements for the space discretization and on finite differences for the time discretization, as explained in Section 5.2) and by implementing it, as described in the following Section 5.1.

Here δ, with two subscripts, stands for the Kronecker delta, λe, μe are the Lamé coefficients for the elastic material occupying Ωe, and λve(t), μve(t) are coefficients of elasticity (for t = 0 ) and viscosity (for t > 0) of the viscoelastic material occupying Ωve (as it is usual in the applications, they are linear combinations of exponentials of time with negative coefficients, such that the memory of the material decreases exponentially when time increases: see Section 6.1 for the exact definition of λve(t), μve(t) used in our experiments). Still in (8) the condition u = 0 in Σ1 means the hexagon boundary is fixed. The condition ij (u) nj = fi , on Σ4, means that a traction force acts on the epithelium in Σ4 due to the pressure of cells inside the crypt. The transmission condition ije (u) nj = ijve (u) nj , on Σ2, means that no discontinuity of the constituent law is accepted across the boundary Σ2, between the crypt region Ωe and the inter-cryptal region Ωve. Finally, the condition [[u]] = 0 on Σ2 means that the displacement u is continuous in Σ2. Finally we remark that the choice of modeling the top of the crypt Ωe (which is a thin ring with thickness equal to the diameter of a cell and filled by a layer of cells) and the region surrounding the crypt’s

5.1. Pseudocode for the biomathematical model Denoting by 0 = t 0 < t1 < t2 <
Mathematical Biosciences 315 (2019) 108221

I.N. Figueiredo, et al.

of proliferation α is assumed, whereas in the normal case it is considered a normal rate of proliferation α*). Then, for each subsequent time tn, with n = 0, 1, 2, …, t final = T , do: Compute in the domain U, by using (5)-(6), the normal pressure p*(x1, x2, tn), by using N*(x1, x2, tn) and the normal rate of proliferation α* in the second equation of (5), the abnormal pressure p(x1, x2, tn), by using N(x1, x2, tn) and the abnormal rate of proliferation α in the second equation of (5), the normal cell density at subsequent time tn + 1, N * (x1, x2 , tn + 1), by using the normal pressure p*(x1, x2, tn) and the normal rate of proliferation α* in the first equation of (5), the abnormal cell density at subsequent time tn + 1, N (x1, x2 , tn+ 1), by using the abnormal pressure p(x1, x2, tn) and the abnormal rate of proliferation α in the first equation of (5). (p p*) in Ω. Compute the force f = The pressures p* and p computed in U, by the previous step, are interpolated in the domain Ω, knowing that the thin ring, between Γ2 and Γ4, shown in yellow in Fig. 4 (b), corresponds to the region between Σ3 and Σ4 in Ω, also shown in yellow in Fig. 4 (c). Orifice deformation: (a) Solve the mixed elastic/viscoelastic problem (8) in Ω to obtain the displacement field u(x1, x2, tn). (b) Use u(x1, x2, tn) to update the geometry of the domain Ω : (x1, x2) (x1, x2) + u (x1, x2 , tn), for all x ∈ Ω. Update the densities N and N*, and both pressures, p and p*, as well as, the coefficients in the current domain. Update n, n n + 1, and go back to the step 2.

2.

the total number of points in the spatial domain U), the basis functions of the finite element space, made of linear functions, and identifying the semi-discrete solutions Nn and pn of (11), with their finite element approximations, denoted by N fen and pfen , respectively, that is,

• • •

m

4.

5. 6.

N fen+ 1

t

, Ntest

Akl

Akl k,l= 1

N n+ 1

xk

2

Akl k,l= 1 2

D

bk

N n+1

k=1 2

Akl k,l= 1

xk

xk

D

N n+1

bk k=1

N n+1 , Ntest xl

pn , Ntest xk

Akl

bk

Akl

x k xl

( N n, ptest ) L2 (U ) ,

2

bk k=1

(12)

xk

,

s L2 (U )

+

s L2 (U )

) N fen+ 1,

bk k=1

pfen xk

s

,

, L2 (U )

s L2 (U )

s ) L2 (U ) .

(13)

by their definitions (12), the system (13) beSubstituting comes an algebraic problem, more exactly a matrix problem, where the unknowns are the coefficients of N fen and pfen , that is, the following two vectors in m

(N fen,

(Nsn )m s= 1

pfen )

(psn )m . s=1

and

(14)

These components represent, respectively, the density of proliferative cells and the pressure, at time tn and in the spatial location defined by node s of the finite element discretization. 5.2.2. Mixed elastic/viscoelastic model The discretization of (8) in Ω is obtained by applying as well the finite element method, and the procedure is similar to that described for (5). In fact, defining the following mixed stress tensor e ij (u) ve ij (u)

=

in

e,

in

ve ,

(15)

the variational formulation of (8), for each time tn, is 2 i=1

L2 (U )

n j ij (u ),

( 2 i=1

(

i (p

n

ui, test ) L2 ( )= p*n ), ui, test ) L2 ( ) .

(16)

n

Here u is the displacement u at time tn (u depends implicity on time), ui,test, for i = 1, 2, is an arbitrary function of H1(Ω) and (.,.) L2 ( ) is the standard inner product in L2(Ω). By decomposing Ω with a finite element mesh, and using linear finite elements, associated to the basis ϕs, with s = 1, …, m , and replacing in (16), un by its finite element apm proximation ufen = s = 1 (u1,ns , u 2,n s ) s , and ui,test by ϕs for i = 1, 2, then (16) reduces to

L2 (U )

pn ,p xk test

,

2

x k xl

k,l= 1

pfen

bk k=1

xl

2pn fe

2

= ( N fen,

xl

+(

xk

k=1

, L2 (U )

2 pn

psn s ,

2

N fen+ 1

N fen+ 1

D

N fen+ 1

2

D

pfen

N fen+ 1

xk

k,l= 1

+

) N n+ 1, Ntest

+(

s=1

L2 (U )

2

=

2

pn xl

s

xk

k,l= 1

ij (u)

2

,

2

L (U )

2

=

N fen t

5.2.1. Colonic cell evolution model To discretize (5)–(7), firstly a variational formulation of the colonic cell evolution model is obtained by multiplying, respectively, the first and second equations in (5) by test functions Ntest and ptest, defined in the Sobolev spaces H1(U) and H01 (U ), respectively, and by integrating the result in the domain U. Secondly, the variational semi-discretization of (5) is obtained by replacing the exact time derivative by an approximate derivative, which using the backward Euler method. Thus, denoting by 0 = t 0 < t1 < t2 <
Nn

pfen =

pn

the complete discretization of (5)–(7) consists in replacing, in (11), (Nn, pn) by (N fen, pfen ), and both Ntest and ptest by ϕs, respectively, which yields

5.2. Discretization of the biomathematical model

N n+1

Nsn s , s=1

• 3.

m

N fen =

Nn

= L2 (U )

(

(11)

where the notation (.,.) L2 (U ) represents the inner product in the functional space L2(U). Thirdly, the complete discretization of (5) corresponds to considering a finite element decomposition of the domain U. Therefore, denoting by ϕs, with s = 1, …, m (m finite and representing

n j ij (u fe ),

s ) L2 ( )

=

n i (pfe

pfe*n ),

s

, L2 ( )

for i = 1, 2.

(17)

This equation yields an algebraic system, whose unknown is the twocomponent vector

(u1,ns , u2,n s )m s= 1 7

Mathematical Biosciences 315 (2019) 108221

I.N. Figueiredo, et al.

representing the plane displacement u of the spatial node s in Ω, at time tn.

[44] and Matlab [45] were used. The normal shape of the crypt orifice is depicted in subfigure (a) of Fig. 5. In each one of these five cases the abnormal assumption consists in assuming that the rate of proliferation α is anomalous. In contrast to the normal case, where α is high at the bottom of the crypt and zero at the top, we now suppose that α is different from zero at the top of the crypt or at the inter-cryptal region. In all the cases we identify the abnormal α as the sum of the normal rate of proliferation, denoted by α* and defined in (20), with a Gaussian function (the perturbation, denoted by αper). We remark that this option for αper is merely mathematical. However, since the quantity to be perturbed is the rate of proliferation of a density of cells, we believe this choice is appropriate and simulates quite well a small alteration in the normal rate. The five cases are associated to a proliferative rate with a different support region for αper as shown in the second row pictures in Fig. 5. The list of the five cases and the resulting deformed shapes of the orifices are as follows:

6. Results 6.1. Parameter values Regarding the values of the Lamé coefficients for the elastic material represented by the domain Ωe we first remark that they are related to the Young’s modulus E and Poisson’s ratio ν, by e

=

E (1 + )(1

2 )

, µe =

E . 2(1 + )

(18)

Therefore we use the values E = 0.936 × 10 kPa (kilopascal) and = 0.3 for the colon measured in the works [38,43]. For the coefficients λve(t) and μve(t), defining the viscoelastic material in the domain Ωve, we consider 9

ve (t )

t

= 2 e (1 + e

),

µ ve (t ) = 2µe (1 + e

t

),

Experiment 1 Support of αper: eight circular regions equally distributed and located in the crypt wall (Fig. 5(b)). Maximum of αper: 10 10 . Radius of each circle: 0.1/32. The shape of the crypt orifice is shown in Fig. 5(l) and corresponds to a Type I normal crypt, with a larger orifice than the normal case. Experiment 2 Support of αper: four circular regions located in the intercryptal region (Fig. 5(c)). Maximum of αper: 10 8 . Radius of each circle: 0.1/8. The crypt orifice is shown in Fig. 5(m) and it exhibits a star like shape very similar to a Type II pit pattern of nonneoplastic, hyperplastic polyp, radially symmetric [1,46]. Experiment 3 Support of αper: one circular region positioned in the inter-cryptal region (outside the crypt at the right) with a high rate of proliferation (Fig. 5(d)). Maximum of αper: 10 8 . Radius of the circle: 0.1/8. Subfigure (n) of Fig. 5 displays the deformation of the crypt orifice and shows an elongated (tubular) shape corresponding to a Type IIIL pit pattern. Experiment 4 Support of αper: two circular regions situated in the intercryptal region (outside the crypt at opposite sides- left and right) with a high rate of proliferation (Fig. 5(e)). Maximum of αper: 10 8 . Radius of each circle: 0.1/8. The deformation of the orifice is displayed in Fig. 5(o) and corresponds to a Type IIIL pit. Experiment 5 Support of αper : two circular regions close together and placed in the inter-cryptal region (outside the crypt) with a high rate of proliferation (Fig. 5(f)). Maximum of αper : 10 8 . Radius of each circle: 0.1/8. Subfigure (p) of Fig. 5 shows the deformed orifice of Type IIIL.

(19)

with = 5.430 as measured in [25]. We note that this choice originates at time t = 0, ve (0) = 2 e and µ ve (0) = 2µe , which means that the material in Ωve is more stiff than that in Ωe. As in [17,18], we use the following expressions for α* (the normal rate of proliferation) and for β (the transformation rate of proliferative cells into the fully differentiated cells) in a normal cell renewal process

* (x3) =

2

2 L 3

x3

0

0 0 (x3) =

0 2 L 3

x3

2 L, 3 elsewhere,

2

x3 <

2 L, 3

x3 <

elsewhere,

(20)

where the value = 10 is obtained by forcing the cell velocity to be equal to 0.85 position per hour at the top of the crypt [11]. In addition, x3 is the third dimension, parallel to the crypt vertical axis (see Fig. 2, right and left) that can be expressed as a function of (x1, x2) in the domain U with values in [0, L] (see as well (4)), as follows 6

L x3 (x1, x2 ) =

L R 0

r

( (x1, x2)

if (x1, x2 )

U1,

r ) if (x1, x2 )

U2,

if (x1, x2 )

U3.

In addition, we nondimensionalize the geometry of the crypt and assume that the edge of the hexagon is equal to 1 and set R = 1/2, r = 1/4 . Then, by using [12] we set the crypt height L = 12.6979 . The initial condition at time t = 0 for a normal cell density, in (7), is N0 (x1, x2 ) = N (x1, x2 , 0) defined by

N0 (x1, x2 ) =

1 if 0

x3 (x1, x2)

2 L, 3

0 elsewhere.

These numerical simulations clearly show a correlation between the location of the abnormal cell proliferation and the shape of the deformed crypt orifice. The difference of pressure depicted in the third (p p*) row of Fig. 5 is responsible through the internal force f = for the deformation observed in the last row of Figure 5. Furthermore, these simulations can also predict the evolution in time of these aberrant orifice shapes. By way of illustration Fig. 6 displays the distribution of the proliferative cell density N, at final time, for experiments 2, 3 and 5, in the domain = e ve , that is, the proliferative cell distribution in the top of the crypt. Note that the proliferative cells accumulate near the crypt orifice in particular in experiments 3 and 5 (see subfigures (b) and (c)). Finally, and as an example, Fig. 7 shows the difference of pressure p p* in the domain U, for experiment 3, at the initial and final times.

(21)

Moreover, in (5) we set the parameter D = 1, and in (8) for the force f= (p p*) we use = 1. We also remark that the normal pressure, denoted by p*, is obtained by solving the elliptic partial differential equation in (5) for a normal cell density distribution (as for instance the prescribed cell density N0 at time t = 0 in (21)), with a normal rate of proliferation α*. 6.2. Experiments – numerical simulations We describe the results obtained with the proposed biomathematical model for five cases, that yield five different (abnormal) patterns of the crypt orifice. For all the experiments the software FreeFem++ 8

Mathematical Biosciences 315 (2019) 108221

I.N. Figueiredo, et al.

Fig. 5. (a) Normal orifice of the crypt at the initial time and the finite element mesh. Column 1 – (b), (g) and (l) – Experiment 1: Type I pit, roundish pit [normal colonic mucosa]. Column 2 – (c), (h) and (m) – Experiment 2: Type II pit, star-like deformation. Column 3 – (d), (i) and (n) – Experiment 3: Type IIIL pit, elongated shape. Column 4 – (e), (j) and (o) – Experiment 4: Type IIIL pit, tubular shape. Column 5 – (f), (k) and (p) Experiment 5: Type IIIL pit, elongated shape. The colours, in the second and third rows, indicate the intensity (the highest intensity is shown in dark red followed, in decreasing order, by red, orange, yellow, green, light blue, blue and dark blue). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

We recall that in this case the abnormal pressure p, which depends implicitly on time (through the density of cells N) is obtained by solving, for each time step, the second equation in (5), assuming that the abnormal rate of proliferation α is the sum of the normal rate of proliferation, denoted by α* and defined in (20), with the perturbation αper shown in Fig. 5(d).

crypts. The proposed model is highly complex, mathematically and computationally, and involves a mechanical model coupled with a mass/transport model, defined in two different domains that have to be appropriately combined. Among the six different types of pit patterns described in Kudo’s classification scheme, and which have a direct correlation with histopatology, the proposed model is able to retrieve three: Type I roundish normal orifice, Type II star like orifice, and Type IIIL orifice with an elongated tubular shape. The overall results clearly indicate that the location of proliferative cells in abnormal positions lead to these three different types. Type I is obtained when an equally spaced distribution

6.3. Discussion and conclusions The purpose of this study was to develop a biomathematical model suited for predicting the appearance of abnormal pit patterns in colonic 9

Mathematical Biosciences 315 (2019) 108221

I.N. Figueiredo, et al.

Fig. 6. Distribution of the cell density N at final time, for experiment 2 in (a), for experiment 3 in (b), and for experiment 5 in (c). The colours indicate the intensity (the highest intensity is shown in dark red followed, in decreasing order, by red, orange, yellow, green, light blue, blue and dark blue). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

order to account for this type of behavior. In addition we note that in our previous work [17] we have studied, using a similar cellular dynamics model, the crypt deformation along the missed dimension of the crypt axis direction. Then, such a model, incorporated in the current proposed model for the crypt orifice deformation, could be used to develop a full three dimensional model of the crypt deformation. We also point out that the proposed biomathematical model depends on several parameters, some associated with the properties of the biological tissues of the crypt and inter-cryptal region, as well as, the cell rates α and β, that have a key role in the obtention of different shapes for the crypt orifices. By relying only on qualitative information for these different parameters (that is the only one available in the literature), which is obviously a strong limitation, it was even though possible to replicate different shapes for the crypt orifices. This fact also shows the right capability of the proposed model. Finally, to the best of our knowledge, this is the first work reporting a consistent mathematical model with numerical simulations that can estimate aberrant patterns for the orifices of colonic crypts.

Fig. 7. Difference of pressure p p* at the initial time (a) and at final time (b), in the domain U, for experiment 3. The colours indicate the intensity (the highest intensity is shown in dark red followed, in decreasing order, by red, orange, yellow, green, light blue, blue and dark blue). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

of abnormal cells is located inside the crypt, along its wall. Type II is obtained when four equally spaced distributions of abnormal cells are located outside the crypt orifice. Type IIIL is obtained when abnormal cells are located in one or two position in the inter-cryptal region. We emphasize that in this work we mention the top-down [29] and the bottom-up [30] phenomena, because these provide biological arguments that support the appearance of abnormal cells either at the top or inside the crypt. These locations for abnormal cells are used in our proposed model and lead to the results described in the paper. However in this paper our intention is not to justify biologically why abnormal cells are located in these irregular locations: we just assume they might exist there because of the top-down or bottom-up phenomena. We note that our main goal was to design an appropriate and somewhat realistic model able to retrieve, in silico, the abnormal shapes of the crypt orifices, observed in vivo in colonoscopy images. On the whole, the proposed model reached this main objective, based on qualitative biological information (for defining the parameters) and also some logical (but biological experiment-free) assumptions. In the future it would be important to promote interdisciplinary collaboration between biologists and mathematicians, aiming at a biological-experimental validation and improvement of the proposed biomathematical model. Despite the complexity of the proposed biomathematical model, it should be noted that it is also a reduced model, from the biological point of view, that intends to mimic the real biological phenomenon in a simplified way. Furthermore, the model is formulated in a two-dimensional framework, in order to imitate the top view of the colon, observed in vivo, in colonoscopy exams. Consequently it can not describe, correctly, from the geometrical point of view, deformations of the inner crypt wall. Moreover we emphasize, that such internal deformations can be responsible for a crypt budding process that can deform or reduce the crypt orifice dimension, and possibly originate patterns of the Types IIIS, IV and V of Kudo’s classification scheme [5], and that were not captured by our proposed model. In fact, our model gives only the deformation of the top of the crypt and of the intercryptal region, and does not give the deformation along the inner crypt walls. Then only deformations occurring outside the crypt orifice can be described by our proposed model. It could be an interesting subject, to extend in the future this model to a three-dimensional framework in

Declaration of competing interest None. Acknowledgment This work partialy was supported by the grant UID/MAT/00324/ 2019 of FCT (Fundação para a Ciência e a Tecnologia, Portugal), and also by FAPESP 2017/05519-5 and FAEPEX-Pesquisa 519.292. References [1] Q. Huang, N. Fukami, H. Kashida, T. Takeuchi, E. Kogure, T. Kurahashi, E. Stahl, Y. Kudo, H. Kimata, S.-e. Kudo, Interobserver and intra-observer consistency in the endoscopic assessment of colonic pit patterns, Gastrointest. Endosc. 60 (4) (2004) 520–526. [2] H. Kang, D. Shibata, Direct measurements of human colon crypt stem cell niche genetic fidelity: the role of chance in non-Darwinian mutation selection, Front. Oncol. 3 (2013) 264. [3] A. Humphries, N.A. Wright, Colonic crypt organization and tumorigenesis, Nat. Rev. Cancer 8 (6) (2008) 415–424. [4] T. Reya, H. Clevers, Wnt signalling in stem cells and cancer, Nature 434 (7035) (2005) 843–850. [5] M. Van De Wetering, E. Sancho, C. Verweij, W. De Lau, I. Oving, A. Hurlstone, K. Van Der Horn, E. Batlle, D. Coudreuse, A.-P. Haramis, et al., The β-catenin/tcf-4 complex imposes a crypt progenitor phenotype on colorectal cancer cells, Cell 111 (2) (2002) 241–250. [6] S.P. Bach, A.G. Renehan, C.S. Potten, Stem cells: the intestinal stem cell as a paradigm, Carcinogenesis 21 (3) (2000) 469–476. [7] H. Clevers, The intestinal crypt, a prototype stem cell compartment, Cell 154 (2) (2013) 274–284. [8] H. Clevers, R. Nusse, Wnt/β-catenin signaling and disease, Cell 149 (6) (2012) 1192–1205. [9] L. Ricci-Vitiani, E. Fabrizi, E. Palio, R. De Maria, Colon cancer stem cells, J. Mol. Med. 87 (11) (2009) 1097–1104. [10] E.R. Fearon, B. Vogelstein, A genetic model for colorectal tumorigenesis, Cell 61 (5) (1990) 759–767. [11] I.M.M. van Leeuwen, H.M. Byrne, O.E. Jensen, J.R. King, Crypt dynamics and colorectal cancer: advances in mathematical modelling, Cell Prolif. 39 (2006) 157–181. [12] D.V. Guebel, N.V. Torres, A computer model of oxygen dynamics in human colon mucosa: implications in normal physiology and early tumor development, J. Theor. Biol. 250 (3) (2008) 389–409. [13] J. Bond, Clinical evidence for the adenoma-carcinoma sequence, and the

10

Mathematical Biosciences 315 (2019) 108221

I.N. Figueiredo, et al.

[14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25]

[26] [27] [28] [29]

management of patients with colorectal adenomas. Proceedings of the Seminars in Gastrointestinal Disease, 11 (2000), pp. 176–184. S.-E. Kudo, S. Tamura, T. Nakajima, H.-o. Yamano, H. Kusaka, H. Watanabe, Diagnosis of colorectal tumorous lesions by magnifying endoscopy, Gastrointest. Endosc. 44 (1) (1996) 8–14. S.-E. Kudo, Early Colorectal Cancer: Detection of Depressed Types of Colorectal Carzcinoma, 166 Igaku-Shoin, Tokyo, 1996. M. Li, S.M. Ali, S. Umm-a OmarahGilani, J. Liu, Y.-Q. Li, X.-L. Zuo, Kudo’s pit pattern classification for colorectal neoplasms: a meta-analysis, World J. Gastroenterol. WJG 20 (35) (2014) 12649–12656. I.N. Figueiredo, C. Leal, G. Romanazzi, B. Engquist, P.N. Figueiredo, A convectiondiffusion-shape model for aberrant colonic crypt morphogenesis, Comput. Vis. Sci. 14 (4) (2011) 157–166. I.N. Figueiredo, C. Leal, G. Romanazzi, B. Engquist, Homogenization model for aberrant crypt foci, SIAM J. Appl. Math. 76 (3) (2016) 1152–1177. R. Lambert, D. Provenzale, N. Ectors, H. Vainio, M. Dixon, W. Atkin, M. Werner, S. Franceschi, H. Watanabe, G. Tytgat, et al., Early diagnosis and prevention of sporadic colorectal cancer, Endoscopy 33 (12) (2001) 1042–1064. C. Albuquerque, M. Cravo, C. Cruz, P. Lage, P. Chaves, P. Fidalgo, A. Suspiro, C.N. Leitao, Genetic characterisation of patients with multiple colonic polyps, J. Med. Genet. 39 (4) (2002) 297–302. H. Greenspan, On the growth and stability of cell cultures and solid tumors, J. Theor. Biol. 56 (1) (1976) 229–242. D. Drasdo, M. Loeffer, Individual-based models to growth and folding in onelayered tissues: intestinal crypts and early development, Nonlinear Anal. 47 (2001) 245–256. C. Gaspar, R. Fodde, Apc dosage effects in tumorigenesis and stem cell differentiation, Int. J. Dev. Biol. 48 (5-6) (2004) 377–386. M. Nelson, J. King, O. Jensen, Buckling of a growing tissue and the emergence of two-dimensional patterns, Math. Biosci. 246 (2) (2013) 229–241. E.L. Carniel, M. Mencattelli, G. Bonsignori, C.G. Fontanella, A. Frigo, A. Rubini, C. Stefanini, A.N. Natali, Analysis of the structural behaviour of colonic segments by inflation tests: experimental activity and physio-mechanical model, Proc. Instit. Mech. Eng. Part H J. Eng. Med. 229 (11) (2015) 794–803. K.K. Chawla, M. Meyers, Mechanical Behavior of Materials, Prentice Hall, Upper Saddle River, 1999. T. Roose, S.J. Chapman, P.K. Maini, Mathematical models of avascular tumor growth, SIAM Rev. 49 (2) (2007) 179–208. A. Di Garbo, M. Johnston, S. Chapman, P. Maini, Variable renewal rate and growth properties of cell populations in colon crypts, Phys. Rev. E 81 (6) (2010) 061909. I.-M. Shih, T.-L. Wang, G. Traverso, K. Romans, S.R. Hamilton, S. Ben-Sasson, K.W. Kinzler, B. Vogelstein, Top-down morphogenesis of colorectal tumors, Proc.

Natl. Acad. Sci. 98 (5) (2001) 2640–2645. [30] S.L. Preston, W.-M. Wong, A.O.-O. Chan, R. Poulsom, R. Jeffery, R.A. Goodlad, N. Mandir, G. Elia, M. Novelli, W.F. Bodmer, et al., Bottom-up histogenesis of colorectal adenomas: origin in the monocryptal adenoma and initial expansion by crypt fission, Cancer Res. 63 (13) (2003) 3819–3825. [31] I.M. Van Leeuwen, C.M. Edwards, M. Ilyas, H.M. Byrne, Towards a multiscale model of colorectal cancer, World J. Gastroenterol. WJG 13 (9) (2007) 1399. [32] M. Hill, Embryology gastrointestinal tract – colon histology. Retrieved (2019, january 17) from https://embryology.med.unsw.edu.au/embryology/index.php/ gastrointestinal_tract_-_colon_histology. [33] G. Duvaut, J.L. Lions, Inequalities in Mechanics and Physics, 219 Springer Science & Business Media, 2012. [34] S. Shaw, J.R. Whiteman, Numerical solution of linear quasistatic hereditary viscoelasticity problems, SIAM J. Numer. Anal. 38 (1) (2000) 80–97. [35] M. Fabrizio, An existence and uniqueness theorem in quasi-static viscoelasticity, Q. Appl. Math. 47 (1) (1989) 1–8. [36] A.E. Bharucha, R.D. Hubmayr, I.J. Ferber, A.R. Zinsmeister, Viscoelastic properties of the human colon, Am. J. Physiol. Gastroint. Liver Physiol. 281 (2) (2001) G459–G466. [37] J. Ranft, M. Basan, J. Elgeti, J.-F. Joanny, J. Prost, F. Jülicher, Fluidization of tissues by cell division and apoptosis, Proc. Natl. Acad. Sci. 107 (49) (2010) 20863–20868. [38] S. Kawano, M. Kojima, Y. Higuchi, M. Sugimoto, K. Ikeda, N. Sakuyama, S. Takahashi, R. Hayashi, A. Ochiai, N. Saito, Assessment of elasticity of colorectal cancer tissue, clinical utility, pathological and phenotypical relevance, Cancer Sci. 106 (9) (2015) 1232–1239. [39] A. Fletcher, C. Breward, S. Chapman, Mathematical modeling of monoclonal conversion in the colonic crypt, J. Theor. Biol. 300 (2012) 118–133. [40] C.M. Edwards, S. Chapman, Biomechanical modelling of colorectal crypt budding and fission, Bull. Math. Biol. 69 (6) (2007) 1927–1942. [41] F.A. Meineke, C.S. Potten, M. Loeffler, Cell migration and organization in the intestinal crypt using a lattice-free model, Cell Prolif. 34 (4) (2001) 253–266. [42] S. Leeman, J. Jones, Visco-elastic models for soft tissues, Acoustical Imaging, Springer, 2008, pp. 369–376. [43] M. Yu, J. Wang, H. Wang, L. Liu, Y. Yan, J. Zhang, Y. Liang, S. Dong, Calculation of the intracellular elastic modulus based on an atomic force microscope micro-cutting system, Chin. Sci. Bull. 57 (15) (2012) 1868–1872. [44] F. Hecht, New development in FreeFem++, J. Numer. Math. 20 (3-4) (2012) 251–266. [45] The Mathworks, Inc., http://www.matlab.com. [46] S. Kudo, C. Rubino, C. Teixeira, H. Kashida, E. Kogure, Pit pattern in colorectal neoplasia: endoscopic magnifying view, Endoscopy 33 (04) (2001) 367–373.

11