Engaging biological oscillators through second messenger pathways permits emergence of a robust gastric slow-wave during peristalsis

Peristalsis, the coordinated contraction—relaxation of the muscles of the stomach is important for normal gastric motility and is impaired in motility disorders. Coordinated electrical depolarizations that originate and propagate within a network of interconnected layers of interstitial cells of Cajal (ICC) and smooth muscle (SM) cells of the stomach wall as a slow-wave, underly peristalsis. Normally, the gastric slow-wave oscillates with a single period and uniform rostrocaudal lag, exhibiting network entrainment. Understanding of the integrative role of neurotransmission and intercellular coupling in the propagation of an entrained gastric slow-wave, important for understanding motility disorders, however, remains incomplete. Using a computational framework constituted of a novel gastric motility network (GMN) model we address the hypothesis that engaging biological oscillators (i.e., ICCs) by constitutive gap junction coupling mechanisms and enteric neural innervation activated signals can confer a robust entrained gastric slow-wave. We demonstrate that while a decreasing enteric neural innervation gradient that modulates the intracellular IP3 concentration in the ICCs can guide the aboral slow-wave propagation essential for peristalsis, engaging ICCs by recruiting the exchange of second messengers (inositol trisphosphate (IP3) and Ca2+) ensures a robust entrained longitudinal slow-wave, even in the presence of biological variability in electrical coupling strengths. Our GMN with the distinct intercellular coupling in conjunction with the intracellular feedback pathways and a rostrocaudal enteric neural innervation gradient allows gastric slow waves to oscillate with a moderate range of frequencies and to propagate with a broad range of velocities, thus preventing decoupling observed in motility disorders. Overall, the findings provide a mechanistic explanation for the emergence of decoupled slow waves associated with motility impairments of the stomach, offer directions for future experiments and theoretical work, and can potentially aid in the design of new interventional pharmacological and neuromodulation device treatments for addressing gastric motility disorders.

Introduction Gastric peristalsis, the coordinated contraction and relaxation of the muscles of the stomach, is a critical phenomenon for food propulsion and waste product elimination [1] and is impaired in motility disorders [2][3][4]. The coordination of the contractions along the rostrocaudal compartments of the stomach causes rhythmic longitudinally travelling aboral muscle contractions [5]. Sometimes in motility disorders, the contractions can occur at a faster rate than the usual (tachygastria) or at a slower rate (bradygastria) [4,6]; whereas in some cases the rostrocaudal coordination of contractions is lost resulting in decoupling (the activity of the caudal end and the rostral end become independent of each other) or functional uncoupling (the activity of the caudal end controlling the activity of the rostral end [7,8]) of contractions in the stomach [9]. Events like these have been associated with gastric motility disorders such as gastroparesis (delay in food transit), functional dyspepsia and gastroesophageal reflux disease [4,6,10].
Peristalsis is governed by diverse and overlapping mechanisms which normally ensure robust motility patterns. Peristalsis emerges from a mutually coupled chain of pacemaker cells called the interstitial cells of Cajal (ICCs). The ICCs can independently generate electrical activity known as pacemaker potentials, which drive rhythmic potentials in the circular and longitudinal smooth muscle (SM) cells embedded in the stomach wall [5,[11][12][13]. The electrical activity propagates through this network of cells in a coordinated manner resulting in gastric slow-wave (GSW) propagation [13,14] that underlies peristalsis.
The propagation is enabled in part by gap junction channels between the ICCs, between ICCs and SM cells, as well as between SM cells [15][16][17]. The SM cells, lacking intrinsic pacemaker capability, do not regenerate the slow-wave. Regeneration of the gastric slow-wave instead occurs within the pacemaker ICCs. Without coordinated ICC-ICC interactions, gastric slow-wave events do not propagate long distances [18]. Therefore, understanding how coordination arises between ICCs for propagating the gastric slow-wave is critical for understanding abnormal peristalsis that is present in gastric motility disorders.
A causal relationship between electrical gap-junction coupling in a network of ICCs and the emergent gastric slow-wave has been unequivocally demonstrated by computational studies (e.g., [7,8,19]). Pacemaker ICCs contain calcium (Ca 2+ ) and inositol 1,4,5-trisphosphate (IP 3 ) within their cytoplasm [20,21]. Gap junctions also facilitate exchange of second messengers such as Ca 2+ and IP 3 between adjacent cells [22][23][24][25]. Such exchange of Ca 2+ and IP 3 can impact intracellular concentration of Ca 2+ and IP 3 within pacemaker cells and therefore, can modulate the oscillation frequency of the concerned cell. In addition, the intracellular concentration of IP 3 is modulated by the enteric neural innervation received by the ICC [26][27][28][29]. Therefore, understanding the role of second messengers and their neural control is likely important for understanding aberrations in slow-wave propagation that result in the emergence of gastric motility disorders.
The longitudinally arranged ICCs along the stomach wall with intercellular gap junctions, resemble a chain of coupled oscillators [30,31]. In such a chain, the rostral or leading oscillator can gradually engage the trailing oscillators along the chain such that, at a steady-state, all the oscillators are frequency and phase-locked, with a constant phase lag between consecutive oscillators [32][33][34]. This phenomenon known as oscillator entrainment is thought to be the basis for the GSW in the stomach [30,31]. The integrative role of neurotransmission and intercellular coupling mechanisms in GSW entrainment remain incompletely understood. To distinguish the intercellular exchange of second messengers from passive electrical conductance is empirically challenging due to lack of pharmacological or genetic tools. The second messengers, including Ca 2+ and IP 3 that can move passively via gap junctions [22][23][24]35,36], are involved in intercellular communication via a spatiotemporal spread of coordinated oscillations in gap-junction coupled networks [22][23][24][35][36][37] and can increase the GSW frequency [38]. Yet, it is not clear whether their intercellular gap junction permeabilities can enable regenerative movement from the leading to the trailing ICCs to entrain cells and thus facilitate propagation of the slow-wave [8,18].
Our objectives were therefore to first develop a computational framework which includes constitutive stomach wall cells, biophysical models of the cells, gap junctions through which electrical coupling exists, second messenger exchange across gap junctions, and modulation of second messenger concentration by an endogenous enteric neural innervation. Second, we utilized this framework to assess the contribution of intercellular electrical coupling and intercellular exchange of second messengers on longitudinal entrainment of the gastric slow-wave.
To develop a computational framework, we created a gastric motility network (GMN). We modeled realistic cell models for ICCs incorporating the intrinsic mechanisms for pacemaker potential generation (or oscillatory) activity and intracellular IP 3 dynamics. To enable inter-ICC coordination, we modeled electrical coupling and second messenger (IP 3 and Ca 2+ ) exchange between adjacent cells in a chain. Each ICC was also coupled to an SM cell, forming an oscillatory unit (OU). The membrane voltage and period (or frequency) of the SM cells were used to assess network entrainment, as in empirical studies. Further, we set a rostrocaudal gradient for the enteric neural stimulus input to the ICCs along the ICC chain.
Since enteric neural innervation influences the intracellular second messenger concentration, we first evaluated the importance of a gradient in the innervation along the chain for development of a caudally propagating gastric slow-wave. We hypothesized that a linear gradient would enable gastric slow-wave propagation along the chain. Next, we examined whether exchange of second messengers alone can generate propagation of the gastric slow-wave. We hypothesized that addition of exchange of second messengers to electrical coupling would enhance the length of the oscillator chain over which entrainment is preserved (referred to as entrainment range hereafter). When the oscillators spanning the rostral end of the chain to the caudal end of the chain are entrained, we refer to the phenomena as rostrocaudal entrainment and if the entrained region spans a certain length from the rostral end, but falls short of the caudal end, we considered it partial entrainment. We also assessed how changing electrical coupling strength and second messenger permeability influenced the pacemaker potential frequency and the time for the slow-wave to propagate from the rostral to the caudal end of the stomach, i.e. the velocity of the slow-wave. We hypothesized that systematic increase in coupling strengths would both increase the pacemaker potential frequency and reduce the time taken for gastric slow-wave propagation. Finally, we hypothesized that variability in coupling strengths from cell-to-cell would not disrupt longitudinal entrainment if both electrical coupling and second messenger exchange mechanisms are present.
A novel computational framework with both pacemaker and muscle cells, that includes intracellular second-messenger dynamics in pacemaker ICCs, intercellular exchange of second messengers between ICC, and where the pacemakers receive a neural innervation was developed to simulate and understand propagation of the slow-wave in the stomach. The results indicate that presence of a gradient in neural modulation along the ICCs supported a rostrocaudal propagation of an entrained gastric slow-wave. Although intercellular exchange of second messengers is not necessary for slow-wave propagation, its presence can enhance the rostrocaudal length of the stomach over which entrainment is preserved and, in its absence, the entrainment range is compromised (partial entrainment is observed). This compromise is reflected by signs of decoupled slow waves and bradygastria. As hypothesized, on an increase of electrical coupling strength and second messenger permeability, the velocity of slow-wave propagation increased while the pacemaker potential frequency increased with the former and decreased with the latter. Our model with the distinct intercellular mechanisms (exchange of second messengers and electrical coupling) in combination with the intracellular feedback pathways and a rostrocaudal neural innervation gradient allows SM cells to oscillate with a moderate range of frequencies and the gastric slow-wave to propagate with a broad range of velocities. Importantly, in the presence of variability in coupling strengths as would occur in biological networks, the existence of intercellular exchange of IP 3 can preserve the longitudinal entrainment to a greater extent along the length of the stomach and eliminate signs of bradygastria and/or tachygastria. Together these results enhance our understanding of the intrinsic and extrinsic mechanisms engaging second messengers for the propagation of a robust gastric slow-wave essential for normal peristalsis.

