Neural regulation of slow waves and phasic contractions in the distal stomach: a mathematical model

Gastric digestion is facilitated by coordinated contractions of gastric smooth muscles, which are governed by a range of regulatory mechanisms. These include rhythmic bioelectrical activity, named gastric slow waves (SWs), as well as neural and humoral factors [1, 2]. At the cellular level, SWs are generated and propagated by a network of excitable pacemaker cells called interstitial cells of Cajal (ICC), which are situated both between and within layers of gastric smooth muscles [3]. In conjunction with innervation from enteric neurons, SWs generated by ICC govern the phasic contractions of gastric smooth muscle cells (SMCs). Populations of excitatory and inhibitory enteric motor neurons release neurotransmitters at both ICC and SMC to regulate gastrointestinal motility [4, 5]. We postulate that an enhanced physiological understanding of the neural regulation of gastric motility will improve neuromodulation therapy for gastric disease, as has occurred in neuromodulation-based treatment of cardiac arrhythmogenesis and overactive bladder among other diseases [6, 7].

For clarity, throughout this work, SW refers to rhythmic electrical events, while phasic contraction refers to the contractions of gastric smooth muscles associated with SWs. ICC intrinsically undergo rhythmic depolarisations giving rise to SWs, while also being electrically connected in a syncytium [4]. ICC can be excited by each other, such that the region with the highest intrinsic frequency, known as the pacemaker region, entrains SWs across the entire stomach, thereby maintaining a sustained, regular propagation pattern. In large monogastric mammals, SWs are observed to emerge in the corpus along the greater curvature and form circumferential bands propagating in an antegrade direction towards the pyloric sphincter [8–11]. In other monogastric species, such as rodents, the exact location of the pacemaker region differs, but distinct functional regions related to motility remain evident [12].

Neural regulation of gastric motility is important, and its roles vary across the stomach [13]. Several neurotransmitters have been identified; acetylcholine (ACh) and tachykinins (TK) have roles in excitatory neuromuscular transmission, while nitric oxide (NO) and purinergic neurotransmission have major roles in the inhibitory pathway [14]. Neurotransmitters influence the frequency and amplitude of SWs and the underlying tonic state of the stomach, thereby affecting gastric volume, mixing, and emptying. The stomach is primarily innervated by the vagus nerve, forming a part of the brain-gut axis [2]. Pre-enteric, vagal efferent neurons innervate enteric nervous system (ENS) neurons at ganglia in the myenteric plexus, between the longitudinal and circular smooth muscle layers, but do not directly innervate the muscle [14, 15]. ENS neurons directly innervate ICC, SMC, and other ENS neurons, while, in contrast, efferent vagal neurons innervate ENS neurons but not ICC or SMC [2].

While experimental work has identified a variety of candidate mechanisms enacting post-junctional neurotransduction, the relative importance of each mechanism in relation to gastric motility control is yet to be clarified. Complexity arises from the interaction and targeting of both ICC and SMC by enteric neurons, and the variety of expressed channels in both these cell types that are involved in post-junctional events. Existing models of gastrointestinal electrophysiology or motility focus on either ICC [16–19] or SMC [20–22]. However, none of these models have incorporated the effects of neural regulation on the properties of the phasic contractions which enable gastric motility. There exists scope to examine the influence that putative mechanisms of inhibitory neural regulation, termed 'effector components', have on SWs and phasic contractions in relation to experimental observations. Furthermore, a modelling study can predict the interaction between inhibitory and excitatory neural regulation. The model developed herein addresses this knowledge gap by incorporating the effects of efferent neural regulation on gastric ICC and gastric SMC. It models how cell level mechanisms interact to control SW activity. This improvement enables the computational study of how one portion of the neural regulation pathways to the stomach control gastric motility [1].

Inhibitory neural regulation via nitrergic and purinergic pathways targets both ICC and SMC in the distal stomach. The nitrergic pathway acts through the activation of cyclic guanosine monophosphate (cGMP) intracellular signalling pathways by NO. Potential targets in ICC include Ano1 [5], transient receptor potential melastatin (TRPM) [23–25], and large-conductance K⁺ (BK) channels [26, 27] though the effects at BK are better documented in intestinal SMC rather than gastric ICC [28, 29]. TRPM channels are non-specific cationic conductance channels (NSCCs). Purinergic neurotransduction via apamin-sensitive small-conductance K⁺ (SK) channels has also been identified in the distal stomach [30], likely via fibroblast-like (PDGFRα⁺) cells conducting to SMC [31, 32]. By observation of fast and slow inhibitory junction potentials (IJPs) it was deduced that ICC are not involved in the purinergic pathway [33]. Additionally, the intracellular Ca²⁺ sensitivity of the SMC contractile apparatus is affected by nitrergic stimulation [34].

Excitatory cholinergic neural regulation affects the frequency of SWs in ICC. Effects on the ICC SW arise from increased IP₃ production of due to muscarinic activation by acetylcholine (ACh). This causes increased Ca²⁺ release from intracellular stores, through the IP₃-activated Ca²⁺ channel (IP₃R), thereby resulting in a higher frequency of SW activation [35]. Experiments have shown that tachykinins, which may be co-transmitted with ACh, impart an additional excitatory effect on SW frequency at stimulation frequencies of 10 Hz or higher [36].

