Layer-Specific Damage Modeling of Porcine Large Intestine under Biaxial Tension

The mechanical behavior of the large intestine beyond the ultimate stress has never been investigated. Stretching beyond the ultimate stress may drastically impair the tissue microstructure, which consequently weakens its healthy state functions of absorption, temporary storage, and transportation for defecation. Due to closely similar microstructure and function with humans, biaxial tensile experiments on the porcine large intestine have been performed in this study. In this paper, we report hyperelastic characterization of the large intestine based on experiments in 102 specimens. We also report the theoretical analysis of the experimental results, including an exponential damage evolution function. The fracture energies and the threshold stresses are set as damage material parameters for the longitudinal muscular, the circumferential muscular and the submucosal collagenous layers. A biaxial tensile simulation of a linear brick element has been performed to validate the applicability of the estimated material parameters. The model successfully simulates the biomechanical response of the large intestine under physiological and non-physiological loads.


Introduction
The primary functions of the large intestine are to digest resistant starch, fiber and undigested proteins, to produce and absorb vitamins, electrolytes, short-chain fatty acids and amino acids, to absorb water, and to transport fecal matters outside the body [1,2]. Extending from the cecum to the anus, the intestinal wall is composed of four layers: outer serosa, muscle, submucosa and inner mucosa [3]. Each layer is a heterogeneous mixture of extracellular matrix and collagen fibers with two additional unidirectional smooth muscle fiber families in the longitudinal and the circumferential muscular layers [4,5]. Such arrangement offers a complex mechanical behavior characterized by a highly nonlinear anisotropic stress-strain relation, hysteresis and rate sensitivity, which have been widely demonstrated in uniaxial tensile, biaxial tensile and inflation experiments [3]. There is a large variability in the experiments and the constitutive models for the characterization of the intestines' mechanical behavior [3], which are logically relevant for approximating its stress-strain relations. However, all the approaches are limited to the lower elastic range or only to the ultimate stress level. Furthermore, Egorov et al. in 2002 [6] observed two maxima in the human colorectal tissue tested in longitudinal direction, postulating the rupture of different layers at different stretch levels. The first higher narrow maximum at around 40% strain is hypothesized due to muscular rupture, and the second lower broad maximum after 80% strain due to the submucosal rupture. For almost two decades, such unique tissue mechanics have been completely neglected in the available characterization approaches, if they have ever been encountered at all. Most importantly, a correct selection of the constitutive model that is correlated with the experimental methodology, and capable of fully characterizing the material properties of such anisotropic soft tissues are necessary [7].
Intestines are subjected to progressive damage by remodeling during diseased condition or by staples/sutures after surgical anastomosis [8]. For example, ulcerative colitis is characterized by depleted mucosal lining, disrupted-thinned muscular collagen fibers [9], thickened submucosa and increased cross-sectional area [10,11]. Such microscopic changes result in different tissue mechanics than the healthier counterparts [10,11], which affects the normal physiological function of the intestine. Following the concept of effective stress reduction factors initially proposed by Kachanov [12], several deterministic, probabilistic and microstructural models have been proposed to describe the damage induced mechanical behavior of fiber-reinforced soft biological tissues [13]. Some of the recent structure-based damage evolution formulations for ligament, rectal sheath and mostly artery have been proposed by Balzani et al. [14], Calvo et al. [15], Volokh et al. [16], Peña et al. [17], Ehret & Itskov [18], Martins et al. [19], Comellas et al. [20], and Noble et al. [21]. The common assumption in all the formulations is that the global damage affects only the isochoric/volume-preserving or deviatoric/distortional part of the deformation as proposed by Simo [22], which can be additively decomposed into contributions from the isotropic ground matrix material and the anisotropic collagen fibers.
Research on intestine biomechanics is aimed to better understand and model anisotropy, viscoelasticity and heterogeneity [4,6,9,[23][24][25][26][27][28][29][30][31][32][33][34][35]. To the best knowledge of the authors, this paper is the first to investigate mechanically-induced damage to the large intestine. Although a phase-field approach to model the fracture in the small intestine wall has been recently presented by Nagaraja et al. [36], the study is limited to uniaxial tension experiments, nonhysteresis mechanics, and overall damage parameters rather than a layer-or-fiber specific fracture analysis. The continuum damage evolution proposed by Comellas et al. [20] has been adapted in this paper to constitutively describe the true stress-strain response of the porcine large intestine during biaxial tensile loading until complete failure. Knowing the fact that mucosal and serosal collagen fibers do not contribute to the in-plane stiffness significantly [6,37], they have been therefore ignored in this study. Therefore, the overall tissue stiffness is provided by the extracellular matrix, the submucosal collagen, the longitudinal and the circumferential muscular collagen fiber families. For this, the fresh porcine intestine specimens have been biaxially stretched. The simple neo-Hookean and the well-known Holzapfel-Gasser-Ogden (HGO) constitutive model [38] have been used to describe the isotropic ground matrix and different collagen fiber families, respectively. The eventual damages on the load paths are defined by separate contributions of damage from each load-bearing collagen fiber family.

Equibiaxial Tensile Specimen Preparation
The large intestines of four pigs were obtained from the Institute of Laboratory Animal Science at the University Hospital RWTH Aachen, Germany. The experimental pig was primarily used for the medical training purposes and euthanized and dissected. Thus, the obtained intestine as a whole organ (Figure 1a) is covered with cloth soaked in a 0.9% NaCl solution to prevent the tissue from desiccation. Within an hour, the lump of a fresh intestine is brought to the testing facility (Biomechanics Laboratory, FH Aachen University of Applied Sciences, Jülich, Germany), divided into segments based on five anatomical coils, and completely cleaned of feces and other connective/adipose tissue. Each coil is taken in turn from which individual specimens with an edge length (L) of 40 mm × 40 mm are prepared and the remaining tissue is stored in the 0.9% NaCl solution until the next testing. Close inspections have been done during preparation and the specimens with any preliminary damage are excluded for the biaxial tensile experiment. Specimen thickness was measured with a micrometer screw prior to testing.

Test Protocol
The equibiaxial tensile tests of the intestine specimens are carried out on the selfdeveloped biaxial machine "BiAiX" (FH Aachen, Aachen, Germany, patent pending, DE 10 2017 116 067 A1) (Figure 1b). The BiAiX machine has already been used for several experiments, and the results have been published in a peer-reviewed journal [35]. It comprises of four arms with force measuring actuators and five light-weight aluminium cylindrical bars with a sharp needle at their free end for tissue fixation. The needles can move laterally to allow free stretching of each side of the specimen. For 2D planar evaluation, one Allied Vision Prosilica GT2450, 5 Megapixel camera (Allied Vision Technologies GmbH, Statroda, Germany) with a GigE Vision Gigabit Ethernet interface and a maximum frame rate of 15 f/s with AQUAMARINE 2.0/28 C lens (Schneider Kreuznach, Bad Kreuznach, Germany) lens is mounted above the specimen position in 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 digital image correlation (DIC) with a calibration target. The specimens are then mounted in the apparatus ensuring the orientation of the sample always the same: longitudinal direction along the L-L axis and circumferential direction along the C-C axis (Figure 1b). To ensure a uniform evaluation of the specimens, the inner mucosal layer of the specimens is always oriented facing up and a graphite pattern (black random spots in Figure 1c,d) is sprayed for DIC. Polarized light with a polarizing filter in front of the camera lens is used to avoid reflections on the wet specimens. The measurement protocols for the tensile test are as follows: five cycles of preconditioning followed by a monotonous stretching both at the speed of 10 mm/min until failure and a frame rate of 1 Hz. Figure 1, panels c and d, show specimen images in the undeformed and highly stretched configurations, respectively with fixed camera-lens position.

Multilayer Passive Anisotropic Strain Energy Function
Recent advancement in the histological techniques such as nonlinear second harmonic generation (SHG) has greatly improved the architectural description of the multilayered intestine [5,37,39]. Maier et al. [5] have found the orientation and thickness of the collagen fibers in the submucosa layers to be ∅ 6.5 µm, aligned in a crosswise arrangement at α = ±30 • to the longitudinal direction ( Figure 2). Moreover, collagen fibers (∅ 1.9 µm) are also present in the longitudinal and the circumferential muscular layers. The collagen fiber orientations in different layers can be defined by the angles measured between the fiber families and the Cartesian coordinate axis (e 1 ) as: where the collagen fibers in the longitudinal and the circumferential muscular layers are perpendicular to each other ( Figure 2). Within the framework of nonlinear elasticity, such fiber-reinforced composites can be mechanically described by highly deformable hyperelastic constitutive model expressed in terms of the quasi-isochoric (negligible volume change under loading) strain energy functions (SEF). The decoupled representation of SEF adopted from Holzapfel et al. [38] ψ(J, C) := ψ • vol (J) + ψ isch (Ī 1 , a 01 , a 02 , a 03 , a 04 ), is based on the kinematic multiplicative decomposition of the deformation gradient tensor, F = F • F into a spherical or volumetric (F • = J 1 3 I) and a unimodular or isochoric part (F = J − 1 3 F), [40]. J := det[F] is the determinant of the deformation gradient tensor, a measure of the volume change. It allows the definition of the isochoric right and left Cauchy-Green deformation tensors as: is the purely volumetric contribution of the SEF and κ 0 is the initial bulk modulus. The purely isochoric term ψ isch can be additively split into the isotropic part (ψ iso ) from the matrix and the anisotropic part (ψ aniso ) from the collagen fibers. The isochoric SEF in Equation (2) can be expressed as: ψ isch (Ī 1 , a 01 , a 02 , a 03 , a 04 ) = ψ iso (Ī 1 ) + ψ aniso (a 01 , a 02 , a 03 , a 04 ) = ψ iso (Ī 1 ) + ψ aniso (Ī 4,a 01 ,Ī 4,a 02 ,Ī 4,a 03 ,Ī 4,a 04 ), whereĪ 1 := tr[C] is the modified first principal invariant of the symmetric modified right Cauchy-Green tensor andĪ 4,a 0i := C : A 0i = C : (a 0i ⊗ a 0i ) = a 0i · C · a 0i = λ 2 i represent the anisotropic invariant, which is quantitatively equal to the square of the fiber stretches λ 2 a 0i associated with the i th collagen fiber family, a 0i . The deformation gradient F transforms the reference fiber direction into current configuration, a i = F[a 0i ] and the respective structural tensors, At low strain, the tissue mechanical behavior is solely due to the ground matrix and can be defined using the classical neo-Hookean material model: where µ 0 > 0 is the initial shear modulus and has the unit of stress. Siri et al. [37] measured similar mechanical behavior of the submucosa and the muscular layers. Therefore, each layer in this study has been modeled with the same form of the HGO constitutive model, represented by an exponential function [38]: where k 1,a 0i > 0 are stress-like material parameters and k 2,a 0i are dimensionless parameters for the i th collagen fiber family and influence the intestine mechanics at high strain only resulting in a polyconvex SEF. For more detail, refer to Bhattarai et al. [35]. Considering the measurement of the muscular (mus) layer collagen fibers (∅ 1.9 µm) by Feng et al. [39], the longitudinal and the circumferential muscle layer both can be assumed to be transversely isotropic, governed by the common material parameters, k 1,mus and k 2,mus . Likewise, the submucosal (sm) collagen fibers (∅ 6.5 µm) arranged in a crosswise alignment and assuming no fiber interaction contribute to the second part of the anisotropic stiffness for which k 1,sm and k 2,sm are the material parameters. Finally, the anisotropic HGO constitutive model for the intestine (Equation (6)) can be re-written as: where ψ mus (Ī 4,a 01 ,Ī 4,a 02 ) = satisfying the conditionĪ 4a 0i andĪ 4a 0j ≥ 1. Adding Equations (3), (5) and (7), the original intestinal SEF Equation (2) becomes:

Anisotropic Hyperelastic Stress Calculation
The mechanical description of the intestine wall requires the derivation of an appropriate stress as a function of strain. The hyperelastic response of the intestinal tissue defined by the SEF in Equation (8) can be written in terms of the second Piola-Kirchhoff stress (S) with respect to invariants [22]. Using the chain rule, the Cauchy stress (σ) can be computed as: where p = κ 0 (J − 1) is the hydrostatic pressure; I is the second order identity tensor, is the material deviator of the structural tensor A 0i with the right Cauchy-Green tensor C operating as metric tensor and dev ]. The derivative of the j th anisotropic SEF component related to the muscular and submucosal collagen fibers can be computed as:

Nonlinear Damage Constitutive Model
In this section, the basic formulation of the damage evolution developed by Comellas et al. [20] has been adopted and presented for the investigation of the damage in the intestine material. Damage of any material is characterized by the reduction of the stiffness at stretch beyond the physiological range, due to microstructural failures, followed by an eventual fracture of the specimen. Initially postulated by Simo [22], damage phenomena are assumed to affect only the isochoric parts of the constitutive Equation (8). However, the damage in each of the components has to be distinguished in order to provide a complete description of the intestine damage. Weisbecker et al. [41] performed layer-specific damage experiments, and used a pseudo-elastic damage model to find aorta damage primarily induced by the collagen fibers. Furthermore, the non-collagenous matrix can be completely damaged well below the physiological loading range. Hence, only the collagenous layer damage is considered in this paper as: where ψ iso , ψ mus , ψ sm are the undamaged isochoric SEF for isotropic matrix, muscular and submucosal layers, respectively. The Kachanov-like reduction factors (1 − D lm ), (1 − D cm ) and (1 − D sm ) are the functions of the normalized scalars with D k ∈ [0, 1] defined as internal damage variables, where k = lm for longitudinal muscular layer collagen fibers, cm for circumferential muscular layer collagen fibers, and sm for submucosal collagen fibers.
Utilizing the concept of Clausius-Duhem inequality for isothermal cases, the internal dissipation can be generalized as where the rate of the SEF in Equation (11) using the chain rule can be written aṡ Comparing the expressions in Equations (12) and (13), the Kachanov effective second Piola-Kirchhoff stress for finite strain and the non-negative internal dissipation can be deduced as respectively. Inequality (14) shows that damage is a dissipative process which implies nondecreasing damage variables, i.e., D k ∈ [0, 1] &Ḋ k ≥ 0 ∀k = {lm, cm, sm}. Applying the push forward, the Kachanov effective Cauchy stress tensor can be obtained as