A biologically realistic gastric motility network (GMN) model for the stomach
To develop a gastric motility network (GMN) model for the stomach, we considered the length of the stomach spanning the mid-corpus to the terminal antrum (a length of approximately 6 cm in the greater curvature [39]). The schematic of the stomach in Fig 1A highlights the arrangement and diversity of cell-types in the stomach wall [11,12]. We focused on modeling the myenteric ICC and circular muscle layers which are most widely studied experimentally. The computational models for the ICCs and SM cells include a diverse set of biophysical properties reported in these cells [40][41][42] and are formulated as conductance-based models [43,44]. In the pacemaker ICCs, we incorporated widely accepted pathways for producing intrinsic oscillatory behavior. These included a cytosolic mitochondrial-endoplasmic reticulum (ER)based Ca 2+ buffering mechanism, a membrane potential dependent intracellular concentration change in IP 3 , and IP 3 receptor mediated Ca 2+ release from the ER, in addition to the transmembrane voltage-activated ionic currents (see Eqs 3 and 4 in Methods) as shown in Fig 1B  [ 20,42]. Whereas for the SM cell model, we included a muscle-specific sarcoplasmic-reticular mechanism of Na + /Ca 2+ exchange to moderate the intracellular Ca 2+ in addition to the transmembrane currents as noted in Fig 1C. The ionic current equations defining the properties of these cells were derived from published models for the ICCs and SM cells (e.g., [43,44]) and rely on relevant experimental work for the validity of the underlying model assumptions (e.g., [40,41,45]). These cellular features were incorporated in our model to ensure both biological compliance as well as agreement with published computational models (e.g., [7,19,43,44] see S1 Table). Our equations, parameter descriptions, and simulation methods are described in the Methods section. The resulting simulated membrane voltage characteristics bear close resemblance to empirical evidence, as shown in Fig 1D ( The intrinsic frequency of oscillations for both ICCs and SM  cells were tuned to~3.0 cycles per minute to match closely with the membrane potential  recordings from the guinea-pig gastric antrum ICCs (see Fig 2A for ICC recording and Fig 2B from SM cell recording in [5]). Specific preparations to which our results are comparable are provided in S2 Table. The gastric motility network architecture consists of a longitudinal chain of ICCs bordered by a similar inner chain of SM cells, with one-to-one connectivity formed through nearestneighbor coupling as shown in Fig 2. The non-boundary cells are connected to two adjacent cells, whereas the boundary cells are connected to only one adjacent cell (see Eq 5 for ICCs and Eq 10 for SM cells). This closely mimics their biological arrangement in the mammalian stomach [14,39,40]. Endogenously ICCs are coupled via gap junctions [15][16][17] which permit exchange of ions and other small molecules. To develop a reasonable quantitative model for gap junction coupling involving multiple mechanisms, we searched for experimental evidence for two essential aspects of coupling: 1) density of gap junctions along the gastrointestinal (GI) tract reflecting the rostro-caudal distribution of coupling strengths, and, 2) knowledge of which molecules are predominantly exchanged during entrainment. Neither of these two pieces of information are currently clarified for ICCs. Therefore, we made biologically plausible assumptions based on available experimental evidence for gap junction coupling from ICCs as well as cellular networks which are composed of similar gap junction channels. As such, we incorporated electrical conductance [15] and diffusive exchange of Ca 2+ and IP 3 [22,23] and tested both equal permeabilities and randomly distributed gap junction permeabilities. Between each ICC and its connected SM cell, and between adjacent SM cells we incorporated electrical conductance-based coupling [15]. For the overall network, we incorporated 42 ICCs and 42 SM cells (offering a spatial resolution of approximately 7 cells/cm) where the intrinsic frequency of the first cell was 3.6 cycles per minute (cpm) and the last cell was 2.1 cpm. Each OU consisting of one ICC and one SM cell was modeled using ordinary differential equations (ODEs) with 37 state variables; there are 131 algebraic equations for relevant nonlinear terms in those ODEs (92 for ICC, 38 for SM cell, and one for the equation describing coupling between ICC and SM cell). The total number of parameters involved in a single OU simulation was 135 (111 for ICC and 24 for SM cell). As such the network offers a complex, but realistic, framework to evaluate mechanisms involving ICC-ICC, SM-SM, and ICC-SM coupling and their contributions to inter-ICC coordination resulting in the GSW propagation. Fig 2 highlights the framework used to assess whether and how the intrinsic and extrinsic mechanisms engaging second messengers control the inter-ICC coordination necessary for generation of the GSW and thus, the entrainment range in this GMN. The enlarged box inset shows the intracellular pathways of the ICC model essential for its intrinsic oscillatory behavior or pacemaker activity. These consist of the multi-stage feedback pathway between the membrane potential, V m and the intracellular Ca 2+ and IP 3 concentrations. The intracellular Ca 2+ is further divided into two compartments: Ca 2+ in the submembrane space (SS) spanning a narrow region between the ER and mitochondria ([Ca] SS,i ) and near membrane Ca 2+ ([Ca] i ). The bold arrows in the schematic show the inflows and outflows of the key model variables (see Table 1 for description of the membrane ion channel currents). The dashed arrows with annotation highlight the positive feedback pathways important for intrinsic pacemaking in the ICC model. Note that the membrane potential exerts a positive feedback on intracellular IP 3 concentration (see V m −IP 3 feedback in Fig 2, Eq 4 in Methods, also [20,27,42]); The intracellular IP 3 concentration, [IP 3 ] i , in turn affects the [Ca] SS,i concentration through the outward Ca 2 + flux from the ER (J Erout ) [46][47][48], as highlighted by the IP 3 −Ca 2+ feedback in Fig 2. Lastly, the [Ca] SS,i increases a non-selective cation channel current (NSCC) which closes the loop by impacting the membrane potential (the Ca 2+ −V m feedback) [49,50].
An initiating event for pacemaker activity (enteric neural innervation) is assumed to increase the IP 3 production rate and is modeled as a constant rate, β, for each ICC (see Eq 4 in Methods). Biologically, increases in IP 3 levels can originate from endogenous release of There is a negative gradient in enteric neural innervation (β) that modulates IP 3 production rate along the rostrocaudal ICC chain (ICC index). Key variables and their interconnections impacting the functionality of an ICC are illustrated in the enlarged inset. The bold arrows in the inset show the inflows and outflows of key variables. The dashed arrows indicate dependency of the sink variable on the source variable. The dashed arrows with (+) sign highlight the positive feedback pathways important for intrinsic pacemaking by an ICC. See Table 1