A mathematical model was developed to elucidate the cellular mechanisms for enteric neural regulation of gastric SWs and associated phasic contractions. Five key cellular mechanisms were incorporated for inhibitory and excitatory neural regulation of ICC and SMC, as well as their interactions and effects on gastric motility. Four of these were inhibitory regulation mechanisms via NO modulation of Ano1, NSCC, and SK channels, and inhibitory regulation of SMC calcium sensitivity. The fifth was a mechanism of the excitatory regulation of IP₃. The model was parameterised to reflect the properties of the rat distal stomach, for alignment with literature on rat gastric neural regulation [36, 37]. Our model allows gastric motility to be simulated under varied enteric motor neuron stimulation frequency and assess the contribution of different post-junctional cellular mechanisms.

The proposed model comprises an electrophysiological model coupling ICC and SMC via gap junctions, an active tension model bridging electrophysiology and biomechanics, and a model simulating neural regulation through nitrergic, purinergic, and cholinergic pathways, as illustrated in figure 1.

Zoom In Zoom Out Reset image size

The Lees-Green small intestinal ICC model [18] was reparameterised for gastric SWs. Within this model, the Ano1, NSCC, and IP₃R components (figure 1) were modified to incorporate enteric neural regulation. The Corrias and Buist [20] gastric SMC model was used to simulate intracellular SMC Ca²⁺ concentration and new SK and active tension components were formulated to extend this model and incorporate enteric neural regulation of SMC.

2.1. Interstitial cell of Cajal model

2.2. Smooth muscle cell (SMC) model

The SMC model includes a variety of membrane currents to model gastric SMC electrophysiology, with SW events driven by current injected via a gap junction (figure 1). The original implementation included a prescribed ICC membrane potential to set the gap junction current. In the present work, this prescribed ICC potential was instead directly coupled to the ICC model via a gap junction current, I_couple,

$$\begin{equation}{I_{{\text{couple}}}} = {g_{{\text{couple}}}} \cdot \left( {{V_{{\text{ICC}}}} - {V_{{\text{SMC}}}}} \right)\end{equation} \tag{ 1 }$$

where g_couple is the gap junction conductance, V_ICC is the ICC membrane potential, and V_SMC is the SMC membrane potential.

2.3. Active tension model

A Hill-type relationship was used to model the Ca²⁺-force relationship,

$$\begin{equation}T = {T_0} \cdot \frac{{{\text{C}}{{\text{a}}_{\text{i}}}^{\text{h}}}}{{{\text{C}}{{\text{a}}_{\text{i}}}^{\text{h}} + {\text{C}}{{\text{a}}_{{\text{50}}}}^{\text{h}}}}\end{equation} \tag{ 2 }$$

where T is the active tension, Ca_i denotes the intracellular calcium concentration of the SMC model, T₀ is the maximal tension, h denotes the Hill sensitivity coefficient, and Ca₅₀ denotes the concentration of cytosolic Ca²⁺ at half maximal tension. Ca₅₀ was modified by inhibitory neural regulation (see equation (16)) and the values for T₀ and h (313 kPa and 3.4, respectively) were obtained from isometric tension generation measurements using skinned smooth muscle strips from guinea pig taenia coli as reported by Arner [43]. Mechano-sensitive mechanisms were not incorporated into this model.

2.4. Neural regulation model

Efferent neural regulation by the stimulation of intrinsic motor neurons was included in this model, but reflex circuitry was not considered. Neural stimulation was modelled by variables for inhibitory and excitatory stimulation f_i and f_e. These variables represent the frequency of electrical field stimulation between 0 and 10 Hz for inhibitory nitrergic and purinergic stimulation, and excitatory cholinergic stimulation, respectively. The stimulation of intrinsic motor neurons is modelled as having a steady-state response, in concert with results from in vitro experiments where electrical field stimulation pulses of 0.1 ms pulse width were applied at a fixed frequency [36, 37]. Terms were formulated to model inhibitory effects (S_iX) at the effector components, Ano1 (S_iAno1), NSCC (S_iNSCC), Ca²⁺ sensitivity (S_iCa50), and SK (S_iSK),

$$\begin{equation}{S_{{\text{iAno1}}}} = {k_{{\text{iAno1}}}} \cdot {w_{{\text{iICC}}}}\end{equation} \tag{ 3 }$$

$$\begin{equation}{S_{{\text{iNSCC}}}} = {k_{{\text{iNSCC}}}} \cdot {w_{{\text{iICC}}}}\end{equation} \tag{ 4 }$$

$$\begin{equation}{S_{{\text{iCa50}}}} = {k_{{\text{iCa50}}}} \cdot {w_{{\text{iSMC}}}}\end{equation} \tag{ 5 }$$

$$\begin{equation}{S_{{\text{iSK}}}} = {k_{{\text{iSK}}}} \cdot {w_{{\text{iSMC}}}}\end{equation} \tag{ 6 }$$

where k_iAno1, k_iNSCC, k_iCa50, and k_iSK are inhibitory effector component weighting terms. The excitatory stimulation effect, S_e, affecting IP₃ concentration was modelled by,

$$\begin{equation}{S_\text{e}} = {k_{\text{eIP3}{\text{ }}}} \cdot {w_\text{e}}\end{equation} \tag{ 7 }$$

where k_eIP3 is an excitatory effector component weighting term.

Scaling terms w_iICC and w_iSMC, for ICC and SMC based inhibitory effector components respectively, and w_e for the excitatory effector component were calculated as functions of inhibitory (f_i) and excitatory (f_e) stimulation frequency

