Virgin Passive Colon Biomechanics and a Literature Review of Active Contraction Constitutive Models

: The objective of this paper is to present our ﬁndings on the biomechanical aspects of the virgin passive anisotropic hyperelasticity of the porcine colon based on equibiaxial tensile experiments. Firstly, the characterization of the intestine tissues is discussed for a nearly incompressible hyperelastic ﬁber-reinforced Holzapfel–Gasser–Ogden constitutive model in virgin passive loading conditions. The stability of the evaluated material parameters is checked for the polyconvexity of the adopted strain energy function using positive eigenvalue constraints of the Hessian matrix with MATLAB. The constitutive material description of the intestine with two collagen ﬁbers in the submucosal and muscular layer each has been implemented in the FORTRAN platform of the commercial ﬁnite element software LS-DYNA, and two equibiaxial tensile simulations are presented to validate the results with the optical strain images obtained from the experiments. Furthermore, this paper also reviews the existing models of the active smooth muscle cells, but these models have not been computationally studied here. The review part shows that the constitutive models originally developed for the active contraction of skeletal muscle based on Hill’s three-element model, Murphy’s four-state cross-bridge chemical kinetic model and Huxley’s sliding-ﬁlament hypothesis, which are mainly used for arteries, are appropriate for numerical contraction numerical analysis of the large intestine.

