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–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²⁺ buffering mechanism, a membrane potential dependent intracellular concentration change in IP₃, and IP₃ receptor mediated Ca²⁺ 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²⁺ exchange to moderate the intracellular Ca²⁺ 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 (left, ICC pacemaker potential; right, SM cell slow-wave potential). 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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Biophysical models for ICCs and SM cells involved in gastric motility.

Schematics showing (A) anatomical divisions of the stomach and the arrangement of various cell types in the stomach wall, (B) membrane ionic currents and intracellular Ca²⁺-IP₃ components included in an interstitial cell of Cajal (ICC) model, and (C) membrane ionic currents in a smooth muscle (SM) cell model. See Table 1 for symbols and details of ionic currents. [Ca]_i = Intracellular Ca²⁺, ER = endoplasmic reticulum, Na⁺/Ca²⁺ = Sodium/Calcium exchange pump, SR = Sarcoplasmic reticulum, SS = Submembrane space, (D) Simulated rhythmic membrane potential dynamics in an ICC (left panel) and an SM cell (right panel), respectively.

https://doi.org/10.1371/journal.pcbi.1009644.g001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Gastric Motility Network model architecture.

The GMN is constituted of a chain of nearest-neighbor coupled interstitial cells of Cajal (ICC) and associated smooth muscle (SM) cells. The ICCs have electrical and second messenger (Ca²⁺ and IP₃) based coupling, while the ICC to SM and SM to SM couplings are only electrical. There is a negative gradient in enteric neural innervation (β) that modulates IP₃ 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 and text for further details.

