Biomechanical Constitutive Modeling of the Gastrointestinal Tissues: Where Are We?

The gastrointestinal (GI) tract is a continuous channel through the body that consists of the esophagus, the stomach, the small intestine, the large intestine, and the rectum. Its primary functions are to move the intake of food for digestion before storing and ultimately expulsion of feces from the rectum through the anal sphincter. The mechanical behavior of GI tissues thus plays a crucial role for GI function in health and disease. The mechanical properties are typically characterized by a constitutive biomechanical model, which is a mathematical representation of the relation between load and deformation in a tissue. Hence, validated biomechanical constitutive models are essential to characterize and simulate the mechanical behavior of the GI tract under physiological and pathological conditions. Numerous constitutive models have consequently been proposed over the past three decades, mainly inspired by work done in cardiovascular tissues. Here, a comprehensive review of these constitutive models is provided. This review is limited to studies where a model of the strain energy function is proposed to characterize the stress-strain relation of a GI tissue. Several needs are identified for more advanced modeling of GI biomechanics including: 1) Microstructural models that provide actual structure-function relations; 2) Validation of coupled electro-mechanical models accounting for active muscle contractions; 3) Human data under physiological and pathological conditions to develop and validate models. The findings from this review provide guidelines for using existing constitutive models as well as perspective and directions for future studies aimed at establishing new constitutive models for GI tissues.

Its primary functions are to move the intake of food for digestion before storing and ultimately expulsion of feces from the rectum through the anal sphincter. The mechanical behavior of GI tissues thus plays a crucial role for GI function in health and disease. The mechanical properties are typically characterized by a constitutive biomechanical model, which is a mathematical representation of the relation between load and deformation in a tissue. Hence, validated biomechanical constitutive models are essential to characterize and simulate the mechanical behavior of the GI tract under physiological and pathological conditions. Numerous constitutive models have consequently been proposed over the past three decades, mainly inspired by work done in cardiovascular tissues. Here, a comprehensive review of these constitutive models is provided. This review is limited to studies where a model of the strain energy function is proposed to characterize the stress-strain relation of a GI tissue. Several needs are identified for more advanced modeling of GI biomechanics including: 1) Microstructural models that provide actual structure-function relations; 2) Validation of coupled electro-mechanical models accounting for active muscle contractions; 3) Human data under physiological and pathological conditions to develop and validate models. The findings from this review provide guidelines for using existing constitutive models as well as perspective and directions for future studies aimed at

Introduction
The gastrointestinal (GI) tract is a long conduit extending from the oral cavity to the anus and consists of five major segments: the esophagus, the stomach, the small intestine, the large intestine, and the rectum (Figure 1). Each segment is separated by sphincters such as the lower, esophageal sphincter, pylorus, and anal sphincters. Although there are 5 anatomical similarities between these segments of the GI tract, there are also differences and in particular differences in function. The GI tract acts as a digestive mechanical system that transforms food into chyme and ultimately feces during transit trough the GI system. Healthy functioning of the GI tract is crucial for absorbing essential nutrients (carbohydrates, proteins, fats, water, minerals, and vitamins) from food as well as for 10 secretion of chemical substances that aid breakdown of the food constituents [1]. Most drugs are administered orally, and their efficacy also relies on the healthy functioning of the GI [2]. The GI functions are affected by numerous structural and functional GI diseases and disorders such as gastro-esophageal reflux disease (GERD), dyspepsia, irritable bowel syndrome (IBS), diverticulosis, constipation, diarrhea, hemorrhoids, and 15 fecal incontinence. These diseases and disorders have a significant impact on quality of life and impose a major financial burden on patients and healthcare systems [3][4][5][6][7][8][9].
Unfortunately, these GI diseases and disorders are highly prevalent with functional GI diseases alone affecting over 40% of the world population [10].
Most of these diseases and disorders are invariably linked to altered GI biomechanics 20 [11][12][13]. Studying the deformation of GI organs and the forces that causes such deformation is thus crucial for understanding their mechanical behavior in health and disease.
The mechanical properties are also essential for developing computational models of GI organs that can generate hypotheses, support diagnosis, predict treatment outcomes, and virtually evaluate new therapeutic approaches. Mathematically, deformation is typically 25 quantified in terms of strain, which is a measure of stretch relative to a reference state, while force is quantified in terms of stress, which is a measure of force per unit surface area. A mechanical constitutive model is a mathematical formulation that relates stress and strain in a material through parameters [14]. It is essential for characterizing the mechanical behavior of a tissue quantitatively since stress cannot be directly measured. It 30 also constitutes one of the three major elements of biomechanical computational models, along with tissue geometry and boundary conditions.
In the present work, we provide a comprehensive review of the published constitutive models of the five GI organs (GI sphincters were excluded). The focus is on studies where a constitutive model is established based on a strain energy function since this approach 35 has been shown to capture well the non linear and large deformation of soft tissues, including GI tissues. The GI tract consists of four main layers all with different mechanical characteristics: submucosa, mucosa, muscle, and serosa. In the review, we consider studies that model intact GI tissue and/or any combination of the individual layers. Over the last three decades, numerous strain-energy derived constitutive models have been proposed 40 for the GI tissues, primarily inspired by models and methods validated in cardiovascular tissues. After summarizing the most common models and their formulations, we conduct a cross-analysis of all the studies comparing them based on several factors such as organ studied, model type, testing protocol followed, species used for the tissue, etc. The objectives are two-fold: 1) Provide an overview of the validated constitutive models of GI 45 tissues currently established to assist in the selection of a suitable model and 2) Identify gaps in knowledge to provide guidelines for future development of constitutive models.
The review strategy is described in Section 2 and the availability of data is discussed in Section 3. Section 4 provides a brief overview of the constitutive modeling process.
Subsequently, a review of passive and active constitutive models of GI tissues is given in 50 Section 5 and Section 6, respectively. In Section 7, the results from our cross-analysis are presented. Section 8 ends the review with a discussion of the open challenges, perspectives, and directions for future studies based on our findings. A summary of the notation used in the manuscript is provided in Appendix A. 55 A PubMed search was performed for each GI organ using the terms listed in Table   S1 (supplemental material, c.f. Section 3). The search was set to Title/Abstract. An according to the SPARC Data Structure using the SPARC data curation software SODA [17][18][19] and the files are made available under open and permissive licenses (CC-BY for data files and MIT for code files). 90