$$\begin{equation}{w_{{\text{iICC}}}} = \frac{{1 - \frac{{{e^{ - {p_{{\text{iICC}}}}{f_{\text{i}}}}}}}{{{f_{\max }}}}}}{{1 - {e^{ - {p_{{\text{iICC}}}}}}}}\end{equation} \tag{ 8 }$$

$$\begin{equation}{w_{{\text{iSMC}}}} = \frac{{1 - {e^{ - \frac{{{p_{{\text{iSMC}}}}{f_{\text{i}}}}}{{{f_{\max }}}}}}}}{{1 - {e^{ - {p_{{\text{iSMC}}}}}}}}\end{equation} \tag{ 9 }$$

$$\begin{equation}{w_{\text{e}}} = \frac{{1 - {e^{ - \frac{{{p_{\text{e}}}{f_{\text{e}}}}}{{{f_{\max }}}}}}}}{{1 - {e^{ - {p_{\text{e}}}}}}}\end{equation} \tag{ 10 }$$

where f_max is the maximum modelled input stimulation frequency fitted from experimental values [36, 37]. For fitting purposes, the frequency of electrical field stimulation during continuous stimulation on the order of minutes was used, in conjunction with data from pharmacological stimulation matched in tissue response to electrical field stimulation. The parameters, p_iICC, p_iSMC, and p_e determine the non-linearity of the relationship between stimulation frequency and the effect of stimulation.

The value of k_eIP3 and p_e were selected to reproduce the results from Forrest et al [36]. Values for the inhibitory effector components (k_iAno1, k_iNSCC, k_iCa50, and k_iSK) were selected by an optimisation procedure (see section 2.6). Fitting of the excitatory pathway parameters was excluded from the optimisation procedure because the experimental data for fitting inhibitory pathway parameters used pharmacological methods to block cholinergic neurotransmission. Conversely, data used to fit the excitatory pathway parameters was obtained under pharmacological block of nitrergic neurotransmission, so the excitatory parameters could be estimated in isolation from the inhibitory parameters.

2.4.1. Nitrergic pathway

Inhibition of the Ano1 channel by cGMP occurs due to a decreased current during Ca²⁺-dependent activation. To model this, an inhibition term was formulated as

$$\begin{equation}{h_{{\text{Ano}}1}} = 1 - \frac{{{S_{{\text{iAno1}}}}}}{{1 + {{\left( {d_{{\text{Ano}}1}^{^{^{\prime}}}} \right)}^2}}}\end{equation} \tag{ 11 }$$

where h_Ano1 is the Ano1 inhibition term, d'_Ano1 is the time derivative of the Ano1 gating variable, d_Ano1. This term inhibits activation during the calcium-activated plateau phase of the Ano1 current. Hence, the modified Ano1 current, I_Ano1, in the ICC model was given by

$$\begin{equation}{I_{{\text{Ano1}}}} = {g_{{\text{Ano1}}}} \cdot \left( {{d_{{\text{Ano1}}}} - {h_{{\text{Ano1}}}}} \right) \cdot \left( {{V_{{\text{ICC}}}} - {E_{{\text{Cl}}}}} \right)\end{equation} \tag{ 12 }$$

where g_Ano1 is the maximal conductance of the Ano1 channel, V_ICC is the ICC membrane potential and E_Cl is the reversal potential for Cl⁻.

Additionally, the role of NSCC as a potential effector of nitrergic stimulation was investigated by modulating the channel conductance, g_NSCC, which is applied to calculate the NSCC channel current, I_NSCC, defined by the following equations:

$$\begin{equation}{I_{{\text{NSCC}}}} = {g_{{\text{NSCC}}}} \cdot {P_{{\text{NSCC}}}} \cdot \left( {{V_{{\text{ICC}}}} - {E_{{\text{NSCC}}}}} \right)\end{equation} \tag{ 13 }$$

$$\begin{equation}{g_{{\text{NSCC}}}} = {g_0} \cdot \left( {1 - {S_{{\text{iNSCC}}}}} \right)\end{equation} \tag{ 14 }$$

$$\begin{equation}{P_{\text{NSCC}}} = \frac{{Ca_i^{{n_{\text{NSCC}}}}}}{{Ca_i^{{n_{\text{NSCC}}}} + Ca_{\text{NSCC}}^{{n_{\text{NSCC}}}}{\text{ }}}}\end{equation} \tag{ 15 }$$

where P_NSCC is a Ca²⁺-dependent gating variable, E_NSCC is the reversal potential, g₀ is the maximum NSCC conductance, Ca_i is the intracellular calcium concentration, n_NSCC is the Hill coefficient, and Ca_NSCC is the half-maximal intracellular Ca²⁺ concentration for gating activation.

Additionally, the effect of nitrergic stimulation on Ca²⁺ sensitivity in SMC was also investigated, by shifting the Ca²⁺-force relationship modelled in equation (2), based on the level of inhibitory stimulation. Ca₅₀ was calculated by

$$\begin{equation}{\text{C}}{{\text{a}}_{50}} = {\text{C}}{{\text{a}}_0} \cdot \left( {1 + {S_{{\text{iCa50}}}}} \right)\end{equation} \tag{ 16 }$$

where the baseline calcium sensitivity, Ca₀, was 0.56 μM [43].

2.4.2. Purinergic pathway

A new SK channel was added to the SMC model by adapting a model of ventricular myocyte SK channels developed by Kennedy et al [44]. The SK channel current was modelled by