https://doi.org/10.1371/journal.pcbi.1009644.g002

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 nearest-neighbor 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–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²⁺ and IP₃ [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²⁺ and IP₃ concentrations. The intracellular Ca²⁺ is further divided into two compartments: Ca²⁺ in the submembrane space (SS) spanning a narrow region between the ER and mitochondria ([Ca]_SS,i) and near membrane Ca²⁺ ([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₃ concentration (see V_m−IP₃ feedback in Fig 2, Eq 4 in Methods, also [20,27,42]); The intracellular IP₃ concentration, [IP₃]_i, in turn affects the [Ca]_SS,i concentration through the outward Ca²⁺ flux from the ER (J_Erout) [46–48], as highlighted by the IP₃−Ca²⁺ 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²⁺−V_m feedback) [49,50].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Ionic currents in the ICC model.

https://doi.org/10.1371/journal.pcbi.1009644.t001

An initiating event for pacemaker activity (enteric neural innervation) is assumed to increase the IP₃ production rate and is modeled as a constant rate, β, for each ICC (see Eq 4 in Methods). Biologically, increases in IP₃ levels can originate from endogenous release of 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₃] results in an increase in the [Ca]_SS,i mimicking the endogenous release of Ca²⁺ from IP₃-operated stores in the ER. In response to [Ca]_SS,i increase, the Ca²⁺ uniporter on the mitochondrial membrane is gated open, and Ca²⁺ 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²⁺ ions from the submembrane space than had previously entered from the ER, causing a temporary drop in the subspace Ca²⁺ 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²⁺ 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₃ 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₃ and Ca²⁺ 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 S1 Fig. As we can observe from S1 Fig, the rise in [Ca²⁺]_i results in the depolarization of ICC and a subsequent upstroke of ICC membrane potential. The abrupt change in [Ca²⁺]_SS is due to the sharp rise of Ca²⁺ in the submembrane space because of a Ca²⁺ efflux from ER and the subsequent sharp drop because of Ca²⁺ movement into the mitochondria. Compared to the [Ca²⁺]_SS change, the change in intracellular IP₃ 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 3 highlights the entrainment and the resulting slow-wave in the GMN. In Fig 3A, a spatiotemporal map of the membrane potentials of the 42 SM cells demonstrates the rostrocaudal propagation of the slow-wave in the network. The grey shaded regions indicate the up-swing, and the black regions the down-swing in the membrane voltage of the SM cells. A single GSW cycle consists of successive phases of activity from the rostral to the caudal-most oscillatory units. In Fig 3B, we highlight the membrane potential of every 7^th SM cell in the network (the minimum value of the membrane potential during the last 7 cycles of the simulation for all 42 SM cells of the network is -63.43 (± s.d. = 0.41) mV, which is very close to the value reported in [56]). We define and measure SM Cell Period as the time between two consecutive peaks in the membrane voltage of an SM cell. An increase or decrease in the SM Cell Period reflects bradygastria and tachygastria, respectively. Relative Lag is defined as the time difference between the peak membrane voltage of SM₁ and SM_i, where i = 2,3,4…42, and Total Lag as the time between the peak membrane voltages of SM₁ and SM₄₂. Slow-wave velocity is inversely proportional to the Total Lag, provided that the Total Lag value reaches a constant value. At steady-state, Relative Lags should reach constant values along the length of the network for a normal slow-wave exhibiting entrainment. Any deviation of Relative Lag from constancy reflects decoupling.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Slow-wave propagation and entrainment in the GMN.

(A) Spatiotemporal plot of the SM cell membrane potential along the length of the stomach (vertical axis) and time (horizontal axis) of all 42 SM cells of the network. The direction of slow-wave propagation occurs from the rostral end of the network (representing the mid-corpus) to the caudal end (representing the terminal antrum). (B) The membrane potential of 7 equidistant SM cells in the 42-cell network. (C) Total Lag for the first 7 cycles of simulation and the last 7 cycles of simulation. (D) Relative Lag (RL) between the 1^st SM cell and i^th SM cell in the network for the last 20 cycles, where i = 2, 3, 4,…, 42. (E) Spatiotemporal map of periods of all 42 SM cells in the network for the last 20 cycles. The inset plots for the Relative Lag (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

In our model, the average SM Cell Period over the last 7 slow-wave cycles of the simulation for the boundary cells SM Cell 1 and SM Cell 42 were 17.7 (± s.d. = 0.01) seconds and 17.7 (± s.d. = 0.01) seconds respectively. These periods match closely with those observed in in vitro studies of guinea pig stomach SM cells (see Fig 1C in [14]). The computed Total Lag reached an approximately constant value of 20.8 (± s.d. 0.02) seconds resulting in a slow-wave velocity of ~0.29 cm/sec (Fig 3C), and it was consistent between consecutive cycles. The Relative Lag for each cell reached an approximately constant value (Fig 3D) as well as the periods for all the SM cells in the network (observed from the spatiotemporal map of SM Cell Periods in Fig 3E). These results indicate the entrainment of all the oscillatory units to a uniform gastric slow-wave 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₃ and Ca²⁺ 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 coupling between ICCs through electrical gap junctions alone (Fig 4A1), with coupling through second messenger exchange of Ca²⁺ and IP₃ (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²⁺ and IP₃ 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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Enhanced longitudinal entrainment. Synergy of electrical coupling and second messenger exchange preserve longitudinal entrainment along the entire length of the stomach.

(A1, A2 and A3) Structure of the network, when only electrical coupling is present (A1), when only exchange of second messengers is present (A2), and when both electrical coupling and exchange of second messengers are present (A3). (B) The Relative Lags for the three different cases shown in Fig 4A are measured from the last 20 cycles of respective simulations. (C) Spatiotemporal maps of SM Cell Periods in the three different cases respectively.

https://doi.org/10.1371/journal.pcbi.1009644.g004

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). For the lowest G_ICC−ICC value considered, the Relative Lag diagram illustrated in Fig 5B2, left shows decoupled slow waves. Here, cell #1–32 are entrained though with a period of ~17.7 sec. Decoupled slow waves are observed from cell #33–42. The SM Cell Period diagram for the same network depicted in Fig 5B3, left demonstrates an increase in SM Cell Period from cell #33–42. The increased cell period could be an underlying sign of bradygastria. Rostrocaudal entrainment was achieved only for higher values of the G_ICC−ICC as shown by the levelling out of Total Lag in the color-coded insets in Fig 5C. Overall, an increase in the strength of electrical coupling conductance produced an enhancement in the entrainment range.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Electrical gap junction and second messenger coupling strengths. The inter-ICC coupling strengths impact GMN entrainment, gastric slow-wave velocity and pacemaker frequency.

(A, E) G_ICC−ICC and P_IP3 are altered for the simulations in Fig 5B–5D and Fig 5F–5H, respectively. (B) Spatiotemporal map of membrane potential (B1), Relative Lags (B2), and spatiotemporal map of SM Cell Periods (B3) for the network when G_ICC−ICC = 0.35 nS (left panels, partially entrained) and 1.4 nS (right panels, entrained). (C) The Total Lag for changes in G_ICC−ICC is shown for the last 7 cycles of 900-sec simulations for different conductance values (nS) indicated by the color legend. The two distinct classes of responses are enlarged in the right panels (lower ones entrained, upper ones partially entrained). An increase in Total Lag indicates decrease in slow-wave velocity. (F) Spatiotemporal maps of membrane potential (F1), Relative Lags (F2), and spatiotemporal maps of SM Cell Periods (F3) for the network when P_IP3 = 4.0 sec^-1 (left panels, partially entrained) and 16.0 sec^-1 (right panels, entrained). The spatiotemporal membrane potential diagrams are shown for the last 120 seconds of the respective simulation, whereas the Relative Lag and the SM Cell Period diagrams are shown for the last 20 cycles of the respective simulation. (G) The Total Lag for changes in P_IP3 is shown for the last 7 cycles of 900-sec simulations for different conductance values (nS) indicated by the color legend. (D, H) For several networks, the mean Total Lag (odd numbered panels) and the SM Cell Period (even numbered panels) of the last 7 cycles for each network with respect to its G_ICC−ICC and P_IP3 can be fit by individual exponential function, respectively. For increasing values of G_ICC−ICC, an exponential fit (, where a₁ = 8.65, β₁ = 0.75, a₂ = 18.23, β₂ = 17.87) has been drawn along the mean value of Total Lag and for increasing values of P_IP3, another exponential fit (, where a₁ = 7.44, β₁ = 1.70e5, a₂ = 32.07, β₂ = 9.37) has been drawn along the mean values of Total Lag. Partially entrained networks have not been considered for equation fitting. For SM Cell Period calculation, the last cell (42^nd cell) has been considered as the representative cell of the network. An increase in Total Lag indicates decrease in slow-wave velocity while a decrease in SM Cell Period indicates increase in pacemaking frequency.

https://doi.org/10.1371/journal.pcbi.1009644.g005

Fig 5E–5G present the network outcome for increasing IP₃ permeability (values of P_IP3 two times above and below a default value of 8.0 sec^-1) at the default value for the electrical coupling conductance. The spatiotemporal maps of SM cell membrane potential and the corresponding Relative Lag and action potential periods of SM cells are shown for the lowest (Fig 5F1-3; left) and highest (Fig 5F1-3; right) P_IP3 values considered. Similar to increasing G_ICC−ICC values, the networks with increasing P_IP3 values demonstrated a transition from partial entrainment to complete rostrocaudal entrainment. For the lowest P_IP3 value considered, the network generates decoupled slow waves as shown in the Relative Lag diagram of Fig 5F2, left. 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₃ 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₃ 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²⁺ permeability () to test its effects on network entrainment. We noted that even two orders of magnitude changes above and below the default value of 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²⁺ compared to IP₃ for gap junctions [18]. Due to this lack of effect of on entrainment and gastric slow-wave characteristics, in what follows, we will only focus on changing P_IP3.

Inter-ICC IP₃ 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₃ 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₃ 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, in particular, provide evidence for the decoupled slow waves. Cell #1–28 are entrained in Fig 6E1, cell #1–27 are entrained in Fig 6E2 and there exists no visible entrained region in Fig 6E3. The caudal ends in Fig 6E1 (cell #29–42) and Fig 6E2 (cell #28–42) and the entire network in Fig 6E3 are decoupled. The SM Cell Period diagrams in Fig 6F1–6F3 show a combination of reduction and increase in SM Cell Period. The entrained cells in Fig 6F1 and 6F2 exhibit a reduction in SM Cell Period (~16.4 sec compared to ~17.7 sec in 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₃ exchange, decoupled slow waves and potential signs of bradygastria and/or tachygastria can coexist.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Entrainment in GMN is resilient to variability in inter-ICC coupling strengths.

(A,G) The variability in only G_ICC−ICC (A) and both G_ICC−ICC and P_IP3 (G) (20%, 50%, and 100% variability) are shown around their mean values. (B, H) Total Lag of the network without inter-ICC IP₃ exchange for three different variabilities in G_ICC−ICC (B) and of the network with inter-ICC IP₃ exchange for three different variabilities in G_ICC−ICC and P_IP3 (H). The latency to entrainment cannot be inferred from any of the panels in (B), since the networks were only partially entrained. However, the latency to entrainment can be measured from the first two panels in (H), since the Total Lag has attained constant value in these two cases. (C, I) Total Lag is shown for each case for the last 7 cycles using violin plot. The asterisk (*) symbol represents a statistically significant difference between the corresponding quantities. Spatiotemporal diagrams of membrane potential (D, J), Relative Lags, (E, K) and spatiotemporal maps of SM Cell Periods (F, L) for the network are shown for three different variabilities in G_ICC−ICC for the network without inter-ICC IP₃ exchange and in G_ICC−ICC and P_IP3 for the network with inter-ICC IP₃ exchange, respectively. The spatiotemporal membrane potential diagrams are shown for the last 120 seconds of the respective simulation, whereas the Relative Lag and the SM Cell Period diagrams are shown for the last 20 cycles of the respective simulation.

https://doi.org/10.1371/journal.pcbi.1009644.g006

Next, we included the inter-ICC IP₃ 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₃ 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 that the value of Total Lag for the last 7 cycles of the respective simulations was 21.37 (± s.d. 0.01) seconds for 20% variability and 22.35 (± s.d. 0.1) seconds for 50% variability. 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₃ can preserve the longitudinal entrainment to a greater extent along the length of the stomach, thereby preventing decoupling and restoring the normal behavior.

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²⁺ and IP₃ 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₃ in a variety of pacemaker cells [27–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₃ 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₃-Ca²⁺-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₃ exchange in addition to the gap junction electrical conductance can increase the range of longitudinal entrainment (Fig 4). Interestingly, when IP₃ and Ca²⁺ 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₃ 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₃ permeability (since Ca²⁺ permeability was shown to have little to no impact on GSW entrainment). Surprisingly, the IP₃ exchange affected the features of the slow-wave in a manner that can stabilize the propagating wave (from tachygastric decoupled slow waves for low IP₃ to a normal coupled gastric slow-wave for elevated IP₃). In particular, increased IP₃ 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₃ coupling could flexibly modulate the SM cell frequencies. We suggest that the apparent compensatory effect of IP₃ 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₃ 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₃ 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₃ 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₃ exchange despite variability in its permeability, can restore rostrocaudal entrainment (Fig 6) and therefore, is essential for peristalsis. Here, inter-ICC IP₃ 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 so-called 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₃ and Ca²⁺, which involve chemical coupling via connexin channels [22–24,35,36]. The diameter of gap junctional pores allows a wide enough path for Ca²⁺ and IP₃ 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²⁺ and IP₃) between adjacent ICCs. Second messenger exchange is modeled in such a way that it enables Ca²⁺ induced Ca²⁺ release and IP₃ production. The latter may involve a positive feedback mechanism as with many biomolecule syntheses (enzymes and others). Experimental evidence suggests that diacylglycerol, IP₃, or one of their metabolites can activate phospholipase C (PLC). The activation of PLC facilitates subsequent production of IP₃ via an enzymatic pathway (see Eq 1.2 in [80]). Although Ca²⁺ induced Ca²⁺ release is a widely assumed phenomenon within a biological cell, IP₃ production due to positive feedback has also been suggested as a plausible event [81]. Ca²⁺ permeability was set an order of magnitude lower than that of IP₃ permeability. This assumption was based on the fact that the range of diffusion of IP₃ is orders of magnitude higher compared to free Ca²⁺ (24 μm compared to 0.1 μm). IP₃ also degrades much more slowly than free Ca²⁺ (1 sec vs. 0.00003 sec) [82]. Previous modeling studies of intercellular Ca²⁺ waves in astrocytes [35] and smooth muscle cells [25] support this approach of modeling Ca²⁺ and IP₃ permeability. The exchange of IP₃ was modeled as proportional to the IP₃ concentration gradient between adjacent cells (see Methods, also [83]). Although IP₃ can diffuse within the cell cytoplasm and not necessarily through the gap junctions, we assumed that IP₃ 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²⁺ 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₃ (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₃] production rate. Presently it is unclear what sets this gradient tone in IP₃ production rate. One possibility could be IP₃ 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₃ 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₃ synthesis occurs in pacemaker ICCs in our model. Voltage-dependent IP₃ synthesis has been demonstrated in canine coronary artery tissue where a vasoconstrictor (thromboxane A2 analogue U 46619) was used to stimulate IP₃ production [88]. Hyperpolarization of the tissue resulted in reduced IP₃ concentration, which was then restored by depolarization suggesting a direct coupling between membrane voltage and IP₃ concentration. Depolarization induced Ca²⁺ influx through voltage-gated Ca²⁺ channels can affect enzymatic pathways that in turn enhance IP₃ 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₃ 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.

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: (1) 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²⁺-IP₃ dynamics in ICCs

The equation guiding the intracellular Ca²⁺ dynamics is modeled similar to Fall and Keizer [93] as follows: (2) 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²⁺ flux through the plasmalemmal Ca²⁺ pump, F is the Faraday constant, f_c is a dimensionless constant, V_cyto represents the cytosolic volume fraction within the ICC, and is the Ca²⁺ flux due to inter-ICC coupling (see Eq 6). Here, is an additional term compared to Eq 2 in [7].

The Ca²⁺ flux in a submembrane space triggered by IP₃-operated stores in the endoplasmic reticulum (ER) is important for initiating a gastric action potential. The Ca²⁺ concentration in the submembrane space is modeled by the following equation: (3) where J_NaCa and J_Uni represent mitochondrial Ca²⁺ fluxes, J_ERout and J_SERCA represent outward Ca²⁺ fluxes from the ER and inward Ca²⁺ 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₃ concentration in each ICC are assumed to depend on: 1) voltage-dependent IP₃ increase, 2) inter-ICC coupling-dependent IP₃ increase, and 3) neurotransmitter-induced IP₃ increase, whereas linear and non-linear degradation of IP₃ decrease its concentration. Based on these assumptions, the rate of change in IP₃ concentration is given as follows: (4) where P_MV is maximal rate of voltage-dependent IP₃ synthesis, k_V is the half-saturation constant for voltage-dependent IP₃ synthesis, η is the rate constant for linear IP₃ degradation, V_m4 is maximal value for the nonlinear IP₃ degradation, k₄ is half-saturation constant for the nonlinear IP₃ degradation, β represents an enteric neural innervation that modulates IP₃ production, and J_g,IP3 is the IP₃ flux due to inter-ICC coupling (see Eq 7). Here, J_g,IP3 is an additional term compared to Eq 3 in [7].

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: (5) 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²⁺ and IP₃ 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²⁺ exchange between adjacent ICCs is modeled by a term (6) complementing Eq 2. Here denotes the permeability coefficient for Ca²⁺ and [Ca²⁺]_x are the intracellular concentrations of Ca²⁺ in the corresponding cells where x = i/j. The value has been taken from a theoretical study done in hepatocytes [94].

The inter-ICC IP₃ flux is modeled using a similar formalism: (7)

Here P_IP3 denotes the permeability coefficient for IP₃ and [IP₃]_x are the intracellular concentrations of IP₃ 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₃] 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 conductance-based rate of change in membrane potential is given as follows: (8) (9) (10) 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]: (11) (11.1) where G_NSCC,SM and m_NSCC,SM represent the conductance and voltage dependency of the activation kinetics of the NSCC channel in SM cells. r_lig is a constant representing the level of channel activity proportional to muscarinic stimulation; This value is set uniformly across all the SM cells in the model. The variable h_Ca describes the facilitation effect that intracellular Ca²⁺ has on I_NSCC,SM as formulated in Eq 11.1.

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 step-size. 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.