Basis of constitutive modeling
We provide here an overview of the constitutive modeling process. This section is limited to concepts identified during our review of GI tissues although it could be applicable beyond. For more details related to continuum mechanics or constitutive modeling of soft tissues, we refer the readers to the literature available on these topics [20-22]. 95

Strain energy function
It has been shown that soft tissues can endure large deformation and their non-linear load-displacement behavior is best modeled by defining a Helmholtz free strain energy function Ψ such that the stress and strain in the material are related as follows: where σ is the true Cauchy stress tensor, F is the deformation gradient, J = det(F) refers to the Jacobian of the deformation map, and C = F T F designates the right Cauchy-Green deformation tensor. The constitutive modeling process consists of establishing an explicit formulation of Ψ, as explained in the next section. The formulations established for GI 100 tissues are discussed in detail in Sections 5 and 6. Here, we provide an overview of the structure of Ψ.
GI tissues are modeled as homogeneous material, meaning that Ψ is independent of the location on the tissue. This has been deemed as an acceptable assumption for all soft tissues at the scale of interest. Moreover, a decoupled form is assumed for Ψ such that: where U is a purely volumetric (dilatational) contribution and W is a purely isochoric (volume preserving) contribution. The volumetric term U is only considered in a few of the studies reviewed [23][24][25][26][27][28]. In most studies, GI tissues are typically assumed to be incompressible (i.e., their volume remains constant). Interestingly, no GI-specific studies are referenced to support this assumption. It appears to have been extrapolated from studies on other soft tissues. Under the incompressibility assumption, U = 0. Since J represents the change in volume with respect to the reference configuration, we also have J = 1. The incompressibility constraint J − 1 = 0 is typically imposed directly on the energy function as: Here, p serves as an indeterminate Lagrange multiplier to enforce incompressibility. It can be identified as the hydrostatic pressure and may only be determined from the equilibrium equations and the boundary conditions. For GI tissues, the isochoric contribution W is typically split as follows: where W p represents the passive contribution and W a the active contribution. The passive contribution is associated with the mechanical contribution of fibers (elastin, collagen) in the tissue and the ground matrix embedding these fibers while the active 105 contribution is associated with the contraction of the smooth and striated muscle cells and the electrical behavior of pacemaker cells.

Modeling of passive mechanical properties
A majority of the studies reviewed investigated passive behavior only. The passive contribution W p of the strain energy function can be further split as follows: Here, W e represents the hyperelastic part that is associated with the instantaneous response of the tissue and is only dependent on the strain. W v represents the viscous part 110 associated with the rearrangement mechanisms of tissue micro-components over time.