PLOS COMPUTATIONAL BIOLOGY
Emergence of robust gastric slow-wave during peristalsis neurotransmitters from enteric neurons [26] and subsequent activation of G-protein coupled muscarinic receptors on ICCs at cholinergic synapses along the GI tract [26,51]. In our model, we assume the neurotransmitter release and uptake by receptors as a single event termed as 'enteric neural stimulus'. The enteric neural stimulus driven increase in [IP 3 ] results in an increase in the [Ca] SS,i mimicking the endogenous release of Ca 2+ from IP 3 -operated stores in the ER. In response to [Ca] SS,i increase, the Ca 2+ uniporter on the mitochondrial membrane is gated open, and Ca 2+ ions flow into the mitochondria down the steep electrochemical gradient (the term J Uni in Eq 3 in Methods represents this mechanism). This is thought to remove a larger number of Ca 2+ ions from the submembrane space than had previously entered from the ER, causing a temporary drop in the subspace Ca 2+ concentration [44,49,50]. This activates the non-selective cation channel current leading to membrane depolarization and onset of oscillation in the ICC. Subsequently, further increase in [Ca] i due to opening of voltage-dependent Ca 2+ channels is followed by activation of voltage-and Ca 2+dependent K + currents to cause repolarization that restores V m to hyperpolarized values. The above events repeat to cause the regenerative pacemaker potential activity.
To enable longitudinal entrainment of ICCs along the stomach's length with a rostrocaudal frequency gradient we set a rostrocaudal gradient for the enteric neural innervation (β) that modulates IP 3 production along the ICC chain (see Fig 2). Such a gradient reflects evidence that cholinergic inputs to the ICCs show a rostrocaudal decrement along the GI tract [52,53].
The overall framework allowed us to examine whether engaging adjacent ICCs via two distinct mechanisms would result in similar entrained slow-wave. One mechanism relies on the electrical conductivity of the gap junction coupling between ICCs which directly depolarizes their membrane voltage during entrainment. This is a widely used formalism used in previous modeling studies [7,19]. The second mechanism based on exchange of IP 3 and Ca 2+ between adjacent ICCs has not been explored in previous models to examine the entrainment range (however see [8]). In our GMN model we also assume that exchange of second messengers between ICCs can occur via gap junctions [22,23,25]; the permeabilities are set for these small molecules (see Eqs 6 and 7 in Methods).

Longitudinal entrainment of ICCs produces a slow-wave for normal peristalsis
We tested whether the GMN model can produce a slow-wave of uniform frequency similar to a rostrocaudally propagating slow-wave for normal peristalsis in the intact stomach (e.g., in the cat [39], in the dog [40], and guinea pig and humans [14]). We also examined whether the SM cells driven by the entrained ICCs exhibit the expected positive and constant time lags between adjacent cells from the rostral to the caudal end of the network [30,33,39,54]. The ICC/SM pair function as an oscillatory unit. We plot the key state variables of such an OU in  3 concentration is relatively gradual. Although each OU has its own intrinsic frequency of oscillation, OU intrinsic frequency is linearly dependent on the enteric neural innervation (see S2 Fig). We report the membrane potential of SM cells as a proxy for the OU activity (also see S4 Fig) as is done in empirical studies [39,55]. The network was simulated for 900 seconds (Fig 3), where entrainment was observed approximately after 300 seconds of transient response. A longer 4000-sec simulation also confirmed entrainment.  Fig 3E). These results indicate the entrainment of all the oscillatory units to a uniform gastric slowwave frequency.
For comparison with empirical results, we also performed meta-analysis of in vitro recordings from the cat stomach reported by Xue S et al., [39] and generated the Relative Lag and spatiotemporal map of SM Cell Periods. These are illustrated as insets in Fig 3D and 3E and corroborate our model findings. The spacing between two consecutive Relative Lags relates to the slow-wave velocity in that particular region. Fig 3D shows that the Relative Lag is small for proximal cells and increases along the network. Interestingly, this spacing between increasing Relative Lag is nonlinear. The inset, which is derived from an experimental recording, supports the simulated observations. Thus, we demonstrate that our GMN model is a biologically plausible comprehensive network capable of generating an entrained slow-wave in the stomach with the appropriate rostrocaudal phase lags to support normal peristalsis.