$$\begin{equation}{I_{{\text{SK}}}} = {g_{{\text{SK}}}} \cdot {d_{{\text{SK}}}} \cdot {S_{{\text{iSK}}}} \cdot \left( {{V_{{\text{SMC}}}} - {E_{\text{K}}}} \right)\end{equation} \tag{ 17 }$$

where I_SK is the membrane SK channel current, g_SK is the maximal channel conductance, E_K is the K⁺ reversal potential, and d_SK is a Ca²⁺-dependent gating variable given by

$$\begin{equation}\frac{{{\text{d}}{d_{{\text{SK}}}}}}{{{\text{d}}t}} = \frac{{{d_{{\text{SK}}\infty }} - {d_{{\text{SK}}}}}}{{{\tau _{{\text{SK}}}}}}\end{equation} \tag{ 18 }$$

$$\begin{equation}{\tau _{{\text{SK}}}} = \frac{1}{{0.047{\text{C}}{{\text{a}}_{\text{i}}} + 0.013}}\end{equation} \tag{ 19 }$$

$$\begin{equation}{d_{{\text{SK}}\infty }} = \frac{{0.81{\text{C}}{{\text{a}}_{\text{i}}}}}{{{\text{Ca}}_{\text{i}}^{{n_{{\text{SK}}}}} + {\text{Ca}}_{{\text{SK}}}^{{n_{{\text{SK}}}}}}}\end{equation} \tag{ 20 }$$

where d_SK∞ is the scaling gating constant, τ_SK is the gating time constant, n_SK is the Hill coefficient of the scaling constant, and Ca_SK is the half-maximal Ca²⁺ concentration of the scaling gating constant.

2.4.3. Cholinergic pathway

Cholinergic stimulation was modelled by modification of the IP₃ production rate in the ICC model. Intracellular IP₃ production rate is given by

$$\begin{equation}{\text{I}}{{\text{P}}_3} = {\text{IP}}{3_0} \cdot \left( {1 + {S_{{\text{eIP3}}}}} \right)\end{equation} \tag{ 21 }$$

where IP3₀ is the baseline rate of IP₃ generation. The IP₃ production rate affects the following equations:

$$\begin{equation}{P_{{\text{IPR}}}} = {\left( {\frac{{{\text{I}}{{\text{P}}_3} \cdot {\text{C}}{{\text{a}}_{\text{i}}} \cdot \left( {1 - y} \right)}}{{\left( {{\text{I}}{{\text{P}}_3} + {K_1}} \right) \cdot \left( {{\text{C}}{{\text{a}}_{\text{i}}} + {K_5}} \right)}}} \right)^3}\end{equation} \tag{ 22 }$$

$$\begin{equation}{J_{{\text{IPR}}}} = \left({k_{{\text{IPR}}}}{P_{{\text{IPR}}}} + {J_{{\text{ER}}}}\right) \cdot \left( {{\text{C}}{{\text{a}}_{{\text{ER}}}} - {\text{C}}{{\text{a}}_{\text{i}}}} \right)\end{equation} \tag{ 23 }$$

$$\begin{equation}{\phi _1} = \frac{{{\text{C}}{{\text{a}}_{\text{i}}} \cdot \left( {{k_4} \cdot {K_2} \cdot {K_1} + {k_2} \cdot {K_4} \cdot {\text{I}}{{\text{P}}_3}} \right)}}{{{K_4} \cdot {K_2} \cdot \left( {{K_1} + {K_3}} \right)}}\end{equation} \tag{ 24 }$$

$$\begin{equation}{\phi _2} = \frac{{{k_2} \cdot {\text{I}}{{\text{P}}_3} + {k_4} \cdot {K_3}}}{{{K_3} + {\text{I}}{{\text{P}}_3}}}\end{equation} \tag{ 25 }$$

$$\begin{equation}\frac{{{\text{d}}y}}{{{\text{d}}t}} = {\phi _1} \cdot \left( {1 - y} \right) - {\phi _2} \cdot y\end{equation} \tag{ 26 }$$

$$\begin{equation}\frac{{{\text{dC}}{{\text{a}}_{\text{i}}}}}{{{\text{d}}t}} = {f_{\text{c}}} \cdot \left( {{J_{{\text{IPR}}}} - {J_{{\text{SERCA}}}} + {J_{{\text{SOC}}}} - {J_{{\text{PMCA}}}}} \right)\end{equation} \tag{ 27 }$$

where J_IPR is the calcium release flux due to IP₃-activated Ca²⁺ channels (IP₃Rs), k_IPR is a scaling constant, P_IPR is the Ca_i dependent gating variable for IP₃R, J_ER is the leak rate, Ca_ER is the calcium concentration within the intracellular store, IP₃ is the intracellular IP₃ concentration, y is the binding variable for IP₃ to IP₃R, ϕ₁ and ϕ₂ are the binding and unbinding probability respectively for IP₃ to IP₃R, based on the rate constants K₁ − K₅ (table 1), f_c scales the intracellular calcium level, and J_x are the rates of calcium flux for the channels IP₃R, SERCA, SOCE, and PMCA (figure 1).

2.5. Model simulation and metric quantification

2.6. Parameter optimisation

To select parameter values for the four inhibitory weighting terms k_iAno1, k_iNSCC, k_iCa50, and k_iSK and two non-linear scaling terms p_iICC and p_iSMC, a constrained, non-linear optimisation problem was solved using the interior-point algorithm ('fmincon' in MATLAB 2022b) to minimise the objective function,