Typical constitutive modeling process
Constitutive modeling in biomechanics consists of identifying a suitable mathematical model/formulation of the strain energy function Ψ for the tissue of interest [29].
Numerous standard formulations cover the diverse range of behavior of soft tissues, as we will summarize in the next sections. Each of these formulations contains unknown parameters, commonly designated as material parameters. If an optimal tuning of these is then only a function of the material parameters of the selected constitutive model.
To estimate these parameters, either the experimentally measured force quantities are converted to corresponding stress values or theoretical formulations are established for these quantities as a function of the Cauchy stress components. As an alternative to theoretical formulations, computational estimations based on the Finite Element Method (FEM) are also used in some studies [30,31]. Then, optimal material parameters of W for the tissue are estimated by fitting the theoretical/computational stress/force values to the corresponding experimental ones for each experimentally imposed/measured deformation. During the fitting process, the following cost function is typically minimized numerically as: where, N indicates the number of experimental data points, s refers to the set of material parameters for the strain energy function of choice, "exp" designates experimental values from mechanical tests, and "mod" indicates the corresponding theoretical/computational value based on the selected constitutive model and associated material parameters. The symbol A represents the quantity that is fitted, which depends on the type of mechanical tests conducted. For instance, A is the axial Cauchy stress σ a if only uniaxial tests were conducted and the cost function is expressed as: Parameters quantifying the goodness of fit (e.g., the final cost function value e(s), the coefficient of determination R 2 , and the reduced chi-square χ) are finally computed to decide if the selected strain energy function is suitable for the tissue. If so, the parameters are commonly validated with test data not used for parameter optimization to confirm 145 the suitability of the strain energy function. If not, the process is repeated with another choice for the constitutive model until a satisfactory model is identified [32,33]. Once a suitable model is identified and the material parameters are estimated, an explicit relation is obtained between the stress and strain in the tissue. The formulation can then be used to analyze the mechanical characteristics of the tissue or integrated into a computational 150 framework, e.g., using the FEM, to conduct load/deformation simulation of the tissue.
In the GI studies reviewed, different constitutive models and constitutive modeling processes (mechanical testing protocol, tissue source, etc.) have been proposed. They are discussed in the subsequent sections. 155 We provide in this section an overview of the most common formulation for passive GI tissues encountered during our review.

Phenomenological models
The first model of a hyperelastic strain energy function applied to a GI tissue was Fung's model was originally introduced to characterize arterial tissues [48][49][50][51][52]. It is developed in terms of the components E ij of the Green-Lagrange strain tensor E = 1 2 (C − I) with C = F T F being the right Cauchy-Green deformation tensor. It is expressed in its most general form as follows: The constant c is a stress-like material parameter while a ijkl are dimensionless material parameters, and the notation r, θ, z refer to the radial, circumferential, and longitudinal directions associated with the tubular tissue of interest. None of these parameters relate to the structure of the tissue, and hence the phenomenological nature of this model. In most studies, the thickness of the tissue is assumed to be negligible compared to the radius and the length leading to a plane stress assumption in the radial direction. Moreover, the tissue behavior is assumed to be orthotropic and shear is not considered. These assumptions lead to the following three exponential parameters formulation [12,36,37,40,41,[43][44][45]47]: Such a model was previously proposed for characterizing the behavior of skin [56], arteries

170
[57, 58], visceral pleura [59], and pericardium [60], among others. The term q, referred to in some studies as the quadratic component, characterizes the tissue's linear behavior (typically at small strains), while the exponential term dictates the non-linear behavior (at larger strains). All three GI studies mentioned above used a three-parameter (b θθθθ , b zzzz , and b θθzz ) quadratic term [34, 40,46]. In one of these studies, it was found that 175 the model with a three-parameter quadratic term and three-parameter exponential Fung term provided a good fit against tubular inflation-extension loading in all the ranges tested for the small intestine tissue while the three and six exponential parameter Fung models without quadratic term only provided a good fit against data in physiological loading ranges [46]. In another study [40], it was shown that use of a three-parameter quadratic term with a three-parameter exponential Fung term did not provide a better fit than using a three-parameter exponential Fung term alone for the mucosa-submucosa and muscle layers of the esophagus. In the same study, the authors also tested a neo-Hookean model (c.f. next paragraph) instead of the quadratic term and reached the same conclusion.

