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

Abstract Objective. Neural regulation of gastric motility occurs partly through the regulation of gastric bioelectrical slow waves (SWs) and phasic contractions. The interaction of the tissues and organs involved in this regulatory process is complex. We sought to infer the relative importance of cellular mechanisms in inhibitory neural regulation of the stomach by enteric neurons and the interaction of inhibitory and excitatory electrical field stimulation. Approach. A novel mathematical model of gastric motility regulation by enteric neurons was developed and scenarios were simulated to determine the mechanisms through which enteric neural influence is exerted. This model was coupled to revised and extended electrophysiological models of gastric SWs and smooth muscle cells (SMCs). Main results. The mathematical model predicted that regulation of contractile apparatus sensitivity to intracellular calcium in the SMC was the major inhibition mechanism of active tension development, and that the effect on SW amplitude depended on the inhibition of non-specific cation currents more than the inhibition of calcium-activated chloride current (kiNSCC = 0.77 vs kiAno1 = 0.33). The model predicted that the interaction between inhibitory and excitatory neural regulation, when applied with simultaneous and equal intensity, resulted in an inhibition of contraction amplitude almost equivalent to that of inhibitory stimulation (79% vs 77% decrease), while the effect on frequency was overall excitatory, though less than excitatory stimulation alone (66% vs 47% increase). Significance. The mathematical model predicts the effects of inhibitory and excitatory enteric neural stimulation on gastric motility function, as well as the effects when inhibitory and excitatory enteric neural stimulation interact. Incorporation of the model into organ-level simulations will provide insights regarding pathological mechanisms that underpin gastric functional disorders, and allow for in silico testing of the effects of clinical neuromodulation protocols for the treatment of these disorders.