Since over 10% of overall deaths in the US each year occur due to digestive tractrelated disease only, treatment failures could potentially be due to inadequately studied bowel mechanics, which has been well reviewed by Bhattarai et al. [1]. Not only endoscopic examination, which may injure the vulnerable tissue of the digestive tract, but also excessive compression of the tissue during a surgical procedure, such as stabilization of anastomosis, may increase the risk of patients developing ischemia (vascular damage that restricts oxygen supply to the blood) and subsequently necrosis (cell death) due to mechanical trauma. Since 25% of the total cardiac output is received by the GI tract alone [35], a drop of blood flow by more than 50% induces ischemia. In addition, progressive tissue remodeling may result in abnormal intestinal motility, which may aggravate chronic diarrhea and the overall nutrient absorption process. Most importantly, a healthy mechanical state of the intestine not only maintains food digestion but also allows effective communication with the brain for drug administration and water absorption, avoiding pathophysiological diarrhea, which consequently regulates the patient's wellbeing.
A comprehensive study of the human large intestine decoupled into extensive experimental and computational approach further supported by proper constitutive approximations predicts its relevant state-of-the-art during physiological and non-physiological processes. However, it is ethically impossible to obtain sufficient human intestine samples for statistical studies. Nevertheless, histological findings suggest the structural (inner mucosal, submucosal, circumferential muscular, longitudinal muscular and outer serosal layers), mechanical (progressive reduction in stiffness from proximal towards distal large intestine) and functional similarities of the human large intestine with that of the pig, subjects that are commonly used for the pre-clinical trials of medicine, surgical procedures and medical devices [1,17,[36][37][38].
The aim of this paper is (a) to perform equibiaxial tensile experiments on porcine large intestine, (b) to identify the material parameters for an anisotropic hyperelastic model for describing its passive mechanical behavior and (c) to develop numerical tools for the numerical validation of the estimated material parameters. This paper further provides a comprehensive review of our current understanding of the active colorectal contractile components that are equally important for describing the behavior of the intestine during physiological processes.

Preparation of the Equibiaxial Tensile Test Specimens
The large intestine of the pig is removed as a whole organ from the euthanized experimental animals previously used in medical training at the Institute of Laboratory Animal Science at the University Hospital RWTH Aachen and covered with cloths soaked in 0.9% NaCl solution to prevent the tissue from drying out. Within an hour, the lump fresh intestine is brought to the testing facility and divided into segments based on the individual anatomical loops and completely cleaned of feces and other connective/adipose tissue. Each segment is taken in turn and cut lengthwise. Individual samples are cut out with an edge length (L) of 40 mm × 40 mm and stored in a 0.9% NaCl solution preventing desiccation until mechanical testing. The average thickness (0.68 < T < 2.15) of the large intestine is equal to 1.17 mm measured with a micrometer screw prior to testing.

Test Protocol
The equibiaxial mechanical tests are carried out on a self-developed biaxial tensile testing machine "BiAiX" (FH Aachen, Aachen, Germany, patent pending, DE 10 2017 116 067 A1). For 2D planar evaluation, one Allied Vision Prosilica GT2450, 5 Megapixel camera (Allied Vision Technologies GmbH, Statroda, Germany) with GigE Vision Gigabit Ethernet interface and a maximum frame rate of 15 f/s with AQUAMARINE 2.0/28 C (Schneider Kreuznach, Bad Kreuznach, Germany) lens, is mounted above the BiAiX machine. The two force sensors of the BiAiX tension test machine are calibrated in an industrial manner with calibration weights. The camera is positioned perpendicular to the specimen surface. This guarantees an easy calibration for the 2D DIC with a specified length in the measuring plane. The specimens are then mounted in the apparatus with five needles on each side, taking care to ensure that the orientation of the sample always the same: longitudinal direction along the L-L axis and circumferential direction along the C-C axis (Figure 1a). The needles can move laterally to allow free stretching of each specimen's side. To ensure a uniform evaluation of the specimens, the inner mucosal layer of the specimens is always oriented facing up and a graphite pattern is sprayed on (Figure 1b) for digital image correlation (DIC). Polarized light with polarizing filters in front of the camera setup was used to avoid any reflections on the specimen. The measurement protocols for the tensile test are as follows: 10 cycles of preconditioning followed by a monotonous stretching both at the speed of 10 mm/min until failure and a framerate of 1 Hz (Figure 1c).

Virgin Passive Anisotropic Strain Energy Function Description
The strain energy function (SEF) (ψ) of the incompressible anisotropic large intestine can be decomposed into volumetric (ψ vol ), isochoric isotropic (ψ iso ) and isochoric anisotropic (ψ aniso ) terms [1,24]: where the Jacobian J = detF is the determinant of the deformation gradient F, the isochoric deformation gradient is F = J −1/3 F,Ī 1 = trC is the first principal invariant of the isochoric right Cauchy-Green deformation tensor, (C = F T F), a 0i = a 0ix , a 0iy , a 0iz T is the unit directional vector of the ith fiber family in the undeformed configuration andĪ 4,a 0i = C : A 0i = a 0i × C · a 0i includes the anisotropic invariants, which are quantitatively equal to the square of the respective fiber stretches λ 2 a 0i . At low strain, collagen fibers are not stretched and do not store strain energy. Thus, the tissue mechanical behavior is solely due to the non-collagenous ground matrix and its isotropic response can be defined using the classical Neo-Hookean material model. At large strains, the collagen fibers resist the applied load and show a strong stiffening effect, which can be represented mathematically by an exponential function using the fiber-reinforced hyperelastic Holzapfel-Gasser-Ogden (HGO) model [39]: where κ 0 > 0 and µ 0 > 0 are the initial bulk and shear moduli, respectively, and have a unit of stress. k 1,a 0i > 0 includes stress-like material parameters and k 2,a 0i involves dimensionless parameters associated with the ith fiber family. It is important to select appropriate k 1 and k 2 such that the anisotropic isochoric strain energy function induced by the collagen fibers does not influence the overall tissue mechanics at the low strain regime.

Numerical Simulations
The finite element (FE) analysis has also been performed in the commercial FE software LS-DYNA (Ansys, Inc., Canonsburg, PA, USA) Version: smp d R11.1.0 in 64 bit Windows 10 to validate the estimated material parameters. Although the standard "MAT-295" material model in LS-DYNA is based entirely on the HGO material [39], it is restricted to only three fiber families and is, therefore, not suitable for modeling the intestinal tissue, which is proposed to consist of four collagen fibers (Equation (2)). Therefore, a user-defined FORTRAN material subroutine has been coded inside LS-DYNA.
Two examples are presented in this paper to demonstrate the numerical performance of the implemented user-defined FORTRAN material subroutine using the material parameters estimated from MATLAB's curve fit. Firstly, a cube of edge length 1 mm meshed with one linear brick element and reinforced by four fiber families (Figure 2a) is considered to validate the basic anisotropic hyperelastic mechanisms of a specific intestine specimen derived from the experiment. The four fiber families are represented by unit directional vectors, a 01 = {1 0 0} T , a 02 = {0 1 0} T , a 03 = {cos30 • sin30 • 0} T and a 04 = {cos30 • − sin30 • 0} T . For equibiaxial tension, the boundary conditions are chosen in such a manner that the four bottom face nodes are constrained to move along the Z-axis, while displacements ∆X = ∆Y = 0.3 mm are prescribed on all nodes in the X-Y plane. The fitted material parameters of a specific specimen is to be used. To ensure a nearly incompressible intestine, a large bulk modulus, κ 0 = 5 × 10 4 MPa is taken, which provides the initial Poisson's ratio ν 0 = (3κ 0 − 2µ 0 )/(2(3κ 0 + µ 0 )) = 0.49999999 ≈ 0.5.  The second example models the biaxial deformation state of an experimented specimen (No. 3, Table 1) with dimension 40 mm × 40 mm × 1.25 mm. Due to the symmetry in geometry and boundary conditions, only a quarter of the model, 20 mm × 20 mm × 1.25 mm, is considered for the numerical analysis ( Figure 2b). Compared to the other two dimensions, the intestinal thickness is very small; therefore, only three linear brick elements are taken across the thickness, while thirty linear brick elements are considered along the length and width, maintaining the aspect ratio of 1.6:1 in the undeformed configuration. A total of 2700 linear brick elements are used in this example. The linear brick element is used in this example as well to maintain a consistency in comparing the numerical results in both examples for the same material parameters and the same FORTRAN user subroutine. For the biaxial tension, the bottom face nodes are constrained to move in the Z-axis. The symmetry boundary conditions are imposed to the left face and front face nodes while prescribed displacements are applied to the right face and back face nodes until a desired maximum values.

Equibiaxial Tensile Experiment
Displacement-controlled equibiaxial extension has been applied, which took around 12 min to rupture each specimen. The tractions (P i ) are continuously recorded on actuators of longitudinal and circumferential directions at each load step ranging to a maximum value of 8.65-13.75 N and 6.6-15.2 N, respectively. For the DIC evaluation with the ISTRA4D software, residuum = 40-50, facet size > 50 and grid spacing = 1 3 facet size ≈ 17 are used to obtain a homogeneous strain field without a significant loss of grid points. Most importantly, inhomogenous and anisotropic strain distributions are observed, defined by several in-plane ripples that looks similarly to calm ocean waves directed vertically for the X-axis and horizontally for the Y-axis (Figure 3b,c). No concrete evidence can be found where the underlying collagen fibers withstood the stress. However, a higher magnification and higher resolution camera, if available, may be utilized to resolve such strain distributions.

Intestinal Anisotropy
As shown in Figure 4, the equibiaxially stretched porcine large intestine is stiffer in the longitudinal direction than in the circumferential direction: the additional stiff-ness is provided by thicker submucosal collagen fibers (a 03 , a 04 ). Feng et al. [40] used the second harmonic generation and measured the diameters of the load bearing collagen fibers in the large intestine. The largest (≈6.5 µm) is found in the submucosal layer and the smallest, ≈1.9 µm, is found in the muscularis layer. The non-interacting collagen fibers in the longitudinal and the circumferential muscular layers are orthogonal, i.e., α Mus = cos −1 (a 01 · a 02 ) = 90 • . However, the submucosal collagen fibers are oriented crosswise at an angle of α Sm = cos −1 (a 03 · a 04 ) ≈ ±30 • to the longitudinal direction, ignoring the possible interactions between these fibers in this study. Assuming that the collagen fiber diameters in the longitudinal and the circumferential muscular layer are the same, one isotropic (µ 0 ) and four anisotropic material parameters (k i,a 01 = k ia 02 = k i,Mus and k i,a 03 = k i,a 04 = k i,Sm , ∀i = 1, 2) are to be estimated to fit the experimental stress-strain curves. For this systematic characterization of the virgin passive anisotropy of the porcine large intestine (Equation (2)), only five parameters are sufficient to fit the experimental Cauchy stress-stretch data. For this, the Green-Lagrange strains (E ii ) in both loading directions are evaluated directly from the heat maps ( Figure 3) using the DIC technique in the ISTRA 4D (Limess Messtechnik und Software GmbH, Krefeld, Germany) software for each individual images and are recorded as an averaged value from a chosen deformation region ( Figure 3). From the Lagrangian strain and using the incompressibility condition, the true or Cauchy stress (σ ii = P i λ i /LT) is set as a function of the elongation (λ i = √ 2E ii + 1) and curve fitting is performed in MATLAB (The MathWorks, Inc., Natick, MA, USA) on the basis of a non-linear least-squares optimization method to determine the passive incompressible anisotropic material parameters per SEF described in Section 2.3. The fitted curves are shown in Figure 4, for which the estimated parameters are listed in Table 1. The quality of the fitting is measured by calculating the coefficient of determination R 2 ∈ [0, 1] in the longitudinal (R 2 L ) and circumferential (R 2 C ) directions, respectively, using the following: where n is the number of measured data points, σ exp i is the Cauchy stress measured experimentally, σ fit i is the corresponding stress values predicted by the fitting procedure using the SEF and σ mean is the mean of the experimental stress. The higher the value of R 2 , the better the fit is to predict the experimental data.

Parameter Stability
The estimated virgin passive anisotropic parameters (Table 1) are to be fed in the FE software for the numerical simulation and, thus, need to be numerically stable, fulfilling the polyconvexity of the SEF. Unstable parameters and non-convex SEF may lead to unphysical behavior during large deformations [41]. Convexity of the SEF in Equation (2) guarantees the existence of at least one local minimum energy state for which the eigen- 2 12 ] of the planar Hessian matrix (H ij = ∂ 2 ψ virgin passive /∂λ i ∂λ j , ∀i, j = 1, 2) need to be positive [42,43]. Positive eigenvalues render the Hessian matrix positive definite and the SEF convex at every point in the strain space results in a positive definite elasticity tensor or a positive definite stiffness matrix that guarantees stable material parameters. For thin specimens with homogeneous plane biaxial stress and deformation states, only two components are sufficient because the third principal stress is zero. The numerical stability of the fitted parameters for the virgin passive deformation state ( Figure 5) can be verified by plotting the positive three dimensional surfaces of the eigenvalues (η 1 , η 2 ) as the function of biaxial stretch (λ L , λ C ) (see Figure 5).   Table 1.