185
Various other phenomenological models have also been proposed to describe the hyperelastic mechanical behavior of intact stomach, large intestine, and rectum tissues as well as the individual layers of the large intestine. A list of such studies is provided in Table   S3 along with relevant information. Contrary to Fung's model, these other phenomenological models are isotropic and thus do not account for the variation of tissue behavior depending on the loading/deformation direction. These models have only been validated against tissue testing results in one direction at a time using compression and/or planar uniaxial tests. Among these, invariant-based phenomenological models have been most commonly utilized. They take the form of a polynomial function in the two invariants I 1 = tr(C) and I 2 = 1 2 tr(C) 2 + tr(C 2 ) of the right Green-Cauchy deformation tensor C as: The constants c ij are the material parameters (with c 00 = 0) with unit of pressure. In the particular case where only c 10 is non-zero, the strain energy function is designated as the neo-Hookean model [61]. In the case where only c 10 and c 01 are non-zero, it is referred to as the Mooney-Rivlin model [62,63]. One study proposed an extended Mooney-Rivlin model for the large intestine tissue where c 11 is also non-zero [64]. One study also used a  Besides invariant-based models, Ogden's model [65], which depends on the principal stretches λ i with i = r, θ, z (λ 2 i = eigen(C)), has been used in various studies [66][67][68][69] as listed in Table S3. The general form of the Ogden's model is expressed as follows Here, µ k (unit of pressure) and α k (dimensionless) are the material parameters of the model. The constant N determines the order of the model. For GI tissues, N = 1 [66,67,69] and N = 3 [68] have been used.

Structure-based models
Most of the structure-based models for GI tissues are inspired by the formulation, An additive split of the strain-energy function is suggested in Holzapfel-type models between a part associated with isotropic deformations and a part associated with anisotropic contribution, as follows The isotropic contribution is associated with the mechanical response of the non-collagenous In one study for the esophagus, no isotropic part was included [76]. Rather, the noncollageneous part was attributed solely to the elastin fibers, and an anisotropic pseudoelastic model was used to describe its behavior. The anisotropic contribution is associated with the mechanical response of collagenous fibers. In its most general form, an exponential function of the following form is used to describe the strain energy stored in the collagen fibers of the tissue: with Here, the superscript j designates a family of collagen fibers, typically defined as all collagen fibers aligned along the same mean direction. The coefficient N indicates the total number of fiber families and n is typically equal to N . The stress-like parameter k (j) 1 > 0 and the dimensionless parameter k is defined as where C designates the right Cauchy-Green deformation tensor, and a (j) represents a unit vector indicating the mean fiber direction. In all the GI studies, it is assumed that the collagen fibers are embedded in the tangential surface of the tissue (no components in the radial direction). Thus, in a cylindrical polar coordinate system, the direction vector is expressed as with e θ and e z referring to the circumferential and axial directions in a cylindrical polar coordinate system, respectively and γ (j) denoting the angle between the a (j) and e θ . In some studies, γ (j) is estimated from microscopy imaging while in others it is considered as an unknown material parameter and estimated from fitting against mechanical testing data. In the original model proposed for arterial wall, a two-fiber family model was proposed with the same stiffness properties (k 1 and k 2 ) and such that they symmetrically oriented (γ (2) = −γ (1) ). Yang and colleagues found that this exact model did not provide a good fit for the mucosa-submucosa and muscle layers of the esophagus 220 against planar uniaxial tests [73]. Looking at their mechanical testing data, they postulated that this model was too nonlinear. Accordingly, they formulated two modified versions of W aniso , which they called modified exponential and bilinear models. They found that their bilinear model provided the best fit for both esophageal layers. In a subsequent study where they considered tubular inflation-extension tests in addition to the 225 planar uniaxial tests, they found that their modified exponential model provided a better fit [74]. Later studies have shown that rather than a modification of the original Holzapfel model, an extension of the model that captures better the collagen fiber distribution in the GI tissues provided a good fit to the tissue's behavior. Overall, multiple fiber family configurations have been tested and validated for GI tissues: One family [28,77], two 230 symmetrical fiber families similar to the original Holzapfel model [82], two-fiber families along the longitudinal and circumferential directions [81], three-fiber families along longitudinal and two symmetrical directions [76,78,79], and four-fiber families along the longitudinal, circumferential and two symmetrical directions [75,76,80]. Typically, the set of parameters k 1 and k 2 is assumed to be equal for the two symmetrical directions.  Table S5. This type of model uses a split between isotropic and anisotropic contribution 240 similar to the Holzapfel-type models formulated in Equation 13. Different formulations are, however, used for both. The model of the isotropic contribution is formulated as a non-linear function of the three invariants of the right Cauchy-Green deformation tensor C. The model for the anisotropic part account for an exponential function of the fourth invariant similar to the Holzapfel-type models but also consider additional components.
We refer the readers to the corresponding studies for more details.
One study has proposed a structure-based model with a different formulation for the esophagus tissue, referred to as Microcontinuum Mechanics model [84]. The proposed model is an extension of the Saint-Venant model for classical elasticity to large deformation with an additional parameter that describes the local orientation of fibers (thus 250 making it a structure-based model per definition). We refer the readers to the corresponding study for more details [84].