Damage Evolution
At any current loading time t, Simo [22] proposed φ t,k = G(τ t,k ) − G(τ max t,k ) ≤ 0 as the damage criterion in the strain space with the condition φ t,k = 0 defining the damage surface.
Here, τ t,k = 2ψ k (C(t)) is the equivalent strain as the norm at current time t related to the undamaged deviatoric SEF, ψ k (C(t)) and τ max t,k = max s=0,t 2ψ k (C(s)) is the maximum value up to current time t or the current damage threshold which satisfies τ max t,k ≥ τ d 0,k , with τ d 0,k as the initial damage threshold stress. Damage is initiated, if the current τ t,k ≥ τ max t,k which is equal to the initial damage threshold stress at the beginning of the deformation, i.e., at t = 0, τ max 0,k = τ d 0,k . To complete the constitutive model with failure, the exponential evolution formulation used by Comellas et al. [20] is adapted in this paper: the evolution of the damage variables is given as where τ d 0,k is the initial damage threshold stress; g d f ,k is the fracture energy per unit volume, and are the material parameters to be determined through fitting the experiment data. The damage variables for each component are defined for the interval D t,k ∈ [0, 1] that satisfies D 0,k = (τ 0,k ) = 0 and D ∞,k = G(τ ∞,k ) = 1.