$$\begin{equation}o = \mathop \sum \limits_{n = 1}^4 \left| {\frac{{M\left( {\vec z,{\text{ }}{f_n}} \right)}}{{{M_0}}} - E\left( {{f_n}} \right)} \right|\end{equation} \tag{ 28 }$$

where M is a function representing the simulated value of the optimised quantity, while E is a function representing the percentage change in quantity at a stimulation frequency, f_n , as determined by experimental work. The parameter values for the simulation, were captured in the vector $\vec z$. M₀ is the simulated plateau tension amplitude for f_i = 0 Hz (i.e. no stimulation) with input parameters $\vec z$.

The fitting procedure was performed separately for parameters in the ICC (Step 1) and SMC (Step 2) models, as described in table 2. First, the parameters in the ICC model were fit to data for the membrane potential change during inhibitory stimulation. Then, with the ICC parameters fixed, the parameters for the Ca²⁺ sensitivity, SK channel, and p_iSMC were fit to experimental data measuring the change in tension during phasic contractions under inhibitory stimulation. The optimisation procedure was performed 48 times using initialisations randomly sampled from a Sobol set of the feasible parameter space (see results section). Using 48 initialisations improved coverage of the parameter space and allowed for the identification of potential local minima. The median value of the solutions to Step 1 had a marginal change of less than ±2% in all parameters when more than approximately 24 initialisations were used after randomly permuting the initialisation sequence 100 times. Therefore, twice this number, 48, was deemed to be a sufficient number of initialisations. After 48 initialisations the marginal change in median value was less than 0.35% for all three ICC parameters. The values of k_eIP3 and p_e were selected based on experimental measurements of SW frequency change during cholinergic stimulation (see section 3).

Table 2. Steps in the parameter optimisation procedure.

	Step 1 (ICC)	Step 2 (SMC)
Parameters optimised	${{\text{k}}_{{\text{iAno}}1}},{\text{ }}{{\text{k}}_{{\text{iNSCC}}}},{\text{ }}{{\text{p}}_{{\text{iICC}}}}$	${{\text{k}}_{{\text{iCa}}50}},{\text{ }}{{\text{k}}_{{\text{iSK}}}},{\text{ }}{{\text{p}}_{{\text{iSMC}}}}$
Other parameters	${k_{{\text{iCa}}50}} = {k_{{\text{iSK}}}} = 0$, ${p_{{\text{iSMC}}}} = 1$	${{\text{k}}_{{\text{iAno}}1}},{\text{ }}{{\text{k}}_{{\text{iNSCC}}}},{\text{ }}{{\text{p}}_{{\text{iICC}}}}{\text{ }}$from the solution of Step 1 (see results section 3)
Experimental data	Figures 1(B) and 2(E) from Kim et al [37]	Figure 1(B) from Kim et al [37]
Quantity optimised	Slow wave plateau amplitude	Phasic contraction plateau amplitude
Constraints (further explanation in results)	$1.11{\text{ }}{{\text{k}}_{{\text{iAno}}1}} + {\text{ }}{{\text{k}}_{{\text{iNSCC}}}} \unicode{x2A7D} 1.41$, $0 < {k_{{\text{iNSCC}}}} < 1$, $0 < {k_{{\text{iNSCC}}}} < 1$	$0 < {k_{{\text{iCa}}50}} < 3$, $0 < {k_{{\text{iSK}}}} < 1$

2.7. Sensitivity analysis

The effect of perturbing the effector component weighting parameters (k_iAno1, k_iNSCC, k_iCa50, and k_iSK) between the 10 and 90th percentile of the optimisation procedure solution was investigated. Additionally, the excitatory component weighting (k_eIP3) was also perturbed between the minimum and maximum estimates of the parameter value. The Saltelli method was used to obtain a representative, uniformly distributed sample of this parameter space [45] with 1024 sampled points per analysed variable. Sampling and analysis were performed in Python using the 'SALib' module, while simulations were run and results quantified in MATLAB as described in section 2.5 [46]. Maximal inhibitory and excitatory stimulation (f_i = f_e = 10 Hz) was set for each sampled parameter combination. The first-order Sobol sensitivity (S₁) index and the total Sobol sensitivity (S_T) index for a given variable X_i were calculated as follows

$$\begin{equation}{S_1} = \frac{{{V_{{X_i}}}\left( {{E_{{X_{\sim i}}}}\left( {Y{\text{|}}{X_i}} \right)} \right)}}{{V\left( Y \right)}}\end{equation} \tag{ 29 }$$

$$\begin{equation}{S_{\text{T}}} = \frac{{{E_{{X_{\sim i}}}}\left( {{V_{{X_i}}}\left( {Y{\text{|}}{X_{\sim i}}} \right)} \right)}}{{V\left( Y \right)}}\end{equation} \tag{ 30 }$$

where V represents the variance, E represents the expected value, X_∼i is all variables except the ith variable, and Y is the metric value.