Inter-ICC electrical coupling and second messenger exchange synergize to generate slow-wave propagation
Our in silico model with distinct electrical gap junction conductance and IP 3 and Ca 2+ permeabilities enabled us to study their specific contributions to network entrainment and GSW properties as shown in Fig 4. We examined the behavior of the in silico GMN model with  Fig 3D) and SM Cell Period (Fig 3E) were generated from the in vitro recordings of SM activity at different locations along the length of the cat stomach [45], where location 1 is the rostral-most recording.
https://doi.org/10.1371/journal.pcbi.1009644.g003 coupling between ICCs through electrical gap junctions alone (Fig 4A1), with coupling through second messenger exchange of Ca 2+ and IP 3 (Fig 4A2) alone, or with both ( Fig 4A3) over 900 seconds (as for the default network in Fig 3). When either electrical gap junction coupling or Ca 2+ and IP 3 permeabilities alone were present, we observed partial entrainment along the rostrocaudal chain as demonstrated by the Relative Lags (Cell #1-33 entrained in Fig  4B1 and cell #1-28 entrained in Fig 4B2). The SM Cell Period diagram in Fig 4C1 and 4C2 also supports this observation. Fig 4C1 shows that for the network having electrical coupling alone, cell #1-33 of the network oscillate with a period of~16.46 sec and no uniformity in cell period is observed from cell #34-42. For the network having exchange of second messenger only, Fig  4C2 shows that cell #1-28 oscillate with a period of~17.8 sec and afterwards, no uniformity in cell period is observed. The results with electrical gap junction coupling alone (Fig 4B1 and  4C1) were akin to those observed in motility disorders where slow waves are decoupled and bradygastria is observed in the antrum while the corpus continues to demonstrate normal electrical activity [56]. However, when both constitutive mechanisms were activated simultaneously, an enhanced entrainment range (Fig 4A3, 4B3 and 4C3) was obtained suggesting that addition of exchange of second messengers to electrical coupling has a synergistic effect.

Impact of increasing strengths of electrical gap junction coupling and exchange of second messengers between ICCs
We performed a sensitivity analysis to characterize the effects of increasing strengths of electrical conductivity versus second-messenger permeabilities in inter-ICC coupling on network entrainment and found that they influence the entrainment range, the gastric slow-wave velocity, and the pacemaking frequency. Fig 5A-5C illustrate the network's behavior for increasing inter-ICC electrical coupling conductance, G ICC−ICC two times above and below a default value of 0.7 nS. The spatiotemporal maps of SM cell membrane potential and the corresponding Relative Lag and action potential periods of SM cells are respectively shown for the lowest (Fig  5B1-3; left) and highest (Fig 5B1-3; right) G ICC−ICC values considered. In this range of coupling conductance strengths, the network transitioned from partial entrainment (Fig 5B1-3; left) to complete rostrocaudal entrainment (Fig 5B1-3; right). Here, even though cell #1-33 are entrained, decoupling is observed from cell #34-42. Unlike the network having lowest G ICC−ICC value, the network having the lowest P IP3 value considered demonstrates a reduced cell period of~17.3 sec in the entrained region composed of cell #1-33 (illustrated by the SM Cell Period diagram in Fig 5F3, left). This increased SM cell frequency could be an underlying sign of tachygastria. The increase in SM cell frequency and decoupling can be prevented by increasing exchange of IP 3 across ICCs (Fig 5F3, right) [57].
For increasing G ICC−ICC and P IP3 , in both cases, the effect on Total Lag was qualitatively similar, i.e., decreased with increasing coupling strengths and reached an asymptote as shown in Fig 5D1 and 5H1 (actual values are provided in S3 and S5 Tables respectively). These results suggest that increasing inter-ICC coupling strengths leads to robust rostrocaudal network entrainment wherein the gastric slow-wave velocity approaches a limit.
However, G ICC−ICC and P IP3 increases differentially altered the SM Cell Periods in the entrained networks. This is shown in Fig 5D2 and 5H2 respectively (actual values in S4 and S6 Tables). The SM Cell Period decreased with increasing G ICC−ICC , whereas increasing values of P IP3 increased the SM Cell Period. These results suggest that with an increase in P IP3 , a decrease in the pacemaking frequency occurs whereas the gastric slow-wave velocity increases. Collectively, the results suggest that on an increase of electrical coupling strength and IP 3 permeability, the velocity of gastric slow-wave propagation increased while the pacemaker potential frequency increased with the former and decreased with the latter. We also exclusively increased the Ca 2+ permeability (P Ca 2þ ) to test its effects on network entrainment. We noted that even two orders of magnitude changes above and below the default value of P Ca 2þ produced negligible effect on all the network characteristics analyzed: entrainment range, Relative Lag, SM Cell Periods and Total Lag in the entrained networks (see S3 Fig). This is likely due to an order of magnitude lower permeability of Ca 2+ compared to IP 3 for gap junctions [18]. Due to this lack of effect of P Ca 2þ on entrainment and gastric slowwave characteristics, in what follows, we will only focus on changing P IP3 .