Numerical Simulations
Numerical simulation has been a cost-effective technique to validate the feasibility of the adopted mathematical formulations and the reliability of the estimated parameters that describe the mechanical behavior of any experimented material. In this study, the finite element (FE) analysis is performed in the commercial FE software LS-DYNA (Ansys, Inc., Canonsburg, PA, USA) Version: smp d R11.1.0 in 64 bit Windows 10. A user-defined Fortran material subroutine that takes into account a) the hyperelastic response of the four collagen fibers, and b) the damage evolution in the load bearing collagen fibers (Equation (15)) has been coded inside LS-DYNA. The first part of the subroutine implementation has been successfully tested in Bhattarai et al. [35]. An example is presented in this paper to demonstrate the numerical performance of the damaged induced anisotropic hyperelastic model and the estimated material parameters. A cube of edge length 1 mm meshed with one linear brick element and reinforced by four fiber families (Figure 3) is considered. The four collagen fiber families are represented by unit directional vectors, a 01 , a 02 , a 03 , a 04 . For biaxial tension in the XY plane, four nodes in the bottom surface are constrained to move along the z-axis, while equibiaxial strain is controlled by displacements, which are prescribed on all nodes in the XY plane.

Equibiaxial Tensile Experiments
The mean specimen thickness varied from 0.68 to 2.15 mm, with a mean of 1.188 mm (N = 102). Displacement-controlled (10 mm/min) equibiaxial extension has been applied on the intestine specimens, which took 12 min on average to completely rupture each specimen. The tensile force (P L and P C ) are recorded at each load step on actuators of the longitudinal and the circumferential directions spanning over 3.3-17.0 N and 3.8-19.55 N, respectively.
The experimental methodology reported in this paper endures a technical problem, which limits the number of valid tests. As shown in Figure 4, the tissue tearing at the fixation with increasing stretch is merely inevitable due to early global fracture at the needle insertion points causing a sudden contraction of the specimen in the opposite of the loading direction. This phenomenon is no longer continuum damage and does not provide information about the intended damage in the tissue. In contrast to the expected stress-softening, unloading like behavior, i.e., backward and downward path of force-displacement or stress-strain curves occurs, which is purely due to tissue tearing at the holes making such data invalid for the continuum damage analysis. All the tested specimens in this study that showed global fracture at the needle insertions were therefore excluded from the damage analysis.
For the DIC evaluation with the ISTRA 4D (Limess Messtechnik und Software GmbH, Krefeld, Germany) software, residuum = 40-50, facet size > 50 and grid spacing = 1/3 facet size > 17 are used to obtain the homogeneous Green-Lagrange strain field without a significant loss of the grid points. The Green-Lagrange strain (E LL , E CC ) in both loading directions are evaluated directly from the heat maps using the DIC technique in the ISTRA 4D software for each individual image and recorded as an averaged value from a chosen deformation region (Figure 4). It is interesting that inhomogenous and anisotropic strain distributions are observed. In specimens S1, S3, and S5, several in-plane ripples are observed that look similar to calm ocean waves directed vertically for the longitudinal load direction and horizontally for the circumferential direction. This is not clearly understood, but the underlying collagen fibers might have withstood the applied stress. A higher magnification and better resolution camera, if available, may be used to resolve such strain distribution.  . Green-Lagrange strain heat maps (E LL , E CC ) measured along the longitudinal and the circumferential direction in the specimens. "SiL" represents the i th specimen along the longitudinal direction, whereas "SiC" represents the i th specimen along the circumferential direction.