Equibiaxial Tension of a Linear Brick Element
The equibiaxial stretch simulation result of the linear brick element is shown in Figure 6. Due to the prescribed elongations (∆X = ∆Y = 0.3 mm) along the X and the Y axis, the brick model is compressed (∆Z = −0.61 mm) along the Z axis (Figure 6a). The Cauchy stress and the strain data in the linear element along the longitudinal (X-axis) and the circumferential (Y-axis) directions are extracted from LS-PrePost. The respective stretches in both directions are computed and the Cauchy stress vs. stretch curves are plotted over the fit dataset. From the simulation, a clear anisotropy has been achieved with a longitudinal stiffness larger than the circumferential direction. Furthermore, the material behavior is exactly same as the dataset that characterizes the chosen specimen No. 3 during the equibiaxial tensile experiment. This would confirm that the (a) materials with the estimated material parameters are numerically stable and (b) the implemented FORTRAN user material subroutine is capable for further numerical investigations; as presented in Section 3.4.2.  Figure 7c,d. The contours of the displacements from the FE simulation are quantitatively and qualitatively mapable to that observed in the DIC images from the experiment (Figure 7a,b).
Thus, the implemented user subroutine tool is satisfactorily able to predict the virgin passive mechanical behavior of the multi-layered intestine tissue and similar architecture, provided that the estimated material parameters for the adopted HGO model are numerically stable. However, the physiological response of the in vivo intestine is not limited to virgin passive hyperelastic anisotropy. Intestinal peristalsis due to continuous smooth muscle contraction-relaxation, viscoelasticity due to water content and damage either due to tissue remodeling, growth or injuries are other important aspects of colon biomechanics. Since the damage modeling intestine tissues has not been studied before, the authors would like to review some of the existing literature on (a) the concept of intestinal peristalsis and (b) various constitutive models for smooth muscle contraction and the viscoelasticity of the colorectal tissues.