Inter-ICC IP 3 exchange preserves longitudinal entrainment
Variability in ICC coupling strengths could stem from variable gap junction densities in the stomach [58,59]. We find that inter-ICC IP 3 ensures appropriate slow-wave propagation and stable pacemaking potential generation under variable ICC coupling strengths. We first introduced variability in coupling conductance by random sampling of G ICC−ICC values from a uniform distribution (see corresponding results in Fig 6A-6F) in the absence of inter-ICC IP 3 exchange. The range of conductance values were default ± 20% or ± 50% or ± 100% (see Fig  6A). As shown in Fig 6B, the Total Lags have not reached a constant value in any of the three networks. Consequently, slow-wave velocity cannot be inferred from here and the failure of Total Lag to attain a constant value indicates absence of rostrocaudal network entrainment. Fig 6C shows the distribution of Total Lags computed over the last 7 cycles of simulation, using violin plots. From these plots we note that the variability in Total Lags consistently increases with the variability in coupling conductance strength. The SM cell membrane potential (Fig 6D), Relative Lag, (Fig 6E) and Period (Fig 6F) indicate that variability of G ICC−ICC indeed gave rise to clusters of partial entrainment. The Relative Lag diagrams in Fig 6E, Fig 3). However, the non-entrained cells in Fig 6F1-3 show a predominantly longer cell period. Cumulatively, these changes in cell period are potential signs of tachygastria and/or bradygastria. These results suggest that in the absence of inter-ICC IP 3 exchange, decoupled slow waves and potential signs of bradygastria and/or tachygastria can coexist.
Next, we included the inter-ICC IP 3 exchange along with electrical coupling and varied both the G ICC−ICC and P IP3 values 20%, 50%, and 100% about their corresponding default values of 0.7 nS and 8.0 sec -1 , respectively (Fig 6G). Interestingly, addition of IP 3 exchange enabled rostrocaudal entrainment for networks with 20% and 50% variability in these coupling parameters, but not with 100% variability (see panels in Fig 6H). The distribution of Total Lags in Fig 6I show  The corresponding values of the slow-wave velocity were~0.28 cm/sec and~0.27 cm/sec, respectively. At 100% variability, the network was only partially entrained. Hence slow-wave velocity could not be measured. The SM cell membrane potential (Fig 6J), Relative Lag ( Fig  6K), and Period ( Fig 6L) further confirm these results. The Relative Lag diagrams in Fig 6K1  and 6K2 show absence of decoupling for 20% and 50% variability, respectively (the entire network from cell #1-42 is entrained). In contrast, for 100% variability, Fig 6K3 shows decoupled slow waves as evidenced by temporal changes in Relative Lag (cell #1-26 are entrained, the rest are decoupled). Similarly, the SM Cell Period diagrams in Fig 6L1 and 6L2 indicate elimination of signs of abnormal rhythmicity. However, for 100% variability, signs of increase in SM Cell Period and decoupling persist (Fig 6L3) in the non-entrained caudal cells (cell #27-42). Furthermore, Fig 6H demonstrates that there is no substantial difference in the latency to reach entrainment between the networks having 20% and 50% variability. Together these results suggest that the existence of intercellular exchange of IP 3 can preserve the longitudinal entrainment to a greater extent along the length of the stomach, thereby preventing decoupling and restoring the normal behavior.

Discussion
We have developed a computational framework consisting of a novel non-linear mathematical model for slow-wave propagation in the stomach wall that includes physiologically established intra-and intercellular mechanisms. Utilizing this framework, we have assessed the contribution of intercellular electrical coupling and intercellular exchange of second messengers on longitudinal entrainment of the gastric slow-wave. The intracellular concentrations of Ca 2+ and IP 3 are modulated by intercellular exchange of respective molecules (see Eqs 2 and 4). We find that our model with dynamically coupled nonlinear oscillators fed with second messengers from intercellular exchange and enteric neural innervation, can regulate the frequency of contractions and velocity of the slow-wave, and can enhance the range of longitudinal entrainment of the gastric slow-wave. The combination of electrical coupling and exchange of second messengers provides robustness to the entrainment range in the presence of biological variability. In summary, our detailed analyses of the ICC coupling mechanisms reveal ways in which the constitutive gap junction coupling mechanisms can enhance the network range and stability of the gastric slow-wave essential for peristaltic movement of food and fluid. Furthermore, this model can be used to examine novel hypotheses concerning aberrant mechanisms that may underlie different motility disorders.

Development of Gastric Motility Network (GMN) model
In the present study, we developed a computational model consisting of a longitudinal arrangement of biophysically based ICC and SM cells found in the stomach wall. Since ICCs and their network are regarded as the key players for generation of pacemaker potentials and propagation of slow waves, we investigated whether and how the intrinsic properties of ICCs that are also modulated by enteric neural inputs enable inter-ICC coordination essential for entrainment. We also examined the crucial intercellular pathways that influence the entrainment range, important for long distance GSW propagation necessary for normal peristalsis. For this, we considered multiple interacting feedback pathways involving the intracellular key variables such as membrane potential, Ca 2+ and IP 3 concentration of an ICC as well as intercellular electrical coupling and exchange of second messengers. Further based on corroborative evidence that enteric neural inputs to ICCs show a rostrocaudal gradient [52,53], and indirect evidence that indicates that neurotransmitters enhance IP 3 in a variety of pacemaker cells [27][28][29], we assumed a gradient in the enteric neural stimulus along the ICC chain. Such a gradient was essential for the rostrocaudal entrainment of the slow-wave as also demonstrated in previous computational models [7,19].

Control of slow-wave characteristics by inter-ICC IP 3 exchange
First, our simulation results for single ICC and SM cells agree qualitatively and semi-quantitatively with the experimentally measured results (Fig 1), thus validating the GMN model components. Next, our results show that the GMN with its multi-stage feedback between IP 3 -Ca 2 + -V m combined with the intercellular pathways (Fig 2) can generate slow-wave propagation with a uniform frequency and a uniform rostrocaudal lag, i.e., entrainment (Fig 3). Particularly, the presence of inter-ICC IP 3 exchange in addition to the gap junction electrical conductance can increase the range of longitudinal entrainment (Fig 4). Interestingly, when IP 3 and Ca 2+ exchange is present between ICCs, the spacing of increasing Relative Lag is nonlinear, attributable to the diffusive nonlinearities of these molecules. On the other hand, the increasing spacing of Relative Lag is uniform with the resistive current exchange when electrical coupling alone is present (See Fig 4B1). Importantly, the presence of inter-ICC IP 3 exchange in addition to the electrical coupling minimizes the potential signs of bradygastric decoupled slow waves. Our observations are in accordance with the suggestion that bradygastria can potentially result from a failure of normal entrainment [60]. While the distinct intercellular mechanisms modeled in our network (exchange of second messengers and electrical coupling) are experimentally inseparable, the biophysical formulation included in the GMN model allowed examination of their individual contributions to longitudinal entrainment of ICCs. Second, we systematically evaluated the network behavior for increasing values of electrical conductance and IP 3 permeability (since Ca 2+ permeability was shown to have little to no impact on GSW entrainment). Surprisingly, the IP 3 exchange affected the features of the slowwave in a manner that can stabilize the propagating wave (from tachygastric decoupled slow waves for low IP 3 to a normal coupled gastric slow-wave for elevated IP 3 ). In particular, increased IP 3 permeability of gap junctions resulted in an increase in the SM Cell Period (reflective of bradygastria), whereas increased gap junction electrical conductance caused a decrease in the SM Cell Period (reflective of tachygastria). Thus, in combination with the contrasting effects of electrical conductivity, the IP 3 coupling could flexibly modulate the SM cell frequencies. We suggest that the apparent compensatory effect of IP 3 permeability could act as a brake on preventing runaway tachygastria (elevated SM cell frequencies, which is inversely related to SM Cell Period or pacemaking potential duration) in the event of increasing electrical coupling strengths. Increased gap junction density, which is reflected by increased gap junction coupling strength in our model, has been demonstrated in colonocytes during bacterial infections resulting in diarrhea generation [61]. Whether the same holds true for gastric motility disorders requires further experimental investigation. Even if the electrical gap junction density becomes exceedingly high, the inter-ICC IP 3 coupling makes sure that the SM cell frequency, and therefore, gastric contraction frequency stays in a physiologically plausible range.
Interestingly, increasing values of electrical conductance or IP 3 permeability led to an increased slow-wave velocity. In our results, the increase of slow-wave velocity was demonstrated by a decrease in the Total Lag of the slow-wave (Fig 5). These observations are very similar to the observations in studies of the dog small intestine where a negative gradient of gap junctions exists across the duodenum, jejunum and ileum and so does the slow-wave velocity in these compartments of the small intestine [62]. Electrophysiological experiments with gap junction enhancers (such as in Rotigaptide [63]) also support the observations of our computational study. As per our knowledge, such gap junction enhancers have not been employed as medicinal interventions yet, although they have the potential to elevate the reduced slow-wave velocity as observed in gastroparesis patients [60]. We would have to be very cautious though in introducing such enhancers, since slow-wave velocity has been reported to be already highly elevated in the gastric corpus of aged patients having gastroparesis with impaired peristalsis [64] and somewhat elevated in glucagon induced hyperglycemic dogs [6]. Future experiments could investigate the proper dosage of the gap junction enhancers so that the slow-wave velocity remains in a moderate range. For now, we can only speculate that in an intact stomach, the electrical gap junction coupling strength probably does not reach an exceedingly high value, and even if it has a moderately high value, Fig 5D1 and 5H1 show that the slow-wave velocity reaches a limit, which is within a physiologically reasonable range.

Distinct period-velocity relationships exist for increasing G ICC−ICC and P IP3
According to coupled oscillator theory, the direction of the slow-wave would depend on the frequency gradient of component pacemaker ICCs, whereas the slow-wave velocity would depend on the intercellular coupling strength [65]. Therefore, ideally, there should be a dispersion relationship between slow-wave velocity and frequency (expressed by SM Cell Period, in our case). However, previous studies of the canine gastric antrum [65] and porcine gastric corpus and antrum [66] suggest a dependence of slow-wave velocity on the observed SM Cell Period. The period-velocity relationship as noted in Fig 5H1 and 5H2 (higher the period, higher the velocity) for increasing permeability of IP 3 supports these in vitro [65] and in vivo [66] findings. In contrast, the period-velocity relationship noted in Fig 5D1 and 5D2 (lower the period, higher the velocity) for increasing electrical coupling strength is somewhat surprising. Because, if the period decreases, the velocity of a slow-wave should drop as a response to encroachment on the tail of the previous slow-wave (see Fig 7 in [66] for a more intuitive understanding). Although such a relationship has been observed along the length of the intestine [62] under normal conditions, its occurrence is probably due to the emergence of frequency plateaus, which in fact reflect localized decoupling due to the reduced gap junction density at certain points along the intestinal length [59]. Frequency plateaus are not observed in the stomach under normal conditions. Consequently, localized decoupling due to the reduced gap junction density along the length of the stomach is an unlikely explanation for the observed period-velocity relations in Fig 5D1 and 5D2. We believe that the less pronounced effect of increasing G ICC−ICC on velocity (reflected as Total Lag in Fig 5D1) and period (illustrated in Fig 5D2) compared to the effect of increasing P IP3 on the same quantities (Fig 5H1  and 5H2) makes sure that the trajectory of a slow-wave having increased G ICC−ICC does not appear synchronously with the tail of the previous slow-wave, thus avoiding the positive correlation between period and velocity. In our model results, we are concerned only with the velocity from the rostral end to the caudal end. Measurement of velocity gradient (if there is any) and a non-uniform modeling of gap junction conductance [15] along with the measurement of gap junction density in various animal models could shed further light on the model predictions.

Exchange of second messengers ensures robust rostrocaudal entrainment when biological variability is present
Finally, we showed that while variability in electrical coupling results in loss of rostrocaudal entrainment (with changes in cell periods reflective of bradygastria and/or tachygastria), presence of IP 3 exchange despite variability in its permeability, can restore rostrocaudal entrainment ( Fig 6) and therefore, is essential for peristalsis. Here, inter-ICC IP 3 exchange offers resilience to the increased variability in coupling strength. This is noteworthy because gap junctions in ICC networks can have variable density in different compartments of the stomach [15,16] and the observed consequences of increasing the variability in inter-ICC coupling in our simulations highlight that this might be a feature important for regulating the wave velocity while maintaining the range of entrainment. Our model with the distinct intercellular mechanisms (exchange of second messengers and electrical coupling) in combination with the intracellular feedback pathways and rostrocaudal neural innervation gradient allows SM cells to oscillate with a moderate range of frequencies and the gastric slow-wave to propagate with a broad range of velocities (Fig 5). Thus, the GMN confers high robustness to the rostrocaudal entrainment of ICCs, including under instances of biological variability in coupling pathways (Fig 6). Robustness is a fundamental property of evolvable complex biological systems [67]; a simple mechanism cannot handle extreme changes in physical quantities. Our model offers a new integrative framework for conceptualizing GSW propagation and regulation as a robust system of dynamically coupled oscillators fed with second messengers through intercellular exchange.

Biological and theoretical assumptions of the GMN model
Although numerous biological mechanisms are involved in the orchestration of gastric motility [68], it is well-established that the peristaltic movement of food/liquid is mediated by a propagating wave of smooth muscle contractions along the stomach wall [5,13,14,69]. This socalled gastric slow-wave involves electrical activity transmitted aborally within the SM cells [55,70]. However, the SM cells on their own cannot produce such a regenerative wave [70,71]. The regenerative electrical activity responsible for the GSW is known to originate largely in the intrinsic pacemaker ICCs [5,44]. Different sets of ICCs innervate circular and longitudinal muscle cells (e.g., myenteric ICCs, intramuscular ICCs [5,13]) and all ICCs may not contribute similarly to pacemaker potential generation and GSW propagation [68]. Our ICC cell model is derived primarily from experimental work on myenteric ICCs, which are widely accepted as intrinsic pacemakers involved in the generation of the pacemaker potential required for GSW propagation [5]. We incorporated known membrane properties in the ICCs and SM cell models using the conductance-based formalism and hand-tuned the conductance parameters to match experimental voltage recordings in these cells (Fig 1C and 1D) [5]. Although our model can reproduce rather accurately the established properties of a gastric slow-wave, it should be noted that the description of the underlying mechanisms is by no means exhaustive. For example, we have chosen to use a Hodgkin-Huxley formalism for the ion currents whereas Markovian models would allow for a more detailed description of the complex kinetics of processes (activation, deactivation, inactivation, and recovery from inactivation) that the channels exhibit. However, it can be very challenging to meet the information requirements for defining the transitional rate constants of a Markovian model, in addition to higher computational processing demands [72]. Our assumptions for model construction are based on various experimental models and assays (in vitro whole stomach, cell culture, etc.) as no one experiment or model system is available that encompasses all the biological details laid out in the Methods section. Despite numerous parameters and nonlinear equations in the model which capture the biological complexities at molecular, cellular and network scales, we note that the simulated GSW characteristics closely matched the generally conserved features in the mammalian stomach as shown in S2 Table [5,39]. Due to the nonlinear system of equations describing the model, we also simulated the default network for up to 4000s. Entrainment persisted and no abrupt transitions were observed. However, our simulations only provide an estimate of a theoretical steady-state of the network and we cannot rule out the possibility that very slow transients could affect long-term behavior [73]. Under natural circumstances, several factors not considered here could also influence the network activity.
Although ICCs are intrinsic pacemakers, the GSW can significantly be impacted by the coupling between them [8,74,75]. Gap junctions are well-established pathways for such coupling [76]. The gap junctions carry various ions and molecules. Some of these (e.g., K + ) are passively conducted due to the voltage gradients between adjacent cells (the electrical component), while others spread when there are periodic depolarizations in ICCs that generate these molecules in abundance. The latter are typically second messengers such as IP 3 and Ca 2+ , which involve chemical coupling via connexin channels [22][23][24]35,36]. The diameter of gap junctional pores allows a wide enough path for Ca 2+ and IP 3 to move across the cells through gap junction connexins [77]. To model electrical coupling, a simple ohmic conductance has been widely assumed by previous modeling studies [7,8,19,69,78,79], supported by findings that pacemaker current generated in ICCs is transmitted to the SM cells by gap junction channels located between the ICCs and SM cells [15,17]. We additionally assumed a distinct exchange of second messenger (Ca 2+ and IP 3 ) between adjacent ICCs. Second messenger exchange is modeled in such a way that it enables Ca 2+ induced Ca 2+ release and IP 3 production. The latter may involve a positive feedback mechanism as with many biomolecule syntheses (enzymes and others). Experimental evidence suggests that diacylglycerol, IP 3 , or one of their metabolites can activate phospholipase C (PLC). The activation of PLC facilitates subsequent production of IP 3 via an enzymatic pathway (see Eq 1.2 in [80]). Although Ca 2+ induced Ca 2+ release is a widely assumed phenomenon within a biological cell, IP 3 production due to positive feedback has also been suggested as a plausible event [81]. Ca 2+ permeability was set an order of magnitude lower than that of IP 3 permeability. This assumption was based on the fact that the range of diffusion of IP 3 is orders of magnitude higher compared to free Ca 2+ (24 μm compared to 0.1 μm). IP 3 also degrades much more slowly than free Ca 2+ (1 sec vs. 0.00003 sec) [82]. Previous modeling studies of intercellular Ca 2+ waves in astrocytes [35] and smooth muscle cells [25] support this approach of modeling Ca 2+ and IP 3 permeability. The exchange of IP 3 was modeled as proportional to the IP 3 concentration gradient between adjacent cells (see Methods, also [83]). Although IP 3 can diffuse within the cell cytoplasm and not necessarily through the gap junctions, we assumed that IP 3 only moves to the adjacent cells based on the concentration gradient. Presence of electrical coupling alone between ICCs misses a key mechanism of second messenger exchange which appears essential for robust entrainment of ICCs (Figs 5 and 6). Although Ca 2+ exchange via gap junctions does not play a significant role in ICC entrainment (S3 Fig), the range of ICCs entrained, that is limited by electrical coupling alone, can be extended by inter-ICC exchange of IP 3 (Fig 4). For Figs 3-5, we considered constant values of gap junction coupling parameters and for Fig 6, we assumed uniform distribution of gap junction coupling strength, which resembles the gap junction density. Because of the scarcity of data regarding gap junction distribution, we limited our simulations to only uniform distribution.
In a chain of intrinsic pacemakers, oscillators have their own distinct intrinsic frequency. In an arbitrarily long chain of oscillators, entrainment can emerge provided there exists a linear frequency gradient with fixed frequency difference between the ends [84]. Interestingly in the stomach (also throughout the GI tract), it is known that there is a rostrocaudal frequency gradient in the ICC pacemaker cells. Although such gradients may be achieved by numerous intrinsic and network mechanisms [7,19], we assume that such a gradient in intrinsic frequencies can be achieved due to a maintained gradient in [IP 3 ] production rate. Presently it is unclear what sets this gradient tone in IP 3 production rate. One possibility could be IP 3 produced by activation of muscarinic receptors due to graded distribution of cholinergic neuron inputs in different functional compartments of the GI tract [52,53]. Intrinsic frequency differences could be a consequence of intrinsic oscillatory properties, inter-oscillatory coupling strengths, and extrinsic inputs such as from neural innervation [30,85]. Although peristalsis can occur even without the help of neural excitation [68], under most circumstances, the enteric nervous system, provides excitation of the musculature required for the stomach wall contraction [26,86], primarily via ICCs. Enteric neurons located in the Myenteric Plexus preferentially directly influence ICCs through neurotransmitters released at synapses which then connect to muscle cells via gap junctions. Since ICCs are closer to the nerve terminal endings and have the muscarinic receptors (M2 and M3) that are responsive to the neurotransmitters, the effect of the Myenteric neurons on ICCs is more important (much smaller gap for neurotransmitters to diffuse) than the direct link to the SM cells (far away and with lesser innervation). Acetylcholine, the primary neurotransmitter released from enteric neurons, is broken down by Acetylcholine esterase, thus preventing it from reaching receptors on SM cells. Consequently, we considered the enteric neural stimulus (β) to act primarily on ICCs in our GMN model. The animal models that lack ICCs with muscarinic receptors show that there is little or no cholinergic response in SM cells because Acetylcholine is broken down before it can be taken up by receptors on SM cells [26]. Although excitatory and inhibitory inputs modulate the enteric network including the ICCs and SM cells [26], here we simplified the effect of neural input to the network to single parameters primarily leading to depolarization of ICC via IP 3 production rate, beta (β) and SM cells via NSCC posited due to muscarinic receptor activation (rlig). More detailed representation can enhance the scope of our model in future work. In our model, we maintained a linear gradient of β to ensure that rostral ICCs successively entrained caudal ICCs. This way, we assembled a chain of realistic ICCs which show a rostrocaudal gradient in their intrinsic frequencies like that observed in the stomach. In the heart, fast-pacing sinoatrial node cells entrain the slow-pacing atrioventricular node cells [87], like what is observed in an intact stomach.
In addition, we assumed that voltage-dependent IP 3 synthesis occurs in pacemaker ICCs in our model. Voltage-dependent IP 3 synthesis has been demonstrated in canine coronary artery tissue where a vasoconstrictor (thromboxane A2 analogue U 46619) was used to stimulate IP 3 production [88]. Hyperpolarization of the tissue resulted in reduced IP 3 concentration, which was then restored by depolarization suggesting a direct coupling between membrane voltage and IP 3 concentration. Depolarization induced Ca 2+ influx through voltage-gated Ca 2+ channels can affect enzymatic pathways that in turn enhance IP 3 production [89,90]. Similarly, in guinea pig coronary myocytes [91] and in the circular smooth muscle of the guinea pig gastric pylorus [20], voltage-dependent IP 3 production has been observed. Based on the above information, we followed the formulation in [42] and included a corresponding term in Eq 4.
Besides being intrinsic pacemakers, ICCs are also involved in neurotransmission, setting the membrane potential of SM cells, and in stretch sensing [51]. Future experiments examining effects of neural input on ICCs may shed further light on the neural contribution to pacemaker potential generation and slow-wave propagation. Empirical measurements of the neural innervation on ICCs at different points along the stomach would strengthen the validity of our model, offering further hypotheses for the mechanisms underlying genesis of functional and pathological slow waves in the stomach. The mathematical modeling framework outlined in our study is a step in this direction and provides an exploration testbed for precise modulation of intra/intercellular pathways to examine their role in the above-mentioned phenomena.