Damage Based Anisotropic Intestine Mechanics
Using the incompressibility condition, the true or Cauchy stress (σ LL = P L λ L /LT and σ CC = P C λ C /LT) are set as a function of the elongation, λ L = √ 2E LL + 1 and λ C = √ 2E CC + 1. Figure 5 shows the representative Cauchy stress-stretch curves of the intestines from four pigs in the longitudinal (blue curves) and in the circumferential (red curves) directions under equibiaxial tensile stretching. The porcine intestine showed clear anisotropy with stiffer mechanical response in the longitudinal direction than in the circumferential direction. The standard J-shaped stress-stretch curves are observed for all specimens followed by the reduction of the stress slope induced by microstructural failure of the intestine wall composure. Furthermore, different tissue failure stretches are observed in both directions. Therefore, different values of the damage variables are expected in each direction due to progressive damages of individual fibers during the course of the tensile stretching. Significantly lower (by approximately two orders of magnitude) shear strains than the normal strain are found and are therefore neglected in this study. Cauchy stress (MPa) S1L Exp S1C Exp S1L Fit S1C Fit S1L Simul S1C Simul [S1] [S4] Figure 5. Fitting of the representative stress-stretch curves using damage formulation (Equation (17)). "Exp" represents the experiment data set, "Fit" represents the Matlab curve fit data set, and "Simul" represents the simulation values.