Physiology of Intestinal Peristalsis
The colonic motility or peristalsis is characterized by the coordinated contractionrelaxation of longitudinal and circular smooth muscle cells (SMCs) regulated by the enteric and the parasympathetic nervous system [44][45][46]. On an anatomical level, during peristaltic propulsion, the circular SMCs behind fecal content contract and relax in the receiving segment, whereas longitudinal SMCs relax behind the luminal content and contract in the front (Figure 8). The luminal cross-section in the receiving segment widens to facilitate the forward-movement of the intraluminal content [47]. On the electrical level, the smooth muscle contraction-relaxation is driven by an excitatory-inhibitory response of an electrically coupled multi-cellular SMC/ICC/PDGFRα + cell syncytium (SIP syncytium) connected via gap junctions [48]. Innervated with excitatory enteric motor neurons, the interstitial cells of Cajal (ICC) stimulate unique Ca ++ -activated Cl − channels or K + channels and generate excitatory electrical slow waves or currents to enhance SIP syncytium excitability and to initiate the rapid contraction of SMCs [49,50]. In contrast, platelet-derived growth factor receptor α-positive (PDGFRα + ), the inhibitory neurotransmitter released from enteric motor neurons, expresses small conductance calcium-activated K + channel 3 (Ca ++ -activated SK3) channel to generate inhibitory currents, thereby reducing SIP syncytium excitability and reducing SMC contraction [49,51]. Most importantly, the level of free cytoplasmic or intracellular Ca ++ within SMC determines its cellular contraction-relaxation activity. In response to a slow influx of Ca ++ through leak channels, acetylcholine (AcH) neurotransmitter, H + , hormones, sarcoplasmic reticulum calcium spark, the stimulation of voltage-gated Ca ++ channels, ATPase, sacro/endoplasmic reticulum Ca ++ ATPase (SERCA), stromal interaction molecule (STIM) protein and Orai protein activities rapidly increase Ca ++ concentrations. Such processes are schematically shown by three greendashed regions in Figure 9a. In any such cases, enhanced SMC contraction is depicted by a rapid spike of the membrane potential, known as action potential well above the threshold potential, which is counteracted by repolarization, i.e., outflow of K + ions from the activated K + channels (Figure 9a).   Figure 9. Schematic representation of processes involved in SMC contraction within the intracellular space. (a) Intracellular Ca ++ is increased mainly by ion channels and excitation by ANS. Green and red arrows represent ion inflow and outflow, respectively. (b) Increased four intracellular Ca ++ ions bind with each inactive CaM (brown) to form a CaM-4Ca ++ complex (green) that stimulates MLCK for MLC 20 activation. (c) ATP and MLCK activity initiates cross-bridging between actin and myosin filaments, followed by power stroke and contraction (−δ). Intracellular Ca ++ from the CaM-4Ca ++ complex is supplied back to the sarcoplasmic reticulum and extracellular matrix via ion channels. Green circles with +ve sign represent an excitatory response, whereas red circles with −ve sign represent an inhibitory response.
On the molecular level, smooth muscle contraction-relaxation is regulated by the opposing activity of the myosin light chain kinase (MLCK) and myosin light chain phosphatase (MLCP), thereby modulating the phosphorylation of 20 kDa regulatory light chain of myosin II (MLC20) [52,53]:The latter is Ca ++ independent. Triggered by inositol 1,4,5-trisphosphate (IP3), Ca ++ activates the enzymatic domain of MLCK, which is mediated by the calmodulin (CaM) via CaM-4Ca ++ -dependent complex formation (Figure 9b). This CaM-4Ca ++ complex phosphorylates the Ser19 hydroxyl group at the neck of MLC20 to stimulate myosin Adenosine tri-phosphate (ATP) hydrolysis into Adenosine di-phosphate (ADP), inorganic phosphate (P i ) and energy at its globular head. The inhibition of tropomyosin (Trop), calponin (CaP) and caldesmon (CaD) followed by the release of P i from the myosin head exposes the myosin binding sites to initiate myosinactin cross bridging (Figure 9c). Powered by the chemical energy stored in the myosin's head, actin-myosin slide over one another with the release of ADP to shorten the length between two dense bodies by −δ. Sliding continues until a new ATP molecule is attached to the myosin head, which immediately detaches the myosin from the actin binding sites. At the end of this cycle, Ca ++ -independent MLCP dephosphorylates MLC20 [52] and inhibits CaM-4Ca ++ -activated ATPase activity [53]. Furthermore, the kinases from cyclic nucleotides inhibit incoming Ca ++ channel activities, phosphorylate MLCK to reduce its affinity for the activated CaM-4Ca ++ complex and finally decrease myosin phosphorylation [54]. Furthermore, they may stimulate Ca ++ uptake from the CaM-4Ca ++ complex and the cytosol to the sarcoplasmic reticulum and outside the SMC through Ca ++ /Na + or Ca ++ /H + exchanges (Figure 9c). Thus, reduced intracellular CaM-4Ca ++ MLCK activity and increased MLCP activity relaxes SMCs.