Microstructure-based models
Microstructure-based models account for the microstructural contribution to the macroscopic stress as they incorporate not only the orientation but also the mechanical contri- where, N indicates the number of collagen families. Similar to the Holzapfel-type models, 270 a family of fiber is defined as all collagen fibers aligned along a preferred direction. A dispersion of the fibers is, however, considered around this preferred direction as discussed below. For each family j, the micro-sphere approach was used to express the overall strain energy function W f,j as the average of the contribution, over a unit volume, of each fiber of the family. Specifically, integration of the strain energy of individual fibers w f,j over a 275 unit sphere S was considered such that the contribution of fibers dispersed in all direction is accounted as follows: where r expresses the unit vectors associated with the direction along the spherical integration and λ(r) = √ Fr · Fr indicates the stretch of the fibers in that direction. Here, n is a constant representing the isotropic network chain density [87], ρ j is the orientation 280 density function (ODF) to take into account the fiber dispersion. A numerical Gaussian quadrature was used to estimate the integral leading to the following expression of W f,j : where m represents the total number of Gaussian points, and (ω i , r i ) the associated sets of Gaussian weights and points. This led to the following overall formulation of W aniso as: A Holzapfel-type model was proposed for w f,j accounting for the tension-compression switch formulation (i.e., fibers do not support compressive loading). In this study, N = 2 was chosen based on structural observation of the tissue. However, only one set of parameters (k 1 , k 2 , θ) was considered indicating that the authors assumed a symmetrical orientation of the preferred directions for the two col-290 lagen families. Two different ODFs ρ j , respectively known as Microfiber von Mises model and Microfiber Bingham model, were used in the two microstructural models that were tested in the study. We refer the readers to the corresponding study for the formulation of the ODFs and additional details. Overall, the study found that the Holzapfel-type four-fiber families model provided the best fit to biaxial testing data, even in comparison 295 to the miscrostructure-based models tested. The microstructure-based models still provided a good fit (0.83 < R 2 < 0.95 for model fitting across specimens for the different locations considered in the study along the colon) with fewer parameters (5 vs. 8 for Holzapfel-type four-fiber family model).

Viscous and visco-hyperelastic models 300
Overall, three studies investigating the visco-hyperelastic behavior of GI tissues were identified, as summarized in Table S6. We outline below the formulation utilized in these studies. We refer the readers to the respective papers as well as classical literature on the topic [88,89] for more details.

In the first study published in 2007, Higa and colleagues [90] proposed a Mooney-
Rivlin incorporated convolution integral model to characterize the stress relaxation in the colon tissue of the following form: where t is the time measured with respect to a reference time (e.g., beginning of the experiment), the integral variable s represents an infinitesimal time within time 0 to t, and the term under the brackets represents the Mooney Rivlin model presented previously in Equation 11. Here, g represents a normalized relaxation function expressed in terms of the generalized Maxwell's model. In this study, n = 2 Maxwell elements were selected leading to the following expression of g as: where g 1 , g 2 (normalized modulus) and τ 1 , τ 2 (relaxation times) are material parameters.

305
The model parameters were estimated from in vivo compression tests at different velocity and resulting model-estimated stress values provided good agreement with experimental data.

Fontanella and colleagues [83] modeled the visco-hyperelastic behavior of the stomach.
A Natali-type model, presented in section 5. 1.2, was used for the hyperelastic component.
A viscous variable-based model was used to characterize stress relaxation in the tissue as: This model was previously utilized to characterize the viscous behavior of periodontal ligament [91,92] and heel pad region [93], among others. Here,Ċ = ∂C ∂t designates the time derivative of the right Cauchy-Green deformation tensor C, t is the time measured with respect to a reference time (e.g., beginning of the experiment), and the integral variable s represents an infinitesimal time within time 0 to t. The parameters n indicate the number of viscous elements in the model [88,89]. It was assumed in this study that n = 2. The evolution of the viscous variables q is determine by the following differential where γ i and τ i (i = 1, 2) are material parameters of the model. The parameters of the hyperelastic model were estimated using previously published planar uniaxial tests [94]. 310 Viscous parameters were identified from relaxation curves obtained from ex vivo tubular inflation tests at various inflation/deflation rates and relaxation times. The model was subsequently implemented into a FEM framework and model parameters were further fine tuned to validate the predictions of the FE model against the tubular inflation data.
Panda and Buist [95] investigated the visco-hyperelastic behavior of the small intestine, large intestine, and rectum within a framework proposed previously as an extension of the original work of Huber and colleagues [96,97]. A multiplicative decomposition of the deformation gradient was assumed as: where subscripts e and i indicate elastic and inelastic parts of F . The following decomposition of the strain energy function was also introduced as: where C e = F T e F e , and Ψ E and Ψ OE are strain energy functions associated with the inflation data aimed at studying stress relaxation [98] and demonstrated that their model was able to reproduce the experimental data with a neo-Hookean formulation for Ψ E and Ψ OE in conjunction with either linear or non-linear viscous models (it was not specified which performed better). They used colonic inflation data aimed at studying creep [99] and demonstrated that their model was able to reproduce the experimental data with 325 a Humphrey model (we classified as Fung-type) for Ψ E and Ψ OE along with a linear viscous model. They used rectal inflation data aimed at studying the behavior of the rectum under normal conditions and Hirschsprung's disease [100], and demonstrated that their model was able to reproduce the experimental data for healthy and diseased tissue with a neo-Hookean formulation for Ψ E and Ψ OE along with a linear viscous model.