Stable solutions were obtained for V_ICC and V_SMC with a self-oscillatory SW frequency of 3.2 cpm after simulating for 300 s beginning at the initial conditions, with the SW interval varying by no more than 3 ms. This resulted in phasic contractions at the same frequency (figure 2(A)). The resting membrane potential, initial peak depolarisation amplitude potential and plateau amplitude potential of each SW event was consistent over time for both the ICC (−65.3 mV, 26.2 mV, and 36.6 mV respectively) and SMC (−65.7 mV, 16.0 mV, and 17.6 mV respectively). The phasic contraction amplitude was also consistent over the 300 s window, with a basal tension of 1.6 kPa, an initial peak tension amplitude of 43.3 kPa and a plateau amplitude of 41.7 kPa. The resting membrane potentials and basal tension varied by less than 0.0001 mV and less than 0.0001 kPa respectively over a 300 s simulation. Similarly, over the same period the initial peak amplitude of V_ICC, V_SMC, and peak tension varied by 0.010 mV, 0.002 mV, and 0.007 kPa respectively. SW events were synchronised but not necessarily time-aligned between tension, SMC membrane potential, and ICC membrane potential.

Zoom In Zoom Out Reset image size

Of the four inhibitory effector components, only two of these affect membrane currents in the ICC: Ano1 and NSCC. Varying the respective weightings of these components (k_iAno1, k_iNSCC) in the inhibitory stimulation pathway showed that when both components were non-zero, the effect on the ICC SW plateau was much greater than if only one was non-zero (figure 2(B)). However, the initial peak depolarisation was not affected by k_iAno1 and k_iNSCC. Furthermore, when both components had a high weighting, a non-physiological SW signal morphology exhibiting a second depolarisation resulted (figure 2(B) right). By performing a parameter sweep for the values of k_iAno1 and k_iNSCC (between 0 and 1 with step 0.05) the constraint 1.11 k_iAno1 + k_iNSCC ⩽ 1.41 was identified for Step 1 of the optimisation problem. This excludes the combinations of parameter values which gave a non-physiological SW signal morphology.

Inspecting the change in plateau tension for varying weightings at maximal stimulation (figure 3(A)) showed that the Ano1 and NSCC inhibitory effector components, had a much smaller impact on tension compared to the effect of the SK channel (k_iSK) and Ca²⁺ sensitivity (k_iCa50). Additionally, excitatory stimulation in isolation had relatively little effect on ICC SW amplitude regardless of k_eIP3. In contrast, excitatory stimulation via the IP₃ effector component was the only component that had a large effect on SW frequency (figure 3(B)). The inhibitory effector components of Ano1 and NSCC both had small depressive effects on SW frequency, while other components had no impact on SW frequency.

Zoom In Zoom Out Reset image size

The value of k_eIP3 was based on the results reported by Forrest et al [36]. During in vitro electrical field stimulation at 5 Hz and 10 Hz the observed SW frequency was 5.2 cpm and 5.4 cpm, respectively. The similarity of the effect at different neural stimulation intensities was modelled by setting p_e to 5. The selected k_eIP3 value was 1.0 to match the experimentally reported 62% increase in frequency attributable wholly to cholinergic electrical field stimulation at 5 Hz. On the other hand, at 10 Hz the response was not wholly attributable to cholinergic stimulation, and the value matching this response (0.82; figure 3(B)) was treated as the lower bound of the k_eIP3 variable for sensitivity analysis.

Solving Step 1 of the parameter optimisation procedure yielded median values of 0.31, 0.78, and 3.16 for k_iAno1, k_iNSCC, and p_iICC respectively across the 48 initialisations, with two clusters of solutions, as shown in figure 4. To define the threshold between clusters in an unbiased manner a Gaussian kernel probability density estimate of the parameter value for k_iAno1 was calculated with a half-width of 0.05. The parameter value with a minimal probability between the modal values was selected as the threshold (k_iAno1 = 0.20). One of the clusters gave Ano1 almost no weighting (k_iAno1 = 0.09, k_iNSCC = 0.84, p_iICC = 3.34) and the other gave Ano1 an effect, but the effect of NSCC was slightly diminished (k_iAno1 = 0.33, k_iNSCC = 0.77, p_iICC = 3.15). The median (lower quartile, upper quartile) of the objective function (equation (28)) for the percentage change in ICC SW plateau amplitude was 2.66% (2.53%, 3.37%) overall, 2.45% (2.43%, 2.49%) for the low k_iAno1 cluster, and 2.66% (2.65%, 2.66%) for the high k_iAno1 cluster. Due to experimental evidence demonstrating the role of the Ano1 channel in IJPs the high k_iAno1 cluster (k_iAno1 = 0.33, k_iNSCC = 0.77, p_iICC = 3.15) was selected and carried forward to the next step [5, 42]. Solving Step 2 of the optimisation procedure resulted in median values of 0.40, 0.30, and 0.12 for k_iCa50, k_iSK, and p_iSMC respectively. The median (lower quartile, upper quartile) of the objective function (equation (28)) for the percentage change in phasic contraction plateau amplitude was 0.18% (0.16%, 0.20%). The large positive value of p_iICC indicates that the marginal effect of stimulation decreases with stimulation intensity for effector components in ICC. Conversely, the small value of p_iSMC indicates that a nearly linear shift in Ca₅₀ and SK channel conductance modelled the response of SMC to inhibitory stimulation. Conducting Step 2 of the optimisation procedure with the low k_iAno1 cluster values resulted in fitted similar values for k_iCa50, k_iSK, and p_iSMC.

Zoom In Zoom Out Reset image size