Introduction
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][9][10][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 postjunctional events.Existing models of gastrointestinal electrophysiology or motility focus on either ICC [16][17][18][19] or SMC [20][21][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][24][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 apaminsensitive 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 2+ 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 3 production of due to muscarinic activation by acetylcholine (ACh).This causes increased Ca 2+ release from intracellular stores, through the IP 3 -activated Ca 2+ channel (IP 3 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 3 .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.

Methods
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.
The Lees-Green small intestinal ICC model [18] was reparameterised for gastric SWs.Within this model, the Ano1, NSCC, and IP 3 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 2+ concentration and new SK and active tension components were formulated to extend this model and incorporate enteric neural regulation of SMC.

Interstitial cell of Cajal model
The ICC model incorporates four ionic currents modelled by Hodgkin-Huxley formalisms with Ca 2+ or voltage dependence: T-type Ca 2+ channel (CaT); a lumped background K + channel (Kb); background Na channel (Nab); and NSCC.Additionally, a model of store operated calcium entry is modelled within a Ca 2+ microdomain for the Ca 2+ -dependent and voltage-dependent Ano1, Cl − anion channel.Note that myenteric ICC (ICC-MY) and intramuscular ICC (ICC-IM) are modelled as a single entity.
The ICC model was modified to simulate the frequency of SW activity in distal stomach rather than the small intestine.Lees-Green et al adopted calcium dynamics from a model of murine airway SMC [38], which were themselves adapted from the de Young-Keizer model of single pool calcium dynamics [39].As indicated in table 1, time constants for the IP 3 R channel were modified to simulate the SW frequency of the rat stomach.Additionally, to maintain physiological calcium levels, the maximal rate of Ca 2+ sequestration by sarcoendoplasmic reticulum calcium ATPase (SERCA) was increased.It was verified that the effect of Ano1 knockout, as simulated by the original ICC model, was maintained after making the specified changes to calcium dynamics.In the modified model, Ano1 knockout (modelled as Ano1 conductance set to 0 nS) results in small amplitude oscillations of ICC membrane potential (5 mV).Furthermore, E Cl was set to −11 mV to more closely align with experimentally measured values in cultured murine myenteric ICC [40].The role of a Ca 2+ -dependent chloride current (CaCC), such as that facilitated by Ano1, in maintaining the plateau potential of the ICC SW was identified by Kito et al [41].A decrease in Cl − conductance has been proposed as a mechanism for nitrergic inhibition in ICC due to its observed involvement in the IJP [42].

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 , where g couple is the gap junction conductance, V ICC is the ICC membrane potential, and V SMC is the SMC membrane potential.

Active tension model
A Hill-type relationship was used to model the Ca 2+force relationship, where T is the active tension, Ca i denotes the intracellular calcium concentration of the SMC model, T 0 is the maximal tension, h denotes the Hill sensitivity coefficient, and Ca 50 denotes the concentration of cytosolic Ca 2+ at half maximal tension.Ca 50 was modified by inhibitory neural regulation (see equation (16)) and the values for T 0 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.

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 Table 1.Parameters modified in the Lees-Green et al (2014) ICC model to adjust the Ca 2+ oscillation frequency.These parameters are based on the de Young-Keizer model [39] of single pool calcium dynamics.Note that K i = k i− /k i in equations ( 22), ( 24) and (25).SERCA: sarcoendoplasmic reticulum calcium ATPase, IP3: inositol triphosphate.

Variable description Variable Value
Receptor binding (µM Maximal SERCA rate (µM s −1 ) Ve 160 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 2+ sensitivity (S iCa50 ), and SK (S iSK ), where k iAno1 , k iNSCC , k iCa50, and k iSK are inhibitory effector component weighting terms.The excitatory stimulation effect, S e , affecting IP 3 concentration was modelled by, 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 ) 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 nonlinearity 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.

Nitrergic pathway
Inhibition of the Ano1 channel by cGMP occurs due to a decreased current during Ca 2+ -dependent activation.To model this, an inhibition term was formulated as 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 calciumactivated plateau phase of the Ano1 current.Hence, the modified Ano1 current, I Ano1 , in the ICC model was given by 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: where P NSCC is a Ca 2+ -dependent gating variable, E NSCC is the reversal potential, g 0 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 2+ concentration for gating activation.Additionally, the effect of nitrergic stimulation on Ca 2+ sensitivity in SMC was also investigated, by shifting the Ca 2+ -force relationship modelled in equation ( 2), based on the level of inhibitory stimulation.Ca 50 was calculated by where the baseline calcium sensitivity, Ca 0 , was 0.56 µM [43].

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 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 2+ -dependent gating variable given by 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 2+ concentration of the scaling gating constant.

Cholinergic pathway
Cholinergic stimulation was modelled by modification of the IP 3 production rate in the ICC model.Intracellular IP 3 production rate is given by where IP3 0 is the baseline rate of IP 3 generation.The IP 3 production rate affects the following equations: ) where J IPR is the calcium release flux due to IP 3activated Ca 2+ channels (IP 3 Rs), k IPR is a scaling constant, P IPR is the Ca i dependent gating variable for IP 3 R, J ER is the leak rate, Ca ER is the calcium concentration within the intracellular store, IP 3 is the intracellular IP 3 concentration, y is the binding variable for IP 3 to IP 3 R, φ 1 and φ 2 are the binding and unbinding probability respectively for IP 3 to IP 3 R, based on the rate constants K 1 − K 5 (table 1), f c scales the intracellular calcium level, and J x are the rates of calcium flux for the channels IP 3 R, SERCA, SOCE, and PMCA (figure 1).

Model simulation and metric quantification
The model was implemented in CellML, and solved using the 'ode15s' , variable step, variable order integrator in MATLAB 2022b (The MathWorks, Natick, MA, USA).The absolute and relative tolerances were both 1 × 10 −6 and the maximum step size was 1 ms.
To confirm model stability a 300 s simulation was conducted.Following this, in every case where analysis was performed, the system was solved for 60 s to ensure a steady state oscillatory solution and analysis was performed on a further 120 s of simulation time.Unless specified, stimulation was applied for the entire simulation.Metrics were quantified to characterise the steady state oscillatory solution, frequency in cycles per minute (cpm), initial peak tension (kPa), and plateau tension (kPa).Additionally, initial peak depolarisation and plateau potential of SMC or ICC membrane potentials (mV) were also quantified as required.
For both membrane potential and active tension, the initial peak was defined as the maximum value during a 100 ms period after the peak positive rate of change.The amplitude values are given as a difference from baseline (resting membrane potential or baseline tension).The plateau value was defined as the mean value during which the smoothed absolute derivative (moving mean: 200 ms window) was less than 2.5 mV s −1 during a 3 s window after peak depolarisation.The derivative threshold and window size for identifying the plateau period were selected by inspecting the result for extreme cases of parameter values.

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, 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 ⃗ z.M 0 is the simulated plateau tension amplitude for f i = 0 Hz (i.e.no stimulation) with input parameters⃗ 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 2+ 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  [37] Figure 1(B) from Kim et al [37] Quantity optimised Slow wave plateau amplitude Phasic contraction plateau amplitude Constraints (further explanation in results) 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).

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 firstorder Sobol sensitivity (S 1 ) index and the total Sobol sensitivity (S T ) index for a given variable X i were calculated as follows 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.

Results
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 timealigned between tension, SMC membrane potential, and ICC membrane potential.
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 2+ 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 3 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.
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, Figure 3. Plateau tension (A) and frequency (B) of phasic contractions when weighting parameters were individually varied from 0 to 1 in 0.1 increments, with f i = 10 Hz and fe = 10 Hz, non-varying constants set to 0, and non-linear scaling terms w iICC , w iSMC and we set to 1. Metrics were calculated over 120 s of simulation time, with stimulation applied throughout this period (see section 2.5).When k iAno1 , k iNSCC , k iCa50 , k iSK , and keIP3 were all set to 0 the frequency was 3.2 cpm and the plateau tension was 41.7 kPa, with the labelled percentages being the percentage change from this baseline when the varied parameter has value 1. Inhibition at ICC by Ano1 or NSCC was not sufficient to inhibit tension generation to the extent observed in experiments.Excitatory stimulation alone, via keIP3, did not have a large effect on peak tension generation, but it did have a large effect on frequency.
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 50 and 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).
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 1 = 0.60).Meanwhile the variance in simulated frequency was almost entirely attributable to perturbation of k eIP3 (S 1 = 0.99).Simulating the model at the extreme values of these parameters in the tested sensitivity

Discussion
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 2+ 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 2+ 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 2+ 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 steadystate 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 postjunctional 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 wholeorgan 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 gutbrain 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.

Conclusion
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 2+ 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.

Figure 1 .
Figure 1.Schematic representation of the coupled electrophysiology model of a distal gastric interstitial cell of Cajal (ICC) and gastric smooth muscle cell (SMC), illustrating the modelled neural regulation pathways.The gastric SMC model developed by Corrias and Buist [20] was coupled to a reparameterised model of ICC originally developed by Lees-Green et al [18].The components shaded in orange are novel additions formulated in the present work.Inhibitory effector components are outlined in purple and excitatory effector components are outlined in green.CaL: L-type calcium channel, CaLVA: T-type calcium channel (see also CaT in the ICC model), SK: small conductance potassium channel, BK: large conductance potassium channel, Kb: lumped background potassium current, Kr: delayed rectifying potassium channel, Ka: A-type voltage-gated potassium channel, NSCC: non-specific cationic conductance, Na: lumped sodium current, NCX: sodium-calcium exchanger, SR: sarcoplasmic reticulum, Gαq: G-protein α subunit, ACh: acetylcholine, CaT: T-type calcium channel (see also CaLVA in the SMC model), PMCA: plasma membrane calcium ATPase, IP3R: IP3-activated calcium channel, SERCA: sarcoendoplasmic reticulum calcium ATPase, SOC: store-operated calcium channel, Ano1: anoctamin-1 voltage-gated, calcium activated chloride channel, cGMP: cyclic guanosine monophosphate second messenger nucleotide, NO: nitric oxide, Nab: lumped background sodium current.

Figure 2 .
Figure 2. (A) Membrane potential for ICC (top), SMC (middle) and active tension (bottom) simulated with no inhibitory or excitatory stimulation over 300 s.The resting membrane potential was stable and slow wave frequency was consistent throughout the simulation.(B) ICC membrane potential (top) and active tension (bottom) during a single slow wave cycle with varying k iAno1 alone (left), k iNSCC alone (centre), and k iAno1 and k iNSCC together with k iAno1 = k iNSCC (right) from 0 to 1 in 0.25 increments and f i set to 10 Hz.Arrows indicate increasing weighting (k iAno1 and k iNSCC ) values.Inhibition of Ano1 and NSCC alone does not give an appropriate decrease in plateau membrane depolarisation.Inhibition of both Ano1 plateau current and NSCC is required to yield a closer match to experimental data, but a complete block of both components (k iAno1 = k iNSCC = 1) results in non-physiological behaviour as indicated by the broken lines.

Figure 4 .
Figure 4. Results of the parameter optimisation procedure for effector component parameter values in two steps, showing a combined box plot and jittered scatter plot across the 48 initialisations.The selected parameter values are k iAno1 = 0.33, k iNSCC = 0.77, k iCa50 = 0.40, k iSK = 0.30.

Figure 5 .
Figure 5. Plateau tension (A) and frequency (B) for varying f i and/or fe from 0 to 10 Hz with step 0.05 for the optimised values of the inhibitory effector components ([k iAno1 , k iNSCC , k iCa50 , k iSK ] = [0.32,0.78, 0.4, 0.30]).(C) Traces of the SW (top) and phasic contractions (bottom) for a stimulation duration of 60 s indicated by the red line with only inhibitory (f i = 10 Hz), only excitatory; (fe = 10 Hz), or both inhibitory and excitatory (f i = fe = 10 Hz) stimulation applied.

Figure 6 .
Figure 6.Parameter sensitivity analysis for neural regulation parameters.The first-order (top row) and total (bottom row) Sobol sensitivity indices are shown for the tension plateau amplitude (left column) and frequency (right column) metrics during the stimulation period.Vertical black lines atop bars indicate the 95% confidence interval for the sensitivity index.

Table 2 .
Steps in the parameter optimisation procedure.