Computational model of gastric motility with active-strain electromechanics

We present an electro-mechanical constitutive framework for the modeling of gastric motility, including pacemaker electrophysiology and smooth muscle contractility. In this framework, we adopt a phenomenological description of the gastric tissue. Tissue electrophysiology is represented by a set of two minimal two-variable models and tissue electromechanics by an active-strain ﬁnite elasticity approach. We numerically investigate the implication of the spatial distribution of pacemaker cells on the entrainment and synchronization of the slow waves characterizing gastric motility in health and disease. On simple schematic model geometries, we demonstrate that the proposed computational framework is amenable to large scale in-silico analyses of the complex gastric motility including the underlying electro-mechanical coupling.


INTRODUCTION
Gastric motility is a coordinated dynamics of the stomach wall that typically results in rhythmic contractions, which accomplishes several physiological functions such as mixing, grinding and propulsion of ingested food. [1] Gastric motility is partly governed by the electric activity of the gastric electrophysiological system. [2,3] The electric activity of the gastric tissue controls the contractile stress exerted by the smooth muscle cells (SMC) of the tissue via synchronized propagating bioelectrical waves, known as slow waves. [4] Slow waves are permanently present in most of the regions of the gastrointestinal tract. However, smooth muscle contraction only occurs if certain neuronal and hormonal signals coincide with these electrical signals. [5] Slow waves are generated by specialized pacemaker cells, termed interstitial cells of Cajal (ICC), i.a., located between the longitudinal and circular muscle layer. Gastric slow waves are characterized by four distinct phases: in the initial upstroke phase, the cell depolarizes rapidly from a resting potential to a peak potential. The initial upstroke is followed by a pronounced plateau phase, in which the potential stays close to the peak potential. In the repolarization phase, the membrane potential returns to its resting value and during the refractory phase, it remains there. It has been shown that Ano1 channels, a type of Ca 2+ -activated Cl − channel, play a fundamental role in slow wave generation of ICC. [6] It is hypothesized that calcium-induced calcium release from internal stores into microdomains might be an important mechanism for pacemaker functionality. [7] ICC are spatially distributed in a continuous network-like arrangement. Within this network, slow waves propagate across ICC via voltage-dependent mechanisms. [3,8] Therefore, the pacemaker system of the stomach is spread across the whole organ, unlike in the heart where specific nodes dictate the pacemaker activity. [9] Overall, several details of the complex biophysical mechanisms ruling the electrophysiology of slow waves remain unknown, despite considerable research efforts. [2,3,6,7,[10][11][12][13] While isolated ICC generate slow waves at distinct intrinsic frequencies, they synchronize to the highest frequency in an undamaged network. This process is called entrainment. In the healthy human stomach, the driving pacemaker region is located on the greater curvature of the proximal stomach. There, the ICC with the highest intrinsic frequency of approximately 3 cpm (cycles per minute) can be found. From the pacemaker region, the intrinsic frequency degrades towards the pylorus, a fact known as (intrinsic) frequency gradient. In the distal antrum the instrinsic frequency is between 0.8 cpm to 1.8 cpm. The frequency gradient, combined with the process of entrainment, is essential for the coordinated onset and propagation of slow waves in distal direction resulting in the characteristic shape of ring wavefronts. Under healthy conditions, between three and four slow waves simultaneously travel at a distance of approximately 60 mm towards the pylorus. The fundus has been shown to be electrically quiescent so that slow waves do not spread proximally from the pacemaker region into the fundus. [9,14,15] In the presence of pathologies, abnormal slow waves are observed and associated with different gastric motility disorders. In particular, the physiological entrainment regulated by ICC may degenerate leading to ectopic initiations, conduction blocks etc., [16] usually named dysrhythmias. Most of these phenomena appear astonishingly similar to what is known from cardiac tissue, e.g., in the form of cardiac arrhythmias, [17] further emphasizing the shared biophysical features between different active biological media. [18] The specific electro-mechanical activity of the stomach is realized through the coordinated interplay between ICC and SMC. The slow waves generated and actively propagated by ICC spread passively to the surrounding SMC via dedicated proteic connections, the so-called gap junctions. [19] These localized structures furnish a direct electrical coupling between different cells forming a so-called syncytium. SMC are depolarized by the conduction of slow waves. The depolarization activates L-type voltage-dependent Ca 2+ channels (VDCC); [3] the initial event to start complex multiscale mechanisms, e.g., calcium-induced calcium release, [20,21] which result in the smooth muscle contractions deforming the gastric wall as a whole in a coordinated and synchronous manner.
In the last decades, simplified phenomenological, mathematical approaches, based on FitzHugh-Nagumo-like dynamical systems, have been proposed for modeling the ICC-SMC electrical coupling both in the intestine [22,23] and the stomach. [24] In these models, the coupled ICC-SMC system is represented with a minimal level of complexity by a pair of two state variables for each cell type, namely the normalized transmembrane voltage, , and a local recovery variable, . Unfortunately, this specific modeling approach suffers from an inherent instability which typically arises for simulation times longer than just a couple of seconds. The approach is thus practically not applicable for long-term in-silico analyses. Moreover, despite significant and ongoing research efforts (see Du et al. [16] for a recent review), several mechanisms controlling gastric motility are still not yet fully understood, which makes it difficult to develop detailed mechanistic models. Therefore, there is a pressing need for a simple phenomenological continuum-scale model of gastric electrophysiology that is computationally robust, efficient and stable.
Modeling gastric motility requires, however, not only a model of gastric electrophysiology but also of the active muscular contractions and electro-mechanical coupling giving rise to the peristaltic waves observed during the digestion of food. [16] Active biological tissues are intrinsically characterized by multiphysics and multiscale behavior involving several length and time scales. [25,26] The study of such complex microstructured materials fosters a continuous forefront research in applied mathematics and mechanics, including innovative theoretical formulations [27] and homogenization procedures, [28,29] as well as experimental [30][31][32][33] and numerical tools. [34,35] So far advanced theoretical and computational studies of soft tissue contractility have addressed mainly the nonlinear behavior of skeletal [36][37][38] and cardiac muscles [39][40][41][42] including a wide and variegated literature. [43] Much less attention, though, has been dedicated to modeling the multiphysics of gastrointestinal motility, which distinguishes itself by a particularly high structural and functional complexity. [3,44,45] Only a few contributions have so far addressed the mechanics of active contractions in the gastrointestinal tract. [46,47] For the small intestine, Gizzi et al. [48] recently proposed a generalized thermodynamical framework with a statistical description of the intestinal wall microstructure. [49][50][51] Moreover, some attempts to characterize the altered mechanical contractility of human colon strips in the presence of unbalanced thermal states were proposed over the last years. [23,52] What remains wanted so far is a computational model that combines a simple phenomenological, and computationally robust model of gastric electrophysiology with an electro-mechanical finite elasticity framework, enabling thereby comprehensive studies of the complex multiphysical nature of gastric motility on the organ level.
The present contribution aims at partially filling this gap. We propose a simple phenomenological approach coupling a modified-version of the two-variable Mitchell-Schaeffer electrophysiological model [53,54] with an active-strain electromechanical description of the gastric tissue. [40,55] Our aim is to establish a robust computational framework amenable to large scale in-silico analyses of the complex gastrointestinal motility including the underlying electro-mechanical coupling in a finitestrain framework. We directly relate the contractility dynamics in the stomach with the transmembrane voltage of SMC neglecting more complex biophysical phenomena due to neural, hormonal and mechano-electric effects. Our modeling approach can successfully reproduce several experimentally observed key phenomena as indicated by a number of computational tests. We demonstrate the broad applicability of the model by simulating both normal and dysrhythmic dynamics within a finite elasticity electro-mechanical framework varying both electrophysiological material properties and mechanical boundary conditions. This paper is mainly dedicated to general, phenomenological aspects of the mathematical modeling of the electro-mechanical coupling in the stomach. Therefore, all the computational examples were accomplished on simplified model domains such as rectangular tissue patches or cylindrical tubes, which resemble the gastric geometry but still skip the anatomical complexity of realistic geometries.

Gastric electrophysiology
Across the gastric wall, there is a distribution of ICC and SMC which both respond to electric signals. In particular, the ICC generate the electric signals and therefore have the ability to trigger SMC. The electrical behavior of individual cells is intrinsically coupled to adjacent cells by the transport of electric potential through dedicated proteic structures. The coupling between the electrophysiological processes at the cell level and the transport of electric potential at the tissue scale gives rise then to ordered electric patterns (emerging behavior) known as gastric slow waves at the organ level. In the following two subsections, we will first introduce a simple phenomenological two-variable model of cell electrophysiology and, subsequently, the monodomain model for the transport of electric potential across the tissue. Thereby, we provide a system of coupled partial differential equations that can reproduce the phenomenon of electric propagation of slow waves within the gastric wall.

Cell electrophysiology
To model cell electrophysiology, we adopt a phenomenological approach based on a modified version of the Mitchell-Schaeffer model, [53] originally introduced to describe the cardiac action potential, and here characterized for gastric slow wave propagation. The dynamical features of this model, delineated in the following, justify its usage for the phenomenological description of gastric electrophysiology following a wide literature on the subject. [56] The model, in fact, is capable of reproducing both pacemaker signals, needed to model ICC, and smooth muscle excitation signals, necessary to model SMC, [54] improving and stabilizing the phenomenological approach of Aliev et al. [22] . The mathematical formulation assumes that gastric electrophysiology on the level of single cells can be modeled phenomenologically both in ICC and SMC by two state variables each. These are the normalized transmembrane voltage ∈ [0, 1] and the gating (recovery) variable ℎ ∈ [0, 1], which mimics time constants and morphology of subcellular ionic currents within ICC and SMC, respectively. The index ∈ {i, m} distinguishes herein between ICC (index: i) and SMC (index: m). The evolution in time of these state variables is assumed to be governed by the system of coupled ordinary differential equations (ODEs): describes the conservation of ionic and capacitive currents across the cell membrane, a fundamental physical law of theoretical electrophysiology. [57] Ions are transported into cells by the inward transmembrane ionic current in ( , ℎ ) ≥ 0 or possibly also by some time-dependent external stimulus current stim applied to the tissue (e.g., due to some biomedical device). At the same time, charged ions are withdrawn from the cells by the outward transmembrane ionic current out ( ) ≤ 0. The conservation of ionic fluxes requires that these three currents equal rate of change of the ionic charge across the membrane (capacitive effect). Herein, all currents are normalized with respect to the capacitance of the cell membrane. Therefore the left-hand side contribution is represented by the rate of change ∕ of the transmembrane potential. From biophysical experiments, we know that the in-and outward transmembrane ionic currents depend on the transmembrane voltage . [3,6,12,13,58] To reproduce such complex behavior, we use herein simplified, FitzHugh-Nagumo-like cubic reaction functions where in and out are time constants, and represent excitability parameters. The formulation is able to reproduce ICC selfoscillation (pacemaker activity), i.e., periodic electric signals with a defined intrinsic frequency, as well as the excitable (passive) SMC dynamics via the adjustment of these excitability parameters. In Equation 1b, the voltage-dependent time constant ( ) and the voltage-dependent steady state value ℎ ∞ ( ) are defined by where open and close are time constants phenomenologically related to subcellular ionic dynamics and the parameters gate and gate modify the sign of the steady-state value. We note that this model is well-known to be valid under the assumption in ≪ out ≪ open , close . [53]

Transport of electric potential on organ scale
Electric slow waves result from an interplay between electric oscillations triggered by the ICC at each point in the tissue and a transport of electric potential on the tissue scale through cell-cell communication. We model this transport process on the domain Ω 0 ⊂ ℝ , = 1, 2, 3, in the time interval [0, ]. Let i , m ∶ Ω 0 × [0, ] → [0, 1] be the normalized transmembrane voltage of ICC and SMC, respectively. We assume that the transport of electric potential happens across both types of cells. For each cell type, it is modeled by a so-called monodomain continuum formulation (also known as cable equation) [22,59,60] according to the following reaction-diffusion system: where i , m denote the isotropic diffusion coefficients for electric potential for ICC and SMC, respectively; gap describes the homogenized resistance of gap junction proteins between the two cell types; and the normalized currents i tot , m tot are defined by the right-hand side of ODEs (1), differentiated in the parameter choices in Table 1 and assumed space-dependent. Δ denotes the Laplace operator on Ω 0 , which is linked to the Nabla operator ∇ on Ω 0 via Δ = ∇ 2 . The problem formulation is closed by suitable initial conditions, ( = 0) = 0 , ℎ ( = 0) = ℎ 0 for ∈ {i, m} (see Table 1), and with Neumann zero-flux boundary conditions (representing an insulated tissue) on the boundary Ω 0 of the domain Ω 0 with outward normal , i.e., where (⋅) denotes the scalar product.

Continuum mechanical framework
In this section, we briefly summarize the equations of continuum mechanics governing the finite-strain deformation of the stomach. The kinematics of a homogeneous continuum body, the concept of active contraction, and the associated constitutive behavior are recalled. For the sake of conciseness, the fundamental balance principles of continuum mechanics, i.e., the balance of linear and angular momentum, essential to pose the initial boundary value problem are not stated herein.

Kinematics
Let denote a material point in a body Ω 0 at time = 0 in the three-dimensional Euclidean Space, which is also called reference or material configuration. Herein, we assume that the reference configuration is load-and stress-free. At times > 0, the body may undergo a transformation, i.e., motion and deformation, to the current or spatial configuration Ω ⊂ ℝ . Hereby, the material point is related to the spatial point in the spatial (deformed) configuration by the nonlinear deformation map ∶ Ω 0 → Ω ,  → . We introduce the deformation gradient = ∕ and its determinant = det > 0, which describes the local volume change from the material to the spatial configuration.
Describing the active contraction of the gastric tissue, we follow an active-strain approach. [40] The reliability of such an approach with respect to others (e.g., active stress) has been theoretically and biophysically motivated in previous contributions. In particular, the active-strain approach can be derived from thermodynamical arguments [61] respecting the notion of distorsions in continuum mechanics. At the same time, the multiscale nature of excitable contractile cells, such as SMC, justifies the adoption of such a concept to reproduce calcium and voltage-based intra-cellular mechanisms leading to the sliding of proteic filaments (actin and myosin) generating the motion. In fact, these dynamics occur within the cell at a much smaller temporal and spatial scale with respect to our macroscopic modeling framework. Accordingly, we assume that the active deformation of the tissue due to smooth muscle contraction can be represented by a change of the local traction-free configuration of the infinitesimal volume elements of the tissue. This change of the local traction-free configuration can be represented -similar to plastic deformation in plasticity theory -by a multiplicative decomposition of the deformation gradient into a part a resulting from local activestrain deformation rather than external loading and a complementary elastic part e such that = e a .
If smooth muscle contraction, i.e., a , is known, e can be computed in a standard fashion from two conditions, first the balance of momentum and second the geometric compatibility of the deformed configuration Ω . The effect of the active deformation a is often illustrated via a so-called intermediate configuration Ω on which a maps the reference configuration and which itself is mapped on the current configuration Ω by the elastic e .
Due to the high amount of water in gastrointestinal tissue, we assume the material to be incompressible, i.e., = 1. Also, we require that both the elastic and the active part of the deformation gradient tensor fulfill such a condition, i.e., e = det e = 1, and a = det a = 1. We finally introduce the right Cauchy-Green deformation tensor as well as the so-called elastic Cauchy-Green deformation tensor

Passive material behavior
We assume a passive hyperelastic material behavior with strain energy density Ψ, which is a function of the elastic deformation gradient e : .
Though more general assumptions could be made, for simplicity and in order to focus on key aspects of electromechanical coupling herein, we assume an isotropic, incompressible neo-Hookean material model with scalar stiffness parameter and the first invariant of the elastic right Cauchy-Green tensor e 1 ∶= tr e . According to the active-strain approach, the first Piola-Kirchhoff stress tensor reads in which we do not consider additional mechano-electric feedback. In the previous relation, we have introduced the elastic first Piola-Kirchhoff stress tensor e , which describes the stresses induced by the elastic deformation e on the intermediate configuration Ω . A detailed derivation of Equation 12 is provided in the Appendix. The second Piola-Kirchhoff stress tensor and its counterpart on the intermediate configuration e can be derived as For efficient computational schemes one typically needs the tangent of the passive material ℂ e = 2 e ∕ e . We identify (8) and (13) as the pull-back of strains and stresses from the intermediate configuration into the reference configuration and apply the same pull-back operation on ℂ e . In Einstein's index notation, the resulting fourth order material tangent reads: Note that (14) only holds if the active part of the deformation gradient a is independent of the deformation .

Active material behavior
Gastric tissue typically exhibits two layers of smooth muscle with distinct fiber orientations, circumferential and longitudinal, with respect to the main geometric axis of the stomach. Therefore, modeling the active contraction of the stomach requires the introduction of an anisotropic material behavior. Let c , l be orthogonal unit vectors in the fiber directions of the circumferential and longitudinal muscle layer, respectively, in the reference configuration. n = c × l represents the unit vector in thickness direction of the gastric wall.
Allowing for active deformation in circumferential and longitudinal directions, we assume that the active part of the deformation gradient takes the following form [40] where ⊗ denotes the usual tensor product, = ( m ) ∈ [0, 1] is a smooth activation function depending on the excitation state, and c and l are material parameters controlling the intensity of the active contraction in circumferential and longitudinal direction, respectively. We introduce the additional parameter to enforce the incompressibility condition of the active deformation gradient by a stretch orthogonal to the smooth muscle layers in direction n . Note that this means that we assume that incompressibility manifests solely in a thickness change of the gastric wall. Making use of the orthogonality requirement of the fiber directions, we can define the Jacobian of the active deformation gradient as Solving (16) for results in In order to ensure that the fibers cannot contract to zero or negative lengths, we enforce each factor in the denominator of (17) to be positive. For an arbitrary ∈ [0, 1], this condition is equivalent to Note that negative values of c and l describe lengthening of the muscle fibers.

Excitation-contraction coupling
The generation of active force within gastrointestinal smooth muscle cells depends on the intracellular Ca 2+ concentration. [21] The Ca 2+ concentration itself is regulated by complex spatio-temporal dynamics involving both VDCC and intra-cellular stochastic mechanisms. [58] In order to keep the computational cost of the proposed model minimal, in this first contribution, we neglect such intra-cellular multiscale Ca 2+ -dynamics and propose a direct relation between the transmembrane voltage of the SMC, m , and the active-strain activation parameter in the form Here, t is the normalized opening voltage of the VDCC, 1 is a parameter resembling the homogenized Ca 2+ dynamics, and 2 describes the opening dynamics of the VDCC. ( ) denotes the Heaviside step function.

General
We implemented the coupling of the two-variable Mitchell-Schaeffer electrophysiology model with the active-strain electromechanical description of the gastric tissue as described in Sec. 2 in our in-house research simulation code BACI (written in C++) and we performed a series of computational tests ensuring numerical reliability of the proposed computational model. The scalar transport problem and the structural mechanics problem are solved numerically using the finite element method. The geometry is discretized by scalar transport elements and standard nonlinear membrane finite elements. Our implementation allows for non-matching meshes of both discretizations in order to reduce the overall computational costs. A membrane is considered to be a thin-walled hyperelastic structure without bending stiffness in the surface tangent plane. Plane stress is assumed in the membrane. Nevertheless, a pressure load may act in the normal direction on the membrane surface and be balanced by the in-plane stresses and membrane curvature. In general, membranes exhibit a highly nonlinear behavior due to simultaneous presence of both geometrical and material nonlinearities.

Preliminary study on a one-dimensional gastric electrophysiology model
To verify the correct implementation and the numerical properties of the proposed mathematical model, we set up a onedimensional (1D) gastric electrophysiology model, not including any mechanics yet. The one-dimensional domain with T A B L E 2 Parameters for the spatial distribution functions i ( ) of type (20) and i ( , ) of type (22 ∈ [0, ] and = 250 mm represents a line on the stomach's surface from the pacemaker region to the pylorus along the greater curvature. We apply no-flux boundary conditions at the two end points, = 0 mm and = 250 mm. In this 1D numerical analysis, we consider a time step size of = 3.125 ms and a mesh size of ℎ = 15.625 m, corresponding to 16000 two-noded linear scalar transport elements. The electrophysiology parameters for ICC and SMC together with the initial conditions are collected in Table 1.
In the gastric tissue, only the ICC exhibit the ability to actively trigger electric waves. Therefore, for the SMC the excitability parameter m = 0. For the ICC the excitability parameter, as introduced in Equation 2a, is in general unequal to zero and varies in space. Its exact spatial distribution can be expected to vary patient-specifically. Here, we assume it to be given by the smooth scalar function where i min , i max are the desired minimal and maximal excitability parameter values 1 ; is a normalized function controlling the spatial variation of the excitability in -direction. For simplicity, we herein assume the function to be of the Gaussian-like type where the index indicates the correlation between the normalized variable = ∕ to the physical variable , is the chosen reference length. In the present case, we have = . In (21), the shape parameter is linked to the variance of the underlying Gaussian and the parameter is defined such that (1) = . Note that (0) = 1 holds for all ∈ ℝ ≠0 , ∈ ℝ.
All electrophysiology parameters for ICC and SMC together with the initial conditions are collected in Table 1. The constitutive model parameters for distribution function (20) are provided in Table 1 and Table 2.

Entrainment Analysis
A key feature of gastrointestinal electrophysiology is the so-called entrainment process. This process is a form of emergent behavior. In the upper, proximal part of the stomach the excitability parameter and therefore the intrinsic frequency of electric oscillations on cell level is higher than in the lower, distal part, towards the pylorus. Therefore, if isolated from each other, the ICC in different regions of the stomach would exhibit electric oscillations with different frequencies. However, due to the coupling of the ICC by the transport of electric potential across the tissue on the continuum scale, the ICC in the lower part of the stomach are entrained. That is, the frequency of their electrostatic oscillations on cellular level is forced towards the higher intrinsic frequency of the ICC in the upper part of the stomach. The region containing the ICC with the highest instrinsic frequency, located in the upper part of the stomach, is therefore also referred to as pacemaker region.
With such a phenomenology in mind, Figure 1 shows the temporal evolution of the frequency of the oscillations of the electric potential at each point in the domain as observed in our computational model.
In the beginning, at each point the frequency is given by the intrinsic frequency of the ICC at this point, that is, determined by . Therefore, one observes a significant gradient of the frequency across the domain at = 0. Over time this gradient diminishes. At = 500 , the electric oscillations of the ICC exhibit nearly the same frequency in the whole domain, which is equal to the initial frequency of the ICC in the upper part of the domain (i.e., at = 0). That indicates that the ICC at the opposite end of the domain have been entrained to these pacemaker ICC in the upper part. According to experimental evidence, [14,15] the frequency entrainment values obtained with our phenomenological model fall within the physiological range with good accuracy. As The excitation of SMC follows the slow waves closely but with reduced amplitude (note the overall darker tone) and with a slightly delayed peak (see also Figure 4) indicated by straight horizontal lines, for > 350 the entrainment is stable over the whole domain. This result ensures that our phenomenological modeling is capable of reproducing a physiological entrainment process.
To further confirm such behavior, we provide in Figure 2 the development of the phase portrait of the two electrophysiological state variables from the intrinsic to the entrained state for both ICC and SMC. We highlight in red the last 40 s showing that both cell types develop a stable limit cycle. In addition, Figure 3 shows the space-time propagation pattern of the ICC transmembrane potential i during the last 100 s, where a clear periodic entrained state is recognized. Slow waves (grey to white areas), marked by a steep activation front, are initiated at the pacemaker site and travel along the domain, i.e., from left to right, with a constant conduction velocity. We also observe stable refractory phases at the resting potential of i = 0 (see also Figure 4). The excitation of SMC follows the slow waves generated by the ICC closely. However, the amplitude is smaller and the peak occurs slightly F I G U R E 4 Normalized transmembrane potential i of ICC (depicted in black) and normalized transmembrane potential m of SMC (depicted in blue). Solid lines depict a reference solution computed with = 3.125 ms, ℎ = 15.625 m, dashed lines depict the solution obtained with a much coarser discretization with time step and mesh size = 0.1 ms, ℎ = 0.5 mm. The relative difference is small; for the conduction velocity (CV) it remains below 5% (see also Figure 5) delayed (see also Figure 4). It is worth noticing that such robust periodic behavior is one of the major advantages with respect to previous phenomenological models. [22]

Conduction Velocity Analysis
The conduction velocity (CV) of a slow wave can be estimated from the gradient of the sharp activation lines separating bright and black regions in the space-time plot in Figure 3. At a certain distance from the pacemaker region a constant CV of 3.6 mm s −1 establishes, whereas within the pacemaker region (up to approx. = 50 mm) CV exhibits a negative gradient. Since equal boundary conditions are applied at both ends of the one-dimensional simulation domain, the computed CV variation is due to the spatial gradient of excitability nonlinearly linked to the entrainment process (see Figure 1). Averaging the conduction velocity over the whole simulation domain results in the value 3.8 mm s −1 which, again, falls within the expected pyhsiological range. [15] For illustrative purposes, Figure 4a shows the temporal evolution of the transmembrane potentials of ICC and SMC in the center of the domain, which corresponds to a vertical cut through Figure 3 at = 125 mm. Figure 4b shows the potentials over space in the entrained state (horizontal cut through Figure 3 at = 500 s). The obtained slow waves are stable and entrain to a constant frequency for both ICC and SMC. In fact, the activation of the SMC follows the pacemaker potential of the ICC with a slight delay (modeling the complex intra-cellular calcium mechanisms) and the SMC experience a lower maximal depolarization. This result, as well as the general shape of the slow waves, is in agreement with experimental recordings. [62] Numerical Convergence Analysis It is well known that the conduction velocity of reaction-diffusion systems in general depends on the numerical solution scheme. [63] Therefore, we examined the influence of the numerical parameters, namely the time step size and the mesh size ℎ, on the computed CV. We chose to evaluate the CV in the middle of the simulation domain in order to minimize the influence of the no-flux boundary conditions and after approximately = 500 s of simulated time, ensuring a complete entrainment of the system. The CV is determined by tracking the propagation of the wave front (threshold i = 0.5) for 0.6 s.
In the convergence study, we varied the time step and mesh size in the ranges 3.125 ms ≤ dt ≤ 0.2 s, 15.625 m ≤ h ≤ 1.0 mm, each subsequently by a factor of two, performing a total of 49 simulations in the selected parameter space. As a reference solution, we choose the simulation result achieved with the smallest time step and mesh size (used for the analysis above). The convergence plot in Figure 5 shows the CV for each pair of numerical parameters and ℎ. Distinctly, the CV converges with decreasing time step size and decreasing mesh size and justifies the choice of our reference solution.
The black contour line in Figure 5 indicates a 5% relative deviation of the CV from the one computed with the finest discretization in space and time (i.e., the reference solution used above). To balance computational cost and numerical accuracy, we select = 0.1 s, ℎ = 0.5 mm (with a CV relative error less than 5%, refer to Figure 5) for the solution of the electrophysiological problem in all subsequent examples. The dashed traces in Figure 4 show the solution of the model using the selected space-time discretization compared to the finest solution, depicted in solid lines. The visual comparison of the two solutions reveals that the shape of the simulated slow waves does not change significantly (see also Figure 4a). While a slight variation of CV appears, the spatial offset remains, however, basically constant.

Electro-mechanical coupling on a two-dimensional tissue patch
We continue the set of examples with a two-dimensional tissue patch incorporating an active-strain based electro-mechanical coupling. Our model domain is {( , ) ∈ [0, ] × [0, ]} with = 250 mm and = 160 mm. This domain can be imagined to resemble, for example, the half lateral surface of an idealized cylindrical stomach that has been cut-open along the greater and lesser curvature and flattened in the -plane, neglecting the electrically quiescent fundus area. The -axis is aligned with the longitudinal direction and the -axis is aligned with the former circumferential direction of the half lateral surface of the stomach. We expect structural displacements to occur in such a flattened tissue patch solely in-plane of the patch, while thickness adjustments due to the incompressibility constraints occur normal to the patch. In the reference configuration, the thickness of the membrane is set to = 3.5 mm. [64] Motivated by the previous convergence analysis, the numerical parameters are fixed to = 0.1 s, ℎ = 0.5 mm for the electrophysiological model throughout the remaining examples. For the structural mechanical part, instead, we use a mesh size of = 2.0 mm, with the purpose of reducing the computational costs. Such a choice is in line with similar studies in which a relation between electrophysiological and mechanical meshes has been investigated, setting the maximal mesh resolution for the mechanical problem between twice and eight times the size of the electrical one (see e.g. Cherubini et al. [39] or Colli Franzone et al. [65] ). Our choice results in 160000 four-noded bilinear quadrilateral scalar transport finite elements for the electrophysiological domain and 10000 four-noded bilinear quadrilateral membrane finite elements for the the structural mechanical problem.
As before, we prescribe no-flux boundary conditions on the whole boundary for the electrophysiological model. The structural boundary conditions of the planar patch, instead, are chosen to constrain displacements normal to the boundary whereas displacements tangential to the boundary are unconstrained. Namely at the vertical boundary lines ( = 0, ∈ [0, ]) and ( = , ∈ [0, ]) displacements in -direction are constrained to be zero. On the horizontal boundary lines ( ∈ [0, ], = 0) and ( ∈ [0, ], = ) displacements in -direction are constrained to be zero. These boundary conditions are chosen such that they resemble the ones to which a half-stomach can be expected to be subject in vivo. There, the plane through minor and major curvature is in good approximation a symmetry plane of the stomach. The corresponding symmetry boundary condition is applied in our model on the boundaries in -direction. At the same time, the stomach appears at both ends to be rather tethered in axial direction (by esophageal sphincter and pylorus) but -at least to some extent -unconstrained perpendicular to the axial direction. In our example, on the whole two-dimensional patch out-of-plane displacements of the membrane are constrained to be zero (mimicking the idea of a tissue patch flattened for experimental testing). For the two-dimensional computational models examined below, we assume the spatial distribution function of the excitability parameter to be described by the smooth scalar function where i min , i max are the desired minimal and maximal excitability parameter values; and are normalized functions controlling the spatial variation of the excitability in -and -direction, respectively. We assume each of these functions to be of the Gaussian-like type defined by (21) meaning we introduce two normalized variables , ∈ { , } with distinct sets of parameters (see Table 2).
In the following, we discuss a series of numerical examples with isotropic, anisotropic, and unidirectional distributions of the excitability parameter i . Constitutive model parameters are provided in Table 2 and the corresponding contour plots of the resulting distributions are depicted in Figure 6a, 6b and 7a. All remaining model parameters together with the applied initial conditions are summarized in Table 1 and Table 3. As described in Sec. 3.2, the system needs approximately 350 s to reach the entrained state. To ensure an entrained state, we study in the following subsections electro-mechanical results at = 500 s.

Isotropic spatial distribution of the excitability parameter
First, we investigate the influence of the spatial distribution of the excitability parameter i ( , ) on the two-dimensional propagation pattern of the resulting slow waves as well as on the mechanical response (see Figure 6 (left)).
We begin with an isotropic distribution of the excitability parameter i ( , ) of ICC (see Table 2 for parameter values) illustrated by the contour plot in Figure 6a. We consider the pacemaker region to be located at the origin of the coordinate system ( = 0, = 0) where we prescribe the maximum amplitude i max that isotropically reduces to i min , at a radius of 250 mm from the pacemaker site. Figure 6c shows the spatial propagation pattern of the ICC transmembrane potential i . Slow waves are initiated periodically at the origin. From there they propagate through the tissue in an almost perfect circular pattern until reaching the boundary at the top ( = 160 mm). This indicates that the CV in this set up is isotropic. Due to the influence of the no-flux boundary condition, [66] the slow waves straighten slightly once they have spread over the complete circumferential length but continue to travel in longitudinal direction whereby they heavily bend. Notably, there are always between three and four slow waves on the domain, which conforms with experimental gastric recordings. [67] Figure 6e shows patterns of mechanical activation following the electrical slow waves, implying that also in the twodimensional case the electrical activation of the SMC is synchronized with the ICC pacemaker signal. While the general spatial pattern is, as expected, very similar to that of the ICC transmembrane potential, we recognize complex mechanical deformations far from the pacemaker site (see Figure 6g). These non-trivial patterns are due to the delay of the electrical activation of SMC (see Figure 4b), the threshold behavior of the active-strain electro-mechanical coupling (see Equation 19), the nonlinearity of the adopted material model, and the effect of boundary conditions. The magnitude of the in-plane displacements, in particular, is closely linked to the mechanical activation parameter ( m ) as can be observed in Figure 6g. In general, the maximum structural contractions occur in the region of the maximum activation parameter ( m ). In addition, localized and alternating contracted and relaxed regions follow positive and negative gradients of the electrical activation waves. Note that the effect of the structural boundary conditions is clearly visible as soon as a propagating contraction wave reaches the boundary of the simulation domain.

F I G U R E 6
Collection of spatial plots illustrating the influence of different spatial distributions of the excitability parameter i ( , ) on the gastric electro-mechanical model. The left column shows the case of an isotropic distribution, the right column the case of an anisotropic distribution. The colorbars apply to both columns

Anisotropic spatial distribution of the excitability parameter
In this section, we investigate the influence of an anisotropic distribution of the excitability parameter i ( , ) comparing and contrasting the results with the isotropic case (see Figure 6 (right)). The parameters for the distribution function in x-direction remain unchanged while the distribution function in y-direction is built so as to achieve a far slower descent of the excitability parameter (see Table 2). A contour plot of the anisotropic excitability parameter i ( , ) is depicted in Figure 6b.
As expected from the preliminary 1D analysis, the decreased gradient of the excitability in y-direction results in an increased CV, yielding an elliptical shaped propagation pattern of the slow waves (see Figure 6d). Once again, under the influence of no-flux boundary conditions, the CV of incoming slow waves at the upper boundary is slightly increased in y-direction leading to a lower curvature of the shape of the wave front.

F I G U R E 7
Collection of spatial plots illustrating the effects of gastric dysrhythmias induced by a conduction block on the electro-mechanical model. The left column shows the physiological case, the right column the pathological case, where a spiral developed. The colorbars apply to both columns From the mechanical point of view, similar to the isotropic case, the spatial pattern of the mechanical activation (see Figure 6f) follows the transmembrane potential with a slight delay. However, the magnitude of the in-plane displacements illustrated in Figure 6h appears much less ordered in the presence of an anisotropic excitation. In particular, the highest magnitude of displacements is observed already on the first wave front propagating in y-direction and multiple shear-like bands appear within the tissue. These mechanical bands are of extreme importance for the correct motility and physiological function of the gastric wall.
These results confirm that the model is capable of representing different spatial electro-mechanical propagation patterns that have a large influence on the structural mechanical response.

Unidirectional distribution of the excitability parameter with conduction block
In the literature, it is has been shown, both in experiments [68] and numerical studies, [23,69] that conduction blocks can cause gastric dysrhythmias. In this section, we address such a complex phenomenon through our electro-mechanical framework, by perturbating the electrophysical activity of the tissue during the entrained state. To begin with, we revisit the unidirectional distribution of the excitability parameter of Sec. 3.2 and apply it to the two-dimensional tissue patch. Figure 7a shows the resulting contour plot of the excitability parameter . This leads to a perfectly synchronized activation along the y-axis at = 0 mm. From there, completely straight slow waves propagate in x-direction, as shown in Figure 7b. The Figures 7d and 7f show the spatial pattern of the activation parameter ( m ) and the resulting structural displacements. As expected, the mechanical activation parameter follows the straight slow waves and results in displacements occurring only in -direction with modulated amplitude along the excitability gradient.
Starting from the fully entrained system, we induce a rectangular conduction block (width: = 50 mm, height: = 80 mm, origin: = 125 mm, y = 40 mm) in the simulation domain at = 500 s in the following way: we reset the state variables , ℎ to the fixed values , ℎ for all ∈ {i, m} on the domain of the conduction block for the duration of one single time step (refer to Table 1 for values). This perturbation leads to an emerging retrograde activation wave initiated at the top edge of the former conduction block by the remaining upper half of the slow wave (not shown). The resulting electrophysiological activity is a stable rotating spiral wave, significantly altering the slow wave propagation behavior, as depicted in Figure 7c. Figure 7c shows the established spiral wave at = 545 s, and the nonlinear interaction (annihilation) with the slow waves generated by the pacemaker. An inhomogeneous conduction delay is clearly recognized from the bending of the slow wave in the lower parts while the upper part of the tissue remains less affected at this stage.
The resulting mechanical activation pattern is shown in Figure 7e and Figure 7g illustrates the magnitude of the resulting displacements. Clearly, the displacement magnitude is very elevated in the area where an approaching slow wave encounters the stable spiral. At the same time, the area covered by the spiral itself shows relatively small displacement magnitude in line with pathological scenarios, e.g. the so called paralytic ileus. [23,70] It is important to note that in a clinical setting such extreme perturbations of the displacement field could not be reconstructed nor analyzed in a straightforward manner from a simple measurement of the electrical activity in the stomach. Rather, to achieve a detailed understanding of such a situation one needs a detailed understanding of the electro-mechanical coupling in the stomach, for which in the future reliable computational multiphysics models of the stomach as the one introduced in this paper will be required. Altoghether the example in this section demonstrates that the proposed electro-mechanical computational model is able to reproduce also complex perturbation protocols and dysrhythmias in gastric tissue.

Idealized three-dimensional electro-mechanical model for gastrointestinal motility
We conclude the set of examples with a three-dimensional idealized electro-mechanical stomach model. To focus on fundamental aspects of our computational model rather than geometric details, we greatly simplify the stomach's geometry and assume a straight circular hollow cylindrical shape. We neglect the electrically quiescent fundus area that does not show the typical peristaltic contraction waves. [71] The cylinder has a length of = 250 mm, a radius of = 50.930 mm, and a membrane thickness of = 3.5 mm. The cylinder axis is aligned with the axis. Both ends of the cylinder are open. The greater curvature of the stomach corresponds in our idealized model to the line that lies within the -plane at = , whereas the lesser curvature is thought to be the line in the same plane at = − .
Motivated by the convergence analysis in Sec. 3.2, we chose a time step size of = 0.1 s and a mesh size of ℎ = 0.5 mm for the electrophysiological part of the model. To reduce the computational costs, we use a mesh size parameter of = 2.0 mm for the structural mechanical part. The above choice of numerical parameters results in a total number of 320000 four-noded bilinear quadrilateral scalar transport finite elements for the electrophysiological domain. The structural mechanical domain is discretized by 20000 four-noded bilinear quadrilateral nonlinear membrane finite elements. The thickness adjustments due to the incompressibility constraints occur normal to the membrane surface.
For the electrophysiological problem, we prescribe no-flux boundary conditions at both ends of the cylinder. For the mechanical problem, we apply zero Dirichlet boundary conditions in axial direction at both ends of the cylinder. To suppress rigid body motions, additionally the displacements of three points are constrained appropriately. Namely at ( , , ) = ( , 0, ) displacements in -direction are constrained and at ( , , ) = ( , ± , 0) displacements in -direction. Additionally, we prescribe a pressure load of = 3.0 mm −2 on the inner surface of the cylinder. [72,73] The values of all model parameters together with suitable initial conditions are summarized in Table 3 and Table 1. Approximately 350 s of simulation time are needed to ensure the system's entrainment (see Sec. 3.2). Therefore, we evaluate the results well after this point.
To ensure the rapid formation of ring-shaped contractions as observed experimentally in the stomach, [74] an increased circumferential CV is essential. We therefore resort to the anisotropic excitability parameter distribution (see Table 2). We map the F I G U R E 8 Idealized cylindrical electro-mechanical stomach model at different points of simulation time. The left column shows the normalized ICC transmembrane potential i . The magnitude of the displacements is depicted in the right column together with the resulting deformation. The colorbars apply to the respective column, where i is dimensionless and the displacement magnitude is depicted in mm planar distribution function onto the cylinder surface such that the pacemaker site with the highest intrinsic frequency is located at ( , , ) = (0, 0, ). As described in Sec. 3.3.2, we expect this to result in the desired anisotropic conduction behavior.
The left column of Figure 8 shows the spatial distribution of the transmembrane voltage i at different points in time. At = 475 s the pacemaker site (top, left) is in the refractory phase just before activation. Once initiated at the pacemaker site, the slow waves propagate both axially and symmetrically in radial direction ( = 475). As desired, the CV in radial direction is higher than in axial direction ( = 480). Therefore, the slow waves form a closed ring within approximately 20 s after initiation (see previous wave in ( = 470). The closed slow wave ring travels in axial direction until it phases out at the right boundary, = . This slow wave propagation pattern is stable. The slow waves are generated with a frequency of approximately 3 cpm, resulting in a spatial distance between subsequent waves of 40 mm. Both values are within the experimentally established physiological range. [74] The right column of Figure 8 visualizes the resulting deformation. Compared to its load-free configuration, the cylinder is inflated by the applied inner pressure to a radius of = 70.485 mm. In this deformed state the electric slow waves induce, via F I G U R E 9 Solution with increased values of l = c = 0.5. For comparison, the black line shows the outline of the solution with l = c = 0.25 in the -plane (grey). The colorbar shows the displacement magnitude in mm. By adjusting l and c , the model is capable of representing different intensities of the contractions. The contraction waves become more pronounced for higher l and c values the resulting muscular contraction, ring-shaped contraction waves which propagate axially in a synchronized and coordinated pattern. We investigated the influence of the contractility parameters l and c on the mechanical deformation. As depicted in Figure 9, the model is capable of representing different intensities of the peristaltic gastric contractions by adjusting the values of l and c . The higher these values are, the more pronounced the resulting contraction waves become. This example clearly demonstrates that our electro-mechanical model is in principle able to reproduce the essential features of peristaltic contraction waves as observed in the stomach.

CONCLUSION
In this paper, we propose a computational model of the electro-mechanical coupling in the stomach based on an active-strain finite elasticity approach. Within this framework, gastric electrophysiology is represented by a simplified phenomenological model relying on two internal state variables considering two coupled modified Mitchell-Schaeffer models. Using a series of computational examples, we demonstrated that this multiphysics approach is sufficient to capture essential phenomena of gastric electro-mechanics such as slow wave generation and entrainment across the organ, gastric dysrhythmias as well as the propagation of stable ring-shaped peristaltic contraction waves along the stomach.
Our computational analyses reveal that the spatial propagation patterns of electric and mechanic activity in the gastrointestinal tissue strongly depend on the underlying distributions of the excitability parameter which controls the intrinsic frequencies at which ICC trigger electric signals within SMC (in the absence of any entrainment due to coupling with surrounding ICC).
As a novel contribution, the present study links gastric electrophysiology and dysrhythmias to the mechanical properties of the stomach highlighting the strong mechanical nonlinearities that characterize the gastric wall. Though several experimental evidences have shown the tight coupling between dysrhythmias and gastroparesis (e.g., irregular initiation, aberrant conduction, low amplitude etc. [75] ), this link remains largely unexplored so far from a modeling and computational viewpoint and we aim at stimulating further investigations in the future.
Altogether, the computational model of gastric electro-mechanics presented in this paper can be understood as a computationally efficient and, at the same time, robust workhorse for future computational studies of the role of the electro-mechanical coupling in the stomach in health and disease.
Nevertheless, several limitations of the proposed model should be mentioned here. In our model the problem of electrophysiological wave generation and propagation is solved on the reference configuration under incompressibility constraints. This implicitly assumes a one-way coupling between electrophysiology and mechanics, neglecting any feedback from the mechanical deformation on the electrophysiology, known as mechano-electric feedback (MEF). [76] Such feedback might, however, play an important role in reality, e.g., via stretch-activated currents [77,78] or nonlinear [79,80] and stress-assisted diffusion. [42] We aim to include the effects of MEF in a future contribution generalizing the present theoretical and computational framework and to validate it against in vivo electro-mechanical data. [73,81,82] Another rather obvious limitation of our model is the purely phenomenological approach to model the spatio-temporal dynamics of cell electrophysiology. In fact, incorporating a detailed multiscale, biophysical description of the intra-cellular processes governing gastric electrophysiology might be a valuable goal for future work. In particular, we aim at validating the present active-strain electro-mechanical framework against experimental evidences in terms of biophysical electrophysiological couplings. [13,83,84] From the computational point of view, a deeper investigation of the CV in the circumferential direction using a more realistic gastric geometry should be accomplished, considering also the presence of the fundus, which is electrically quiescent but responsible for other, specific mechanical effects, namely the storage function of the stomach. Also, more efficient numerical schemes [85][86][87] and GPU-based codes [88,89] are foreseen to speed-up the wide numerical analyses necessary to fully characterize a complex organ like the stomach.
Finally, also a more realistic mechanical model of the gastric tissue should be employed. It is well-known that the gastric tissue exhibits inhomogeneous and in general also anisotropic mechanical properties, which may also strongly depend on the age of an individual. [90][91][92][93][94] Exploring the precise effect of the mechanical properties on the electro-mechanical coupling might be a key to better understand the origin of certain gastric pathologies. Also the effect of the amplitude of muscular contraction should be studied in more depth in the future. Overall, the model results should also be validated against electro-mechanical data sets in order to expand the model's predictive significance.
Addressing all these and in fact several other open points proposes a host of interesting avenues of future research in the still emerging field of gastric electro-mechanics. FOR ACTIVE-STRAIN ELECTROMECHANICS The first term on the right-hand side of (12) can be identified as the first Piola-Kirchhoff stresses that arise form the elastic deformation e in the intermediate configuration Ω a . We denote this tensor as e . For the derivation of the second term, we switch to Einstein's index notation: Neglecting mechano-electrical feedback, the active deformation gradient a does not depend on the the deformation , hence the derivative