Inhibitory stimulation with the fitted parameters (figure 5) resulted in decreased plateau tension while excitatory stimulation resulted in a small increase in plateau tension. When both inhibitory and excitatory stimulation were applied equally there was an overall decrease in plateau amplitude, to a similar extent as that with just inhibitory stimulation (−77% at f_i = 10 Hz, f_e = 10 Hz vs −79% at f_i = 10 Hz, f_e = 0 Hz), indicating an overriding effect of the inhibitory stimulation. Frequency was greatly increased by excitatory stimulation, while inhibitory stimulation resulted in a small decrease in frequency (−10% at f_i = 10 Hz, f_e = 0 Hz). Simultaneous inhibitory and excitatory stimulation resulted in an overall increase in frequency, but this was a smaller increase than with only excitatory stimulation (+47% at f_i = f_e = 10 Hz vs +66% at f_i = 0 Hz, f_e = 10 Hz).

Zoom In Zoom Out Reset image size

Finally, the variance-based model sensitivity analysis (figure 6) notes the relative sensitivity of the model outputs to the fitted values of the effector component weights. This depends on the variance of the model outputs within the tested interval around the fitted parameter value. We see the variance in simulated tension plateau amplitude can be primarily attributed to perturbation of k_iSK (S₁ = 0.60). Meanwhile the variance in simulated frequency was almost entirely attributable to perturbation of k_eIP3 (S₁ = 0.99). Simulating the model at the extreme values of these parameters in the tested sensitivity range with other parameters held at their fitted values shows that the absolute change in plateau tension amplitude decreased only 0.47% between for a 19% change in k_iSK from 0.293 to 0.343, while SW frequency increased only 4.9% between for an 18% change in the value of k_eIP3 from 0.82 to 1.

Zoom In Zoom Out Reset image size

A novel model of gastric electrophysiology and motility incorporating neural regulation was developed. In conjunction, the effect of neural regulation on ICC SWs and phasic contractions was simulated. Key experimental literature was curated and used to inform and validate the model. The simulation results reinforce the relevance of SMC for inhibitory neural regulation. As a corollary, the model predicted a diminished role for Ano1 in influencing phasic contraction characteristics during neural stimulation, which contrasts with the clear role of Ano1 in generating slow IJPs. Simulated excitatory stimulation increased SW frequency but had a limited effect on phasic contraction amplitude.

The high value of the fitted parameter for SMC contractility, k_iCa50, is consistent with observations in experimental literature where the change in SW amplitude was not commensurate with the much larger change in tension amplitude, measured simultaneously in vitro by Ozaki et al [47]. The sensitivity of the SMC contractile apparatus to Ca²⁺ is likely to be regulated by effects on myosin light chain phosphatase (MLCP) concentration via the cGMP signalling pathway [34]. A state-based model of active tension generation in gastric smooth muscle, such as that developed by Gajendiran and Buist [48], could be used to investigate this mechanism further.

Excitatory stimulation at SMC was not included in the model because experiments by Ozaki et al [49] showed that despite the presence of mechanisms that could transduce the effects of ACh on MLCP, no such relation was observed in the canine gastric antrum. This is supported by the lack of cholinergic excitatory junction potentials induced in the distal stomach of guinea pigs and rats [50, 51]. Furthermore, electrical field stimulation of mouse antral tissue deficient in ICC-IM did not yield the excitatory response seen in wild-type tissue both when stimulated for minutes, and when stimulated by short trains of stimulation [36, 52]. This indicates that the direct effect of the excitatory neural transmission is exerted at ICC and effects at SMC are indirect.

In the model, inhibitory stimulation at ICC had a relatively small impact on the amplitude and frequency of phasic contractions, compared to effects via SMC. However, experimentally observed changes in ICC SW characteristics would require changes at ICC, so the Ano1 and NSCC effector components were modelled. This is supported by the observation that IJPs involve Ano1, which is present in ICC but not SMC [53]. However, the cellular mechanisms generating IJPs need not be the same as those yielding a change in SW characteristics. The fitted parameters (figure 4(A)) suggest that of the two putative ICC effector components, k_iNSCC was relatively more important than k_iAno1 in explaining the observed change in SW given inhibitory stimulation. The identity of channels contributing to the NSCC is yet to be identified in stomach tissue, though candidates include channels of the TRP family, some of which exhibit regulation by cGMP [25]. TRP family channels have been associated [24, 54] with ICC, and some evidence of TRPM7 has been found in human small intestinal ICC [55], though their role is controversial [24]. On the other hand, Ano1 channel inhibition, particularly in ICC-IM may play a role due to effects of NO on intracellular Ca²⁺ concentration. As discussed further by Suzuki et al [56], and more recently Sanders and Ward [5], the exact mechanism of NO mediated neural regulation that results in the effects on motility, rather than IJPs, is unclear.

The result observed in figure 5, where simultaneous inhibitory and excitatory stimulation resulted in an overall inhibitory effect is in consonance with observations made by Ozaki et al [57] using chemical stimulation with sodium nitroprusside and ACh, as well as the interaction of junction potentials in colonic tissue reported by Furness [58]. There was little change in the amplitude of phasic contractions during purely excitatory stimulation. This contrasts with experimental observations by Zhang et al [59], in which excitatory stimulation resulted in an increase in basal tone and phasic contraction amplitude. The basal tone did not increase during excitatory stimulation in the model. Taken together, these observations and previous studies suggest the inhibitory stimulation had an overriding effect on phasic contractions. This could have implications on the design of future closed-looped neurostimulator protocols [60].