Damage Parameter Estimation
To determine the damage driven passive, incompressible, anisotropic, hyperelastic material parameters as per the SEF in Equation (15), the curve fittings are performed in Matlab (The MathWorks, Inc., Natick, MA, USA). A careful selection of the experimental curves is made: the non J-shaped and the specimens without concurrent softening stress stretch curves in both longitudinal and circumferential directions are excluded for the model parametrization. The coefficients of determination R 2 ∈ [0, 1] in the longitudinal (R 2 L ) and in the circumferential (R 2 C ) directions, respectively, are used to determine the quality of the fit for the estimated material parameters, R 2 is calculated as follows: 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. A value of R 2 closer to unity indicates a good fit to the experimental data. The fitted stress-stretch curves are shown in Figure 5 and the material parameters are listed in Table 1. For each specimen, unique discontinuities or sudden changes in the stress-stretch slopes are associated with the damages in the longitudinal, the circumferential, and the submucosal collagen fibers governed by the initial threshold stresses and the fracture energies for each layer (Table 1). Figure 6 shows the predicted layer-specific damages for all the fitted specimens. In 50% (S1 and S3) of the studied specimens, the circumferential collagen may initiate damage followed by the submucosal collagen and the longitudinal collagen. However, the longitudinal layer shows a brittle-like failure at a lower stretch than the circumferential and the submucosal layers.

Numerical Simulation of the Biaxial Tension
Based on the user defined material subroutine, the biaxial tension of all pigs has been simulated using a linear brick element as shown in Figure 7. Using LS-DYNA with 16 CPUs, the total simulation time is only 2 min. Under prescribed elongations, DX = 0.31 mm and DY = 0.38 mm and using material parameters as given in Table 1, the brick model is contracted by DZ = −0.6493 mm along the z-axis (Figure 7). The Cauchy stress and the strain data in the linear element are extracted from LS-PrePost. The respective stretches are computed, and the Cauchy stress vs. stretch curves are plotted over the experimentfit dataset ( Figure 5). The details of the experimented large intestine mechanics, such as anisotropy, hyperelasticity and damage, have been appropriately generated, i.e., the material behavior for each specimen is qualitatively and quantitatively similar to the dataset. This confirms that the estimated material parameters are numerically stable, and the implemented Fortran user material subroutine is capable for further numerical investigations.