Constitutive Modeling of Large Intestinal SMC Contraction
In most of the available mechanical characterizations and computational studies of the intestine tissue, the tonic active contraction of the SMCs has been ignored [1]. The mechanical description of the contraction of SMCs qualitatively resembles skeletal muscles [55] and cardiac muscles [56,57]; however, SMCs differ from skeletal muscles in having a longer contraction time at a much slower rate. Myosin ATP hydrolysis provides chemical energy to produce cyclic cross-bridges between actin and myosin filaments and to generate a power stroke to result tonic and phasic SMC (Figure 9). HE Huxley, the nobel laureates AV Hill and AF Huxley and their co-authors pioneered the conceptualization of skeletal muscle contraction, and this formulation is applicable to SMCs as well [58].
The existing muscle models are conceptually based on the three-element muscle model [59] and the sliding-filament hypothesis [60][61][62]. SMC mechanics have been estimated by a system of ordinary differential equations (ODEs) based on the electrical, chemical and mechanical processes involved. The electrochemical process involves a regulation of the intracellular Ca ++ concentration due to stimulated nerve action potential for which the conductance-based model by Hodgkin-Huxley is useful [63,64]. The voltagedependent K + and Na + ion channels (action potential) and leakage (L) ion channel (small conductance for resting membrane potential) are regulated by small gates in the intestinal SMCs. Membrane potentials higher than the threshold potential open the Na + channel to depolarize the membrane and to activate voltage-gated Ca ++ channel for rapid inflow, followed by a rapid spike or action potential and the contraction. In a voltage clamp set up (Figure 10a), the Hodgkin-Huxley model uses electric circuit relationships and ODEs to measure ion currents (I i ) through each ion channels (i = K + , Na + , L) and compares the input/output currents at a desired membrane potential (V m ) at any time t: where I m is the input current flowing into the membrane; I C is the capacity current associated with the membrane capacitor C m that stores the imbalanced positive-negative membrane current; V i = V m and g i are the ith ion channel reversal potentials and respective maximum values of conductance per unit area; n, m and h ∈ [0, 1] are dimensionless gating variables that are associated with K + channel activation (gate opened), Na + channel activation (gate opened) and Na + channel inactivation (gate closed) with the powers in Equation (5) representing the number of gates used. α i and β i are the voltage-dependent transition rate constants for the ith ion channels, where α i is the number of times per second that a closed gate opens and β i is the number of times per second that a opened gate closes. These rates for stretch-sensitive colorectal afferent endings are expressed by Feng et al. [64]. Hai and Murphy conducted chemomechanical setups and proposed a four-state crossbridge chemical kinetic model [65] depending on the chemical states of myosin (Figure 10b): actin detached unphosphorylated myosin (A∼M), actin detached phosphorylated myosin (A∼M p ), actin attached phosphorylated myosin (AM p ), also known as cross bridge, and actin attached unphosphorylated myosin (AM), also known as latch bridge. The model is governed by a first-order differential equations system: with an additional constraint on the population fraction (n i ) of myosin heads in each state, n 1 (A∼M) + n 2 (A∼M p ) + n 3 (AM p ) + n 4 (AM) = 1. This means that the total sum of the myosin fraction must remain constant. Here, (k 1 , k 6 ) are the position-independent state transition rates of myosin phosphorylation, (k 2 , k 5 ) are the position-independent rates of myosin dephosphorylation, (k 4 , k 7 ) and (k 3 ) are the position-dependent rates for myosin detachment and attachment, respectively. The interesting aspects of this models are as follows: (a) k 1 and k 6 are regulated by Ca ++ concentrations and (b) AM→A∼M is irreversible. Furthermore, a hyperbolic dependency of steady-state active stress on cross-bridge phosphorylation [65] and a direct correlation of myosin phosphorylation relative to the force generation and tissue shortening velocity have been found [66]. Singer et al. [67] further supported the Ca ++ -dependent latch bridged mechanical differences with phosphorylated cross-bridges and decomposed the total SMCs stress into active and passive stress: the latter by depleting extracellular Ca ++ . Kato et al. [68] followed Hai and Murphy and proposed a kinetic model to correlate the CaM-4Ca ++ -dependent activation of MLCK, the phosphorylation-dephosphorylation of myosin and isometric contractions in SMC at a steady state. Based on the four-state cross-bridge chemical kinetic model by Hai and Murphy, the active contraction of the guinea pig taenia coli SMC under large deformation has been proposed by Kroon [69]. The calcium-regulated active chemical response is due to the contractile apparatus consisting of elastic spring (cross-bridges stiffness) and a dashpot (relative sliding between actin and myosin) in parallel with the passive viscoelastic combination and is mathematically expressed as follows: where ψ act is the active part of SEF, µ f is related to the stiffness of the cross-bridges and the number of myosin and actin filaments per unit area, (n 3 + n 4 ) is the fraction of crossbridges contributing to the filament stiffness, M is the direction vector of the contractile filaments in the undeformed configuration, C is the right Cauchy-Green deformation, λ e = λ f /λ fc is the elastic stretching of the cross-bridges, λ fc is the relative sliding between contractile filaments and λ 2 f = MCM is the tissue stretch along the contractile filament direction. The evolution of λ fc is adopted as ηλ fc = P cc − tialψ act /tialλ fc , where η is the viscous damping coefficient associated with the cycling/contraction process and P cc is the active stress due to attached cross-bridges causing muscle contraction. The process of the filament contraction development in AM p AM states can be accounted for by an active stress formulation: where κ 3 and κ 4 are the material parameters related to the driving force of cross-bridges in state AM p and passive resisting strength of cross-bridges in state AM, respectively. To model taenia coli SMCs, µ f = κ 3 = 1.11 MPa, κ 4 = 0 (for monotonic contraction with inactive latch state) and η = 44.5 MPa, k 1 = k 6 = 0.17 s −1 , k 2 = k 5 = 0.5 s −1 , k 3 = 4k 4 , k 4 = 0.11 s −1 and k 7 = 0.01 s −1 have been used for both isotonic and isometric behavior with an increase in Ca ++ = 2.5 mol/m 3 [69].
Likewise, an anisotropic visco-hyperelastic formulation has been proposed by Gizzi and Pandolfi [70] to model the colon visco-electro-active response and the interaction of slow waves (electrical control activity) and fast waves (electrical response activity). The electrodynamics of the intestinal longitudinal muscular (LM) and the ICC (I) layers, i.e., dimensionless transmembrane potentials (u LM , u I ) and slow transmembrane currents (v LM , v I ) are defined using a pair of partial differential equations per layer: where u LM , u I can be computed using the dimensional transmembrane potentials (V LM , V I ).
For more information on f (u LM ), g(u I ), β LM , β I , γ LM , γ I , F LM (u LM , u I ), F I (u LM , u I ), D LM , D I , D LM, I , D I, LM , readers may refer to [70] and the values of the parameters are provided in I of Table 2. The excitability parameter I is a function of the distance from pylorus (z-axis) that has been approximated from experiments with the following expression. Table 2. Parameters for the porcine intestine [70].  0 is the vacuum permittivity or the absolute dielectric permittivity of the vacuum. Its value is not specifically provided in [70]; however, it may be taken as 8.854 × 10 −12 Farads per meter.