Active models of GI tissues
Active models comprise a large class of biological tissues characterized by subcellular microstructural rearrangement scaling up to whole organ contractility. In the specific case of slow wave gastroenterology activity, such a feature entails a tight interplay between interstitial cells of Cajal (ICCs) and smooth muscle cells (SMCs) (via specialized 335 gap junction proteins), thus implying an intrinsic multi-field coupling that takes into account transmembrane voltage, calcium dynamics, and actine-myosine overlapping [101] (see section S1). The complexity of active tissue modeling is further comprised of the (i) Multiple scales involved both in time and space (from milliseconds to minutes, from microns to centimeters), (ii) Strong directionality due to collagen and SMC reinforcements 340 within the multi-layered GI wall (usually two and three orthotropic reinforcing layers), and (iii) Marked spatial dependence of GI excitability driving slow waves propagation and associated contraction. These overlapping layers can be separated experimentally and studied individually (a procedure not applicable in other organs, including the blood vessels and heart). Therefore, more valid experimental input and validation processes of 345 complex models are foreseen compared to other organs. Accordingly, generalized electromechanical constitutive laws can be formulated requiring the introduction of advanced notions of thermodynamics of continuum media [102] enforcing the need for homogenization techniques within a multiphysics framework [103,104].
A list of studies presenting an active structure-based model is provided in Table S7. 350 We review common formulations of GI tissues' active models according to their historical appearance and mathematical complexity.

Active Stress
The active stress model, initially developed for the cardiac tissue [105], has been typically adopted as a first modeling approximation of the complex excitation-contraction 355 coupling typical of smooth muscle cells, with particular reference to the esophagus [30,106,107]. An example can be found for the stomach [27] and the small intestine [108], while no publications were found for the large intestine and the rectum using an active stress formulation. As discussed in the literature [109][110][111], the active stress approach suffers from a rigorous thermodynamical derivation. Accordingly, advanced approaches 360 are briefly discussed in Sec. 6.2.
The basic idea of the active stress approach is to superimpose an active contribution σ a to the passive (elastic) stress σ p such that the total equilibrium Cauchy stress can be defined as: with the pulled-back counterpart expressed in terms of the first Piola-Kirchhoff stress tensor: In such a way, the constitutive modeling can be framed into an invariant reference system and can be used both in small and large strains frameworks. Accordingly, the passive stress can be described by one of the hyperelastic laws described before, whereas the active part is directly coupled with the electrophysiological variable responsible for tissue contraction, e.g., the membrane voltage dynamics (see next section and section S1). Although the electromechanical coupling is due to complex multiscale intracellular calcium dynamics (known as Calcium-induced, Calcium-released or CICR mechanism), the active stress approach usually adopts a phenomenological volumetric formulation where: Here, the switch function ϵ(V ) together with the material parameter k Ta