Methods
The model network of ICC-SM cells is described in Fig 2. The network consists of 42 ICCs and 42 SM cells implemented and simulated in MATLAB 2019a (Mathworks.inc). Together, the extensive biological details of ICC and SM cell physiology and the assumed intercellular coupling resulted in a large network model with 1554 nonlinear differential equations (23 for each ICC and 14 for each SM cell). The parameter values were chosen from published models [7,43,44].

ICC cell model
We adopted a well-described conductance-based model by Corrias et al., 2008 [44] for the ICCs (also see [7,19]). The rate of change in the membrane potential, V ICC , for each ICC is as follows: where V ICC is the ICC membrane potential, C m,ICC is the ICC membrane capacitance, I ion,ICC represents the summation of all the ionic currents in ICC, and I g,ICC represents the current through the gap junctions between the ICCs. The different ionic currents included in the model are listed in Table 1. The dynamics of the voltage-dependent gating variables for the ionic currents follow the well-known Hodgkin-Huxley formalism, where each current is represented by a battery (electrochemical driving force) in series with a variable resistance and the cell membrane as a capacitor in parallel. The rate equations for the nonlinear voltage-dependencies and kinetics of the membrane ionic currents are provided in the S1 Text. The combined actions of these ion channel currents reproduce the three phases of a gastric action potential: depolarization, plateau phase, and repolarization. Detailed descriptions of these ionic currents are provided in the github repository (https://github.com/ashfaq-polit/Slow_ waves_in_the_stomach) and are based on previous modeling studies [43,44,92].

Ca 2+ -IP 3 dynamics in ICCs
The equation guiding the intracellular Ca 2+ dynamics is modeled similar to Fall and Keizer [93] as follows: where I L−type and I VDDR are described in Table 1, J leak is a leakage flux between the pacemaker region and the cytosol, J PMCA is the Ca 2+ flux through the plasmalemmal Ca 2+ pump, F is the Faraday constant, f c is a dimensionless constant, V cyto represents the cytosolic volume fraction within the ICC, and J g;Ca 2þ is the Ca 2+ flux due to inter-ICC coupling (see Eq 6). Here, J g;Ca 2þ is an additional term compared to Eq 2 in [7]. The Ca 2+ flux in a submembrane space triggered by IP 3 -operated stores in the endoplasmic reticulum (ER) is important for initiating a gastric action potential. The Ca 2+ concentration in the submembrane space is modeled by the following equation: where J NaCa and J Uni represent mitochondrial Ca 2+ fluxes, J ERout and J SERCA represent outward Ca 2+ fluxes from the ER and inward Ca 2+ fluxes through the ER. V mito , V ER , and V SS represent the volume fractions for the mitochondria, endoplasmic reticulum, and pacemaker submembrane space, respectively. Increases in intracellular IP 3 concentration in each ICC are assumed to depend on: 1) voltage-dependent IP 3 increase, 2) inter-ICC coupling-dependent IP 3 increase, and 3) neurotransmitter-induced IP 3 increase, whereas linear and non-linear degradation of IP 3 decrease its concentration. Based on these assumptions, the rate of change in IP 3 concentration is given as follows:

Intercellular coupling between ICCs
Adjacent ICCs are assumed to be coupled via gap junctions with electrical conductance for passive ionic movement between cells [17]. The inter-ICC gap junction current is: where i denotes the index for the source ICC and j denotes the index for adjacent ICCs; G represents electrical conductance of gap junctions.

Second messenger exchange between ICCs
An exchange of second messengers, namely, Ca 2+ and IP 3 can occur via gap junctions [22,23,25,83]. This phenomenon has not been considered in previous computational models of ICCs. The flux describing Ca 2+ exchange between adjacent ICCs is modeled by a term complementing Eq 2. Here P Ca 2þ denotes the permeability coefficient for Ca 2+ and [Ca 2+ ] x are the intracellular concentrations of Ca 2+ in the corresponding cells where x = i/j. The P Ca 2þ value has been taken from a theoretical study done in hepatocytes [94]. The inter-ICC IP 3 flux is modeled using a similar formalism: Here P IP3 denotes the permeability coefficient for IP 3 and [IP 3 ] x are the intracellular concentrations of IP 3 in the corresponding cells where x = i/j. An experimentally determined value of P IP3 is rare if not non-existent due to the technical difficulties of measuring [IP 3 ] in tissue preparations. Hence, the permeability values used in model simulations were adjusted similar to [25,35].