Intestinal Active Viscoelasticity
Furthermore, Gizzi et al. [70] extended the SMCs in anisotropic tissue with two preferential directions a 1 and a 2 to describe the viscoelasticity of the intestine, and SEF takes the following form: ψ = ψ elas (F elas , a 1 , a 2 ) + ψ act (F, E, a 1 , a 2 ) + ψ visc (F), (12) where the elastic (ψ elas ), viscous (ψ visc ) and active (ψ act ) SEF are expressed as follows: where κ 0 , µ 1 , µ 2 , k 1,a 1 , k 2,a 1 , k 1,a 2 and k 2,a 2 are the material parameters, J elas ,Ī 1,elas ,Ī 2,elas , I 4,a 1 ,elas , andĪ 4,a 2 ,elas are the Jacobian of F, which are the first and second isotropic principal invariants and anisotropic principal invariants of fiber directions, a 1 and a 2 , respectively. The material parameters fitting the viscous active electromechanical model for porcine intestine obtained from the biaxial tensile experiments in [70] are shown in II of Table 2.

Discussion
Computational methods have proven to be a cost-effective technique to approximate almost every real physical problem for which appropriate mathematical formulations and computer models are necessary. Various phenomena associated with the large intestine have been studied computationally so far, such as smooth muscle contraction [64,69,70], active viscoelasticity [70], multi-layered passive hyperelasticity [10][11][12][13][14][15][16][17][18][19][20][21][22][23][24], post-surgical stapled anastomosis simulation [25] and colon peristalsis [71]. Although numerous studies had adapted simple isotropic constitutive models, attention was obtained only after the formulation of the HGO model, which was originally focused on the arterial material definition. Very few researchers [15,22,24] have studied intestinal mechanics biaxially. Howes and Hardy [15] performed equibiaxial tensile tests on human colon specimens using cruciform specimen and found higher average longitudinal second Piola-Kirchhoff failure stress (3.2 ± 1.51 MPa) than in the circumferential directions (2.35 ± 1.37 MPa). Furthermore, the Green-Lagrange failure strains in the longitudinal direction (0.139 ± 0.039) were lower than in the circumferential direction (0.158 ± 0.036). Compared to the uniaxial tension on human specimens, a late failure (λ ≈ 3) was achieved by Egorov et al. [11]. Likewise, the use of the cruciform tends to early failure, for which square specimens pulled by hooks are a better choice [72]. Siri et al. [22] were the first to succeed in separating the mice intestine layers and to find the submucosal layer as the prominent load bearing component, which is with stiffer in the longitudinal direction than in the circumferential direction. However, the study is limited to a very small strain (12%) and may not completely describe the tissue mechanics under higher loading conditions. Proximal to distal porcine colon specimens have been biaxially stretched by Puértolas et al. [24]. A great effort has been placed to characterize the studied tissues using the classical HGO model, the four-fiber model (Equation (2)), the microfiber von Mises model and the microfiber Bingham model. However, the study is also limited to the lower strain regimes (15%) and a further limitation of using frozen specimens may compromise true tissue mechanics [73].
In this experimental study, equibiaxial tensile experiments are performed on the porcine large intestine that showed highly nonlinear, anisotropic, hyperelastic tissue behavior. One of the important aspects of this study is the use of fresh tissue specimens. In contrast to the existing studies, efforts have been placed to achieve the complete failure of the specimen during experiments in this study. However, only data set up relative to the maximum stress are considered for investigating virgin passive mechanics. Therefore, the obtained experimental results are greatly consistent with the existing literature in important aspects. The intestine is stiffer in the longitudinal direction than in the circumferential direction [15,22] with larger average maximum stress in the longitudinal direction (Figure 4f). Four fiber families have been considered to best predict the intestine mechanics [24]; hence, they have been used to characterize and estimate the corresponding material parameters. Although not significant, the average submucosal layer stiffness (0.0136 MPa) is also found to be greater than the average muscular layer stiffness (0.0133 MPa) and the finding is qualitatively consistent with Siri et al. [22].
Another important aspect of this paper is numerical stability and the computational validation of the estimated virgin passive material parameters. Each parameter set passes through the positive eigenvalues of the positive definite Hessian matrix, satisfying the convex SEF at every point in the strain space. In addition, FE simulations are performed to validate the parameter, where equibiaxial tension of a linear brick element validates its numerical stability and provides a computational tool for further investigations. Lastly, the biaxial tension of an experimented specimen has also been computationally validated, with the obtained displacement contours agreeing qualitatively and quantitatively with DIC measurements.
The second part of this paper reviews non-passive large intestine mechanics. During the physiology of the large intestine to temporarily store and transport fecal material for defecation, passive mechanics are not sufficient, for which the active contractile elements in the intestinal wall have been carefully reviewed. The electro-chemical-mechanical interactions between contractile actin and myosin filaments regulated by intracellular Ca ++ ions are undoubtedly essential for the numerical study of the intestine in vivo. There are few models of active contraction in colon tissue, which are reviewed in [64,69,70]. However, there are several other approaches to model SMC contraction, which have been mainly used for arteries [74,75]. Stålhand et al. [76] formulated a chemomechanical finite strain model for SMC contraction based on a spring element and a contractile element in parallel controlled by Ca ++ -dependent state variables, as proposed by Hai and Murphy. The one and two-dimensional structural constitutive models by Tan and De Vita [77] and the three-dimensional phenomenological constitutive model by Schmitz and Böl [78] show a good approximation of the passive, active and total arterial mechanics. Chen and Kassab [79] described the three-dimensional microstructural response of coronary artery differentiating into elastin, collagen and SMCs. Dependent on myosin phosphorylation and cross-bridge kinetics, a dynamic muscle model has been proposed by Gestrelius and Borgström to estimate active force generation in SMCs [80]. Recently, Brandstaeter et al. [81] modeled the pacemaker's electrophysiology and SMCs' contractility in gastric tissue using an electro-mechanical constitutive model. By a multiplicative decomposition of the deformation gradient F = F elas F act with additional tissue incompressibility, J act = 1, the active part F act is computed for SMC. Similarly to intestines, the gastric tissue also consists of circumferential and longitudinal SMCs; therefore, the active gastric mechanics presented in Brandstaeter et al. possess relevance in modeling other GI tracts in the future.
This study is limited to experiments on the first four coils of the porcine large intestine. Studies have showed that the mechanical behavior of the large intestine varies throughout its length from the proximal colon to the rectum. In this regard, the investigation should be extended in the future to distinguish the anatomical sections and their respective mechanics. Furthermore, only five specimens have been characterized in this study that technically lacks statistical significance. However, the finding of this study is still relevant and the objective to present an appropriate numerical tool for the computational analysis of the multi-layered tubular architecture has been met.

Conclusions
The highly nonlinear, anisotropic hyperleastic material response of the large intestine has been found, which is consistent with the existing findings. The average failure stretch for the longitudinal direction is lower than the circumferential direction, while the maximum longitudinal stress is much higher than in the circumferential direction. Premature tissue tearing at needle locations has been a crucial problem in our biaxial system, which needs to be improved in the future. The authors would like to conclude that the adopted constitutive framework in this paper is able to characterize the multi-layered intestine tissue. Furthermore, the presented numerical tool is appropriate for predicting the underlying virgin passive colon mechanical response and can be extended in the future to investigate other underlying colon mechanics, such as active contraction, viscoelasticity, growth and mechanical damage.