Active Strain
The active strain approach, originally developed for the cardiac tissue [109,[112][113][114][115], has been recently applied to the esophagus [82], the stomach [116], and the colon [117] in order to enforce active contractility guided by ICC electrophysiology (see next section and section S1). Adopting the notion of distortion (active deformations [109,118]) in continuum mechanics, the approach is based on the multiplicative decomposition of the deformation gradient (similar to well-known plasticity theories [119]) into a passive elastic, F e , and an active, F a , part, as: In this way, it is possible to enforce multiscale and multiphysics couplings in a continuum homogenized framework through the local active strain deformation map. In particular, intracellular calcium and transmembrane voltage dynamics occur at the cell or subcellular level in much smaller space and time scales with respect to the tissue or 370 organ levels. Accordingly, the local traction-free configuration of the continuum body is herein constitutively related to smooth muscle contractility in terms of infinitesimal volume elements. A representative scheme of such a mapping is provided in Fig. 3, in which an intermediate non-compatible configuration arises due to the presence of the active map F a .

375
From a computational point of view, once the active map is known because of muscle contractility, one can compute the elastic deformation gradient by imposing both the balance of linear momentum and geometric compatibility. Moreover, as usual for soft biological media, the tissue is considered incompressible, i.e., J = detF = 1, and further assuming J a = detF a = 1 it follows that J e = detF e = 1. Finally, in view of obtaining frame invariant constitutive formulations, the right Cauchy-Green deformation tensors can be introduced from (34) as: and the total first and second Piola-Kirchhoff equilibrium stress tensors read, respectively; as: In the previous equation, P e , S e represents the classical elastic first and second Piola-Kirchhoff stress tensors, respectively.
The active strain approach constitutive law is completed by the introduction of the specific formulation for the active map F a . As detailed in [116], considering longitudinal and circumferential principal directions of GI muscle contractility, we can write: where N c , N l are the orthogonal unit vectors in the circumferential and fiber direction, respectively, defined in the reference configuration and within the plane encompassing SMCs, while N n = N c × N l represents the orthogonal to the plane thickness direction.
In Equation (37), α c , α l represent material parameters ruling the amount of contractility in the corresponding directions 1 , while γ n enforces incompressibility as: Finally, the excitation function γ(V ) couples the mechanical problem with the electrophysiological one (see next section), in particular via a smooth activation function defined in [116] in terms of the transmembrane voltage of SM as: where β 1 , β 2 , V thr are material parameters linked to the intracellular Ca2+ dynamics, while H(V − V thr ) is a Heaviside step function switching on the active contraction whenever the threshold V thr is reached.

Active Electromechanics
Active stress and active strain approaches have been recently collected under a sound thermodynamical framework [120] generalizing the constitutive modeling of active media electromechanics in the case of the small intestine [121,122]. To this end, Maxwell electrostatics theory has been introduced, dealing with the material electric field E = −∇ X V defined as the gradient of the electric potential (i.e., the transmembrane voltage V ), the material electric induction D, and the polarization tensors Π linked through the following constitutive relation: where ϵ 0 is the vacuum dielectric constant and Π accounts for electric distorsions due to materials deformations. Following previous works [120,[123][124][125], the local thermodynamics state is defined through an electrical Helmholtz free energy A characterized by the functional dependence on F, E and a collection of internal variables Q, such that: where P is the equilibrium stress and Y the thermodynamic forces work-conjugate to Q. According to classical theories of multiphysics couplings in finite elasticity, an additive decomposition of the free energy is assumed in conjunction with the multiplicative decomposition of F (see Eq. (34)), such that: and the associated stress assumes an additive composition similar to the active stress approach, (31) and (36). In other words, P a assumes an explicit expression when a constitutive prescription of F a is provided together with the active inelastic potential.
In this regard, the directional active Helmholtz free energy can be decomposed into isotropic and anisotropic contributions of the GI wall based on the functional dependencies on the electric field E, the fiber direction n, and associated structure tensor G = n ⊗ n (by means of the fourth invariant I 4 (n) = C : G) as: Respecting the usual linear isochoric and quadratic directional anisotropic contributions (see original paper [120]), the active energy formulations state: with χ iso , χ aniso the associated permittivity tensors. Accordingly, the final expression of the active stress, in terms of the second Piola-Kirchhoff stress tensor, is: Visco-electro-elastic distributed models. In Refs. [122,124,125], the active electromechanics setting was further extended to include viscous deformation and stochastic material properties [126,127] to model the colon wall. Assuming a generalized multiplicative decomposition of the deformation gradient (see Fig. 3): the viscous deformation can be considered as an internal variable, F v = Q, governed by suitable kinetics. Accordingly, the Helmholtz free energy density becomes: from which the total first Piola-Kirchhoff stress tensor reads: where P E = P p +P a represents the equilibrium stress, e.g., from (41), and the constitutive prescription of the viscous stress is provided by a suitable pseudo potential Ψ * as: Furthermore, the active deformation gradient, F a , incorporates directional-dependent stochastic material properties, as key features of GI wall, by means of the functional dependence on a distributed angular variable β, i.e.: Here, β refers to the angle between the aleatoric fiber direction, a, and the electric field unit vector m: such that the expected value (integral average) of the active deformation gradient becomes: according to the usual structural tensor approach. Here, H is the average fourth order structural tensor associated with the second order approximation worked out for the passive contribution [126,128,129]: 6. 3

.1. Multiphysics mathematical modeling of GI tissues
Few attempts have been made to study multiphysics couplings in GI tissue. An excellent example is provided in [117,130], modeling the thermal coupling within the intestine tissue because of its critical role in paralytic ileus disease. Adopting the approach proposed by Bini et al. [131], intestinal excitability was coupled to the thermal transport problem by using the Pennes bio-heat equation. Accounting for the ability of tissues to remove heat both by passive conduction (diffusion) and blood perfusion, these contributions were combined in a generalized reaction-diffusion PDE of the form: We refer to the original work [130] for details on model formulation and implementation.
The role of ion channels mechanosensitivity [137] and recent advances in the electrophysiological mathematical modeling of GI is referred to recent reviews on the subject [101].
We present a cross-analysis to provide an overview of the existing studies for modeling GI tissue mechanics and to help identify remaining gaps.

Overview
A total of 54 studies were found to fit our review criteria (when counting individually for each organ). Among them, 40 (80%) investigated a hyperelastic model, 5 (9%)

Species
Tissue source for mechanical testing was clearly identifiable in 50 studies. Seven Among those studies, 6 studies investigated the effect of a disease on the mechanical behavior of GI tissue. Five of the studies were on the small intestine and investigated 420 the hyperelastic behavior of the tissue in diabetic rats [37], active behavior in patients with systemic sclerosis [108], hyperelastic behavior in obstructed guinea pigs [42], hyperelastic behavior in obstructed rats [47], and visco-hyperelastic behavior in patients with Hirschsprung's disease [97]. One study investigated the hyperelastic behavior of the esophagus tissue in diabetic rats [38]. 425

Mechanical testing condition
The mechanical testing condition and protocol used to conduct the testing necessary for estimating the material parameters of the constitutive models and/or validating the models were identifiable in 50 studies (7). In 48 studies (96%), the model parameters were estimated based on testing conducted in ex vivo conditions. A total of 14 different

Studies on human tissues
An overview of the studies proposing a model based on human tissue testing is provided in Table S8. A total of six such studies have been identified in our review. They proposed models for the intact small intestine [95,108], large intestine [63], and rectum tissue [62,67]  Only three manuscripts have investigated modeling of the visco-hyperelastic behavior of the GI tissues. They studied the stomach, large intestine, and rectum. Given the scarcity of such models, selection of a model is simply guided by availability for the tissue of interest.
Although the critical role of smooth contractile cells into the modeling part is well-495 established, studies on active models development for GI tissues are still very limited. The complexity of the phenomenon and its intrinsic multiscale nature require further effort to produce reliable and robust modeling tools for device design and optimization.

Suggestions for future studies
Microstructure-based hyperelastic models are now commonplace for cardiovascular 500 tissues [140,141]. To date, only one such model has been established (for the large intestine). It is thus recommended that future hyperelastic studies investigate such models.
Visco-hyperelastic models are overall scarce, and especially lacking for the esophagus.
Active models have attracted more attention in recent years, but none have been vali-

Closing comments
Notwithstanding the substantial importance of GI mechanics in our everyday life, the 525 research and industry efforts do not compare with the level of modeling, analysis, and investment associated with cardiovascular mechanics. The development of a comprehensive constitutive and computational model of GI mechanics, parametrized for each different section, can lead to patient-specific in-silico tools unveiling the principles governing food digestion and disease mechanisms at different scales. Therefore, such a paradigm stands 530 as an unprecedented opportunity for computer-aided medicine, innovative device design, and optimal treatment approaches. Although modeling approaches can be adapted from the cardiovascular field, GI-specific technologies are needed to gather mechanical behavior data needed to inform and validate such models, especially in humans.

535
We thank funding from NIH SPARC OT2OD025308 and OT2OD028203.

Author contributions
All authors have participated sufficiently in the work. B. Patel and A. Gizzi contributed to conception and design of the study, acquisition/analysis/interpretation of data, and drafing/revising the manuscript. H. Gregersen and G.S. Kassab contributed to concept 540 and design of the study, interpretation of data, and revising the manuscript.
smooth boundary ∂Ω, ∂Ω s and outward unit normal N, n, respectively. A material point in Ω, Ω s is denoted by X, x, respectively, with motion x = x(X, t). Material and spatial gradient, divergence, and determinant operators are identified with ∇ X , ∇ x , ∇ X · , ∇ x · , and det X , det x , respectively. Cross, vector, and tensor products are denoted ·, ×, ⊗. If not explicitly stated, we denote second order tensors in matrix notation with capital bold