SM cell model
The SM cell model was adopted from Corrias et al., 2007 [43]. Like ICCs, the conductancebased rate of change in membrane potential is given as follows: C m;SM dV SM dt ¼ À P I ion;SM þ I g;SM À I coup � � ; ð8Þ where V SM is the SM cell membrane potential, C m,SM is the SM cell membrane capacitance, I ion,SM represents the summation of all the ionic currents in the SM cell, I g,SM represents the current through the gap junctions between the SM cells, and I coup represents the coupling current from ICC to its corresponding SM cell. In Eq 10, G SM represents the electrical conductance of gap junctions between adjacent SM cells, i and j represent the same entities as described in Eq 5. Although a bidirectional movement of ions can occur via gap junctions, here we assume that during entrainment depolarization spreads unidirectionally from ICC to SM and that V ICC is always more positive than V SM [70]; hence, the coupling current, I coup is present only in Eq 8. The various ionic currents necessary to generate the SM action potential are listed in Table 1.
Although the enteric neural stimulation primarily modulate ICC activity, the muscarinic neurotransmitters can modulate the non-selective cation current in SM cells [95]. The NSCC current in SM cells is modeled by the following equations, similar to [43]:

Model simulation and analysis
The ICC-SM model network consisting of 42 ICCs and 42 SM cells was implemented and simulated in MATLAB 2019a (Mathworks.inc) using built-in function ODE15s and variable stepsize. The total runtime for each simulation was 900 seconds (each simulation lasted approximately 8 hours on an Intel Xeon (R) CPU, 32 GB RAM Desktop computer). The outcome of the network was interpreted from the spatiotemporal membrane potential, Relative Lag and period of SM cells (explained in Results section). Statistical tests (t-test, one-way ANOVA) were performed as required and an α cutoff of 0.05 was chosen for statistical significance.   simulation (Fig 5H1). The asterisk represents the P IP3 value used in the default network in