Two key limitations of the developed model are related to the lumping of ICC subgroups, and the assumption of steady state effects. The populations of ICC in the stomach have been further sub-classified as myenteric (ICC-MY) or intramuscular (ICC-IM). Experimental evidence suggests that ICC-MY are responsible for the generation of regular SWs, while ICC-IM are neuroeffector cells interspersed between ICC-MY and SMCs [61]. Electrical coupling of ICC-MY and ICC-IM to SMC is responsible for initiating the influx of Ca²⁺ which leads to active tension generation. ICC-IM were not explicitly modelled here, but the electrical coupling input to SMC models both the pacemaking and neurotransduction aspects of ICC function. Explicitly modelling ICC-IM and ICC-MY as separate but coupled cells would improve our ability to distinguish the cellular mechanisms required for gastric neural regulation. Furthermore, the intricate structural arrangements and biophysical relationships of different cell sub-classes of ICC and PDGFRα⁺ cells would need to be considered across extended spatial scales. This is a challenging task that requires sophisticated whole-mount imaging of transgenic animals or specific labelling of ICC sub-types accompanied by investigations of specific channels and channel control mechanisms [62].

Additionally, a candidate for the incorporation of further detail into the model is the time-course of neural regulation. In the present work, the interaction between motor neurons and gastric motility were simulated by modelling cellular mechanisms in ICC and SMC with a frequency-dependent, steady-state motor neuron stimulation input. The experimental data against which the model was evaluated used protocols involving either electrical or chemical stimulation at this time scale [36, 37, 57, 59]. However, some post-junctional effects of neural regulation, such as junction potentials, occur at a much shorter time scale. The stimulation intensity terms, f_i and f_e, used in the present model assume that prolonged stimulation at a given frequency results in a steady-state of functional response. In biological settings, this intensity is time-varying at a shorter time-scale and determined by neurotransmitter release characteristics including release frequency, release volume, and the number of endings as well as the physical transport of neurotransmitters across the junction. Consideration of shorter time scales would refine this aspect of the model and could help to further distinguish the cellular mechanisms responsible for post-junctional effects. To this end, it will be valuable in future to obtain sufficient experimental data to be able to model the transient aspects of transmission and events occurring post-stimulation.

The experimental parallel of this model lies in the data used for fitting, wherein electrical field stimulation was applied for many minutes at a fixed frequency, with pulse widths of 0.1 ms [36, 37]. Effects at a longer time-scale, as modelled in the present work, are particularly valuable for extension into whole organ modelling applications related to clinical applications of neuromodulation. Direct electrical pacing stimulation of the porcine stomach has shown that spatiotemporal entrainment of gastric bioelectrical activity at the organ-scale occurs over many SW cycles, on the order of minutes [63].

Simulation of neural regulation at the whole-organ level is particularly relevant in the context of peripheral nerve stimulation research. Models of intestinal motility at the organ level have been used to successfully simulate the transit of colonic and small intestinal luminal contents in conjunction with autonomic nervous system regulation [21, 64]. The neurons that influence gastric function operate as components of reflex circuits and work in conjunction with other influences, such as muscle stretch and hormones [1]. A comprehensive model of the gut-brain axis would be ideal but not tractable in the present study due to its complexity and insufficient understanding of its mechanisms.

Although pathological states were outside the scope of the present work, we believe that our model will provide a basis to investigate pathological changes in gastric electrophysiology and function caused by perturbations of specific mechanisms in our model. Future work bridging the gap between the peripheral nerves, such as the vagus nerve, and enteric neurons will be valuable for applying this model in the context of peripheral nerve stimulation. Particular phenomena, such as vago-vagal reflexes, operate through these pathways and control gastric function unrelated to the phasic contractions of the distal stomach [2]. The vagus nerve influences the gastric portion of the ENS, and recent interest in vagus nerve stimulation as a clinical therapy for a variety of peripheral organ ailments raises the question of how gastric motility is affected by vagus nerve stimulation [15, 65]. Whole organ modelling, incorporating regionally divergent structural and functional properties of the stomach [66, 67], would result in tangible clinical outcomes stemming from the developed mathematical model.

A novel mathematical model of enteric neural regulation of gastric motility was developed and used to determine the mechanisms through which enteric neural influence is exerted. The parameter optimisation procedure predicted that the inhibition of SMC Ca²⁺ sensitivity is an essential part of the functional inhibitory response in this model. Additionally, the inhibition of Ano1 alone during inhibitory stimulation was unable to explain experimentally observed effects on ICC SW amplitude. Furthermore, it was determined that excitatory stimulation affecting ICC alone can explain experimentally observed changes in frequency but not in changes in amplitude. This model is primed for use in further investigations which extend to include vagus nerve influence on gastric ENS regulation and downstream effects on gastric motility at the tissue to whole organ scale.

This work was funded, in part, by the NIH SPARC: Virtual Stomach Grant (OT2OD030538), the Ministry of Business, Innovation and Employment's Catalyst: Strategic fund, and the Marsden Fund Council managed by Royal Society Te Apārangi. ONA was supported by a University of Auckland Doctoral Scholarship.

The data that support the findings of this study are openly available at the following URL/DOI: 10.6084/m9.figshare.23624622 [68].

Additionally, an interactive model simulation is available on o²S²PARC: https://osparc.io/study/1170826c-9a64-11ee-8b3c-02420a0b1983.

ONA, RA, ARC, LKC, and PD initially conceived the research problem. ONA and PD developed and implemented the mathematical model. ONA analysed the model, drafted the manuscript, and generated figures. All authors reviewed the mathematical model formulation, critically reviewed the manuscript, and approved the final manuscript.

The authors declare no competing interests.