Large Intestine Mechanical Behavior
Abundant resources on the large intestine biomechanics can be found to describe its mechanical behavior that are used for various computational applications [3,8,25,27,28,37,[42][43][44][45]. However, these studies are limited to the uniaxial tensile tests, non-fresh experimental specimen, and undamaged passive anisotropic hyperelastic material modeling. For a complete and realistic material description, appropriate experimental techniques and computational models are necessary. Focused on such aspects, equibiaxial tensile experiments of fresh porcine colon, and the layer-specific orthotropic material description beyond the supraphysiological loading have been performed in this study. In general, the large intestine is found to be clearly orthotropic, with the standard J-shaped nonlinear stress-stretch paths under tensile loading ( Figure 5). Furthermore, the stiffness in the longitudinal direction is found to be larger than the stiffness in the circumferential direction, with greater average failure stretch in the latter direction (1.37 vs. 1.52). Thus, the results obtained are consistent with the previous studies. For example, dynamic equibiaxial tension tests were performed on the human cruciate colon specimens by Howes and Hardy [27] and observed a greater average failure 2 nd Piola-Kirchhoff stress (3.20 ± 1.51 MPa) and a lower average Green-Lagrange failure strain (0.139 ± 0.039) in the longitudinal direction, than in the circumferential direction (2.35 ± 1.37 MPa) and (0.158 ± 0.036). Likewise, Sokolis and co-authors performed inflation/extension tests on rat proximal colon and rectum specimens [4,45]. Each section of the intestine possessed greater longitudinal stiffness than the circumferential stiffness for which the Fung exponential [46] and the HGO constitutive models [38] were used. Similar differences between two directions have been obtained in biaxial tension by Siri et al. [32,37] and Puértolas et al. [34]. However, these studies are limited to the lower strain regime (≈15%) and may not provide complete information of the tissue mechanics.
Other aspects in the constitutive modeling of the fiber-reinforced anisotropic soft tissues that are not discussed in this study are in-plane fiber dispersion, interaction between two intertwined fibers, and the fiber shear contribution [47][48][49][50]. For any two intertwined fiber families, e.g., submucosal collagen fibers, a 03 , a 04 , the anisotropic invariants I 5 = a 03 · C 2 · a 03 and I 7 = a 04 · C 2 · a 04 , contributes to the fiber shear deformation, which have a strong correlation with I 4 and I 6 . However, they may lead to an ill-posed parameter estimation problem and are usually not included in the constitutive formulation [51], whereas the invariant, I 8 = (a 03 · a 04 )a 03 · Ca 04 represents the interaction between the in-plane intertwined fiber families [52]. Likewise, an additional tissue stiffness can be provided taking into account the fiber dispersion within the tissue by modifying the fiber family function, Equation (7) based on the Generalized Structure Tensor approach as: where κ (not to confuse with κ 0 = initial bulk modulus) is a measure of dispersion in the fiber orientation. For κ = 0, the model reduces to perfectly aligned fibers and the case associated with κ = 1/3 is isotropically distributed fibers. For more information on the derivation and its application, see [47,48,53].

Damage Evolution of the Intestine Layers
Existing studies on the large intestine are limited to the undamaged passive anisotropy and do not include the damage evolution in the mathematical modeling. To the best of the authors' knowledge, Nagaraja et al. [36] are the only work so far that models the fracture in the intestinal wall, but uses it for the small intestine. With a phase-field modeling approach, the failure of the rectangular notched strip is governed by the fracture toughness (G C > 0) with two additional anisotropy material parameters. Using the neo-Hookean and the four-fiber family based HGO strain energy function, the standard intestine anisotropy (stiffer longitudinal direction) has been well predicted. However, the uniaxial tensile testing and the non-layer-specific fracture parametrization do not fully differentiate the damage progressions in all load-bearing layers.
In order to fill this gap, we emphasize the characterization of the equibiaxially stretched large intestine based on the damage evolutions in all the load-bearing collagenous layers. Herein, the collagen fiber families in the longitudinal muscular layer, in the circumferential muscular layer, and in the submucosal layer contribute to the applied load. Therefore, damage (D k , ∀k = lm, cm, sm) is calculated separately for these fiber families depending on their stiffness related material parameters (µ 0 , k 1,i , k 2,i , ∀i = mus, sm), initial damage threshold stress (τ d 0,k ), and fracture energy (g d f ,k ). The influence of the damage related parameters is such that, the lower the value of threshold stress (τ d 0,k ) for a particular layer, the earlier the damage initiates in that layer. Likewise, the lower the value of fracture energy (g d f ,k ), the steeper the negative slope of the damage stress-stretch curve becomes. From the curve fit using Equation (16) in Figure 6, the following sequences of the damages in the intestine layers have been predicted. The specimens S1 and S3 may show early circumferential layer damage, followed by the submucosal and the longitudinal muscular layer. In specimen S2 with the lowest longitudinal stiffness among all, the longitudinal layer damage is immediately followed by the submucosal layer and the circumferential layer. However, in specimen S4, the submucosal layer damage may have followed by the longitudinal and the circumferential. In all the specimens, the instantaneous (brittle-like) failure of the longitudinal muscular layer is typically observed right after the damage initiation. Furthermore, the circumferential layer and the submucosal layer failure show slower failure that progresses with increasing stretch. Although the experimental data-set fit well in all the presented specimens, damage sequence and evolution are not consistent. The biological variance in animal could be the reason why the specimens from each subject are unique. Therefore, a larger number of specimens is necessary from more animal subjects and is essential to conclude the exact event of the layer failure.
In modeling damage of the fiber-reinforced anisotropic material, localization occurs due to spurious mesh dependency of the softening constitutive model [54]. It means that, when failure initiates in a particular finite element, depending on the material properties and boundary conditions, the deformation increases rapidly and the element collapses. At the same deformation state, the neighboring elements may still be intact to resist more deformation. This generates rapid stress oscillations between such adjoining elements, and a homogeneously distributed stress cannot be achieved. Such pathological mesh dependent effects can be limited by the nonlocal strategy: introduction of an additional non-local variable in the constitutive behavior at a material point, which is dependent on the neighboring material point. Therefore, instead of using the local damage parameter, an integral average damage can be calculated within the interaction domain or the radius of influence, commonly called an internal length that regulates the size of the damaging zone. The non-local formulation based on Lemaitre's model may greatly benefit simulating damage based impact and crash simulations. For more information, see [55][56][57].

Limitations of This Study
This study possesses numerous advantages such as appropriate biaxial tensile experiments, realistic supraphysiological loadings for damage evaluation, fiber-reinforced anisotropic mathematical modeling, and numerically stable material parameters validated by the finite element simulations. Each specimen has shown a strict repeatability and consistency in relation to maximum load, and higher stiffness in the longitudinal direction than in the circumferential direction. However, this study is still limited to only four representative experiments on four pig subjects, which is not statistically significant. Studies have shown that the mechanical behavior of the large intestine varies throughout its entire length from the proximal colon to the rectum. In this regard, the investigation direction should be further extended in the future to an adequately large number of pig subjects throughout its entire length from the proximal colon to the rectum for better statistical comparison and biological variance. Furthermore, one of the major difficulties in the biaxial tensile test using needles is the localized premature tearing, for which a better fixation technique would be helpful. Finally, in order to predict the layer-specific damage evolution, experiments on the separated intestinal wall would be more realistic. However, with our current experience, the separation has not been straightforward because the tissue was traumatized in the process before the experiment, rendering the specimen useless for the mechanical damage analysis. Therefore, a better technique is desired if the objective is to be achieved. Nevertheless, the numerical approach presented in this study is still relevant and the objective to present an appropriate numerical tool for the computational analysis of the damage in the multi-layered tubular architecture has been met.

Conclusions
Fiber-reinforced gastrointestinal biomechanics is not a very new subject of study, but damage has not yet been explored. Damage mechanics is important, especially for the lower bowel diseases and associated treatments. In such a situation, the load bearing collagen families across the thickness are greatly remodeled, altering the stiffness in both longitudinal and circumferential directions. Consequently, the load-stretch resisted prior symptoms will no longer be safe or damage-free. Likewise, for symptoms that cannot be easily treated with medication, anastomosis is commonly performed in which the diseased part is dissected, and two tissue sections are rejoined. On doing so, the tissue junction can be severely traumatized, both globally and locally. In order to achieve a successful treatment technique that is as safe as possible, an appropriate experimental methodology, and a realistic material description is necessary. With a focus on these aspects, the equibiaxial tensile experiment of the fresh porcine large intestine and the layer-specific anisotropic material description beyond the supraphysiological loading have been performed, and the results are presented in this study. Institutional Review Board Statement: Informed consent was obtained from all authorities to use the intestine tissues from dissected pig subjects involved in the study, especially the University Hospital RWTH Aachen, Germany for the usage and the disposal of the experimented tissue and the FH Aachen, Germany to perform the experiment in the University facility. According to the European Directive 2010/63/EU on the protection of animals used for scientific procedures and the German Animal Welfare law, ethics approvals are not required if the samples, cadavers or specimens are being supplied from an already approved research, teaching or breeding project. This procedure was approved and documented by the Animal Welfare Officer of the RWTH Aachen University and the responsible IACUC.