Physiological insights into electrodiffusive maintenance of gastric mucus through sensitivity analysis and simulations.

It is generally accepted that the gastric mucosa and adjacent mucus layer are critical in the maintenance of a pH gradient from stomach lumen to stomach wall, protecting the mucosa from the acidic environment of the lumen and preventing auto-digestion of the epithelial layer. No conclusive study has shown precisely which physical, chemical, and regulatory mechanisms are responsible for maintaining this gradient. However, experimental work and modeling efforts have suggested that concentration dependent ion-exchange at the epithelial wall, together with hydrogen ion/mucus network binding, may produce the enormous pH gradients seen in vivo. As of yet, there has been no exhaustive study of how sensitive these modeling results are with respect to variation in model parameters, nor how sensitive such a regulatory mechanism may be to variation in physical/biological parameters. In this work, we perform sensitivity analysis (using Sobol' Indices) on a previously reported model of gastric pH gradient maintenance. We quantify the sensitivity of mucosal wall pH (as a proxy for epithelial health) to variations in biologically relevant parameters and illustrate how variations in these parameters affects the distribution of the measured pH values. In all parameter regimes, we see that the rate of cation/hydrogen exchange at the epithelial wall is the dominant parameter/effect with regards to variation in mucosal pH. By careful sensitivity analysis, we also investigate two different regimes representing high and low hydrogen secretion with different physiological interpretations. By complementing mechanistic modeling and biological hypotheses testing with parametric sensitivity analysis we are able to conclude which biological processes must be tightly regulated in order to robustly maintain the pH values necessary for healthy function of the stomach.

Introduction

The surface of the gastric epithelium is covered in a continuous layer of mucus, an entangled, hydrated gel of polymeric proteins. For more than six decades, the physiology community has accepted that this gastric mucus layer provides a protective barrier, shielding the gastric wall (mucosa) from the stomach interior (lumen). In addition to preventing epithelial infection by most bacteria and other pathogens, the gastric mucus layer appears to protect the mucosa from the high acidity and digestive enzymes within the gastric cavity, preventing auto-digestion of the epithelium [1, 2]. In numerous mammalian species, the gastric mucus layer has been observed to support a massive (105-fold) hydrogen gradient across a relatively small length scale (several hundred μm) [3-5]. Failure to maintain this pH gradient across the layer (due to infection or misregulation) is associated with numerous pathologies, including gastric ulcers and cancer [6].

There are numerous unanswered questions regarding the mechanisms that give rise to the protective function of the mucus layer. Competing hypotheses exist regarding how acid (which is produced within crypts on the epithelial surface) is transported to the stomach lumen without acidifying the space immediately adjacent to the epithelium [5, 7, 8]. Unfortunately, the experimental evidence is inconclusive, and no consensus has been reached [9]. It has been known for some time that the the gastric epithelium also secretes bicarbonate, which can neutralize dissolved hydrogen ions and contributes to the protective nature of the mucus layer [10, 11]. However, there has been some debate if hydrogen, once transported to the stomach lumen, diffuses normally back down its own concentration gradient towards the epithelium, and if not, what mechanisms retard its diffusive flux [12]. Finally, there is no clear understanding of how the secretory processes within the epithelium are coordinated with the physical processes governing molecular diffusion through the mucus layer in order to maintain healthy gastric function. In this regard, mathematical models provide an invaluable tool, allowing one to analyze the interplay of physical and biological processes that give rise to the gastric pH gradient and often providing insights that may be inaccessible by experimental techniques.

It is known that the ionic composition of gastric secretions can vary in response to numerous stimuli, including histamine levels, the recency of a meal, and the contents of the gastric lumen [10, 13, 14]. At least one theoretical study has analyzed a model of how acid secretion is controlled by nervous and endocrine stimuli [15]. However, this model lacks any description of spatial variation in hydrogen concentration and therefore cannot explain the hydrogen gradient seen in vivo. Other, more physical, models have attempted to quantify the rate of hydrogen flux through gastric mucus and use this to calculate the secretion of bicarbonate required to maintain a healthy pH gradient [16]. However, experiments attempting to measure rates of bicarbonate secretion disagree with the model predictions [17].

To date, few modeling works have attempted to couple a physical description of diffusion through the mucus layer to a biologically informed model of secretion at the epithelial surface. One such model was put forth by Lewis, et. al in [18]. The underlying assumption made in this model is that hydrogen is transported from epithelium to lumen sequestered on mucin polymers (and therefore is not a dissolved, diffusive species during this process), and then released in the lumen by auto-digestion of the mucus network, as suggested by Schreiber, et. al [5]. A physical model of electro-diffusion through a complex gel-like material was coupled with non-linear boundary conditions that mimic the anti-port ion exchange proteins expressed by gastric epithelial cells [11]. In particular, bicarbonate secretion in this model functionally depends on the ionic concentration (chloride, specifically) in the space immediately adjacent to the epithelial surface. Analysis showed that such a model was capable of reproducing the pH gradient across the mucus layer observed in vivo. In fact, the pH gradient was maintained regardless of nearly all variations in model parameters, with one exception: knocking out the hydrogen/sodium exchange proteins in the epithelial surface (by setting a parameter to zero) destroyed or even inverted the pH gradient across the epithelial surface and caused the predicted wall pH to depend “strongly” on the other model parameters.

The sensitivity of model predictions to variation in parameters is an important issue in the realm of biological modeling. Lack of accuracy or precision in experimental measurements can lead to uncertainty in model parameter values. Quantitative assessment of the impact that this uncertainty has on model predictions is therefore a valuable tool. Furthermore, the physical quantities which model parameters represent may not be fixed in the biological system of interest (varying temporally, for example). The sensitivity (or insensitivity) of model output to variations in these values gives meaningful insight into the underlying biology. The analysis of [18] represents what is often termed “one-at-a-time” sampling in the field of sensitivity analysis: restriction to a lower dimensional hyperplane of parameter space. Not only do such lower-dimensional analyses risk missing important features of the system behavior, in some cases they may even lead to misleading results [19]. Because of this, we will analyize a crucial prediction of the mathematical model from [18] using a more robust, global measure of sensitivity.

In this work, we use Sobol’ Indices (SIs) to perform a global sensitivity analysis of the model of gastric epithelial pH regulation first put forth in [18]. In Section 2, we outline the mathematical model as well as the relevant parameters. Section 3 provides a discussion of the theory and calculation of SI. In Section 4, we report the SI calculated for the model parameters and discuss the dependence of epithelial pH distribution on various parameters. In particular, we see that the rates at which ion exchange events take place at the epithelial surface exert the strongest control over the distribution of epithelial pH, followed by the strength of hydrogen secretion. To further explore the dependence of mucosal pH on these “important” parameters, we perform several restricted Sobol’ analyses where specific parameters are fixed at extremal values. These parameter regimes replicate physiologically relevant situations (stimulated acid secretion, for example). The functional dependence of wall pH distribution on the remaining parameters is quantified, and their SI are calculated. Finally, we discuss the physiological significance of these results in Section 5.

Mathematical Model

Here, we utilize a mathematical description of diffusion through a gastric mucus layer that was first put forth in [18]. The model was developed based on the experiments of Schreiber et al. to mathematically investigate the proposed mechanism of hydrogen transport proposed in [5]. To our knowledge this represents the only physics-based model of hydrogen diffusion through mucus which has successful predicted the enormous pH gradient (from lumen to stomach wall) observed in mammalian stomachs. It is based upon the Nernst-Planck equations of electrodiffusion in conjuction with a “two-phase gel” representation of mucus. The two-phase gel framework has been used extensively to model gel-like substances in biological contexts [20, 21]. At each point in space there simultaneously exist two materials: the “gel” (cross-linked mucin polymer network), and the “solvent” (interstitial hydrating fluid). The local composition of the combined material is described by the dimensionless quantities θ and θ ≔ 1 − θ which represent the “volume fraction” of solvent and network, respectively (i.e. the proportion of local volume that is made up of each constituent material). In this work, we are primarily concerned with the diffusion of ionic species dissolved within the solvent phase. Therefore, the phases θ and θ) will be fixed in time, and their spatial profiles specified to represent the inhomogeneous material through which ions are moving (discussed at more length below).

The state variables of the model are the concentration of four ionic species and the electric potential gradient within the domain. Because the ions are dissolved in (and diffuse through) only the solvent phase of the gel, we define the concentration variables as moles per volume of solvent, not moles per total volume. Hydrogen and bicarbonate are denoted CH and CB, both with units of moles per liter (Molar). All other cations and anions are represented by CI and CA, respectively, measured in units of Equivalent per liter (which we denote M for notational convenience). The most common ions in the gastric interior are sodium (cations), potassium (cations), and chloride (anions) [22]. By including sodium and potassium in a single variable, we assume that both cations have the same chemical properties and diffusion constants. The spatial variable x measures distance from the epithelial wall of the stomach towards the lumen, and thus the location x = 0 corresponds to the epithelial surface. The concentration of ionic species is governed by a Nernst-Planck-like equation, which accounts for the flux of ions due to advection, diffusion, and electric potential gradients, as well as any reactions and sources that impact ionic concentration. A time-dependent version of this model was used in [18], and a careful derivation may be found in [23], however the present work is only concerned with steady state behavior. The equation which governs each ionic species may be described as the balance of fluxes and sources:

Here, the index j may take the values H, B, A, or I for hydrogen, bicarbonate, anions, and cations respectively. Terms I, II, and III represent fluxes of ionic concentration due to advective, diffusive, and electric effects repsectively. Term IV captures any reactions and/or sources that impact the concentration of species j. The quantity u(x) is the velocity of the fluid solvent (in cm/sec), which results from secretion of fluid at the stomach wall. The parameters D and z are the diffusion coefficient (cm2/sec) and valence (±1) of species j respectively. The variable Ψ is the non-dimensional electric potential (which may be converted to voltage by multiplying by RT/F, where R is the ideal gas constant, T is absolute temperature, and F is the Faraday constant). Note that Eq. (1) is similar to a standard steady state Nernst-Planck equation, with the exception of two appearances of θ(x) on the right hand side. These terms account for the fact that the ionic species are transported within the solvent phase of a mixture whose composition varies spatially (see [23]). The electric potential Ψ is determined by enforcing electro-neutrality, which is the principle that ionic concentrations may not result in a net charge at any location within space. Mathematically, this is expressed as

Thus, Ψ does not have its own equation but may be viewed as a Lagrange multiplier enforcing the constraint Eq. (2) on the set of evolution equations given by Eq. (1). We note here that the electroneutrality constraint (Eq. (2)) does not account for any charge on the mucus network itself. Though mucus polymers carry negative charge groups, the work of Schreiber, et al. indicates that in vivo mucus carries a large concentration of bound hydrogen ions, effectively neutralizing many of these negative charge groups [5]. Thus, the model implicitly assumes that approximately all of the negative charge on mucus is neutralized by bound cations.

Other works have explored the dynamic rearrangement of mucus in response to ionic diffusion using a similar modeling framework [24]. However, as we wish to focus on the steady state behavior of ionic species, we will treat both the solvent velocity and volume fraction as given model inputs. We assume that gastric mucus is constantly produced at the epithelial surface, moves away from the wall, and is degraded as it approaches the lumen in a manner that maintains a time-independent profile θ(x). We choose θ(x) to have a spatial profile that reflects a “standing front” defining the transition from the mucus gel layer to the stomach lumen. Figure 1 shows the spatial form of θ(x) (recall that θ(x) = 1 — θ(x)). The transition is located 500 μm away from the wall, and the mucus layer is indicated by the grey shaded region in Fig. 1. We define the “mucus layer” as the region where mucus network occupies more than than 1% of the local volume (i.e θ > 0.01, θ < 0.99). The “lumen” is therefore the region where θ < 0.01 and θ > 0.99. We note here that this model does not explicitly describe the swelling of mucus that happens on very short length and time scales immediately after secretion, nor does it incorporate the (relatively modest) changes in mucus volume fraction which may result from changes in monovalent cation concentration (as all ions in the model are monovalent) [25, 26]. Finally, the profile of u(x) is chosen to ensure that the quantity u(x)θ(x) is constant in space (as required by conservation of solvent mass) and approximates observed rates of fluid secretion from the gastric mucusa (see [18]). Choosing u(x) in this manner also ensures that Eq. (1) (in the absence of sources and reactions) is conservative in the quantities θ(x)C, which are the total ionic concentrations (moles per total volume).

Fig. 1:

Illustration of gel volume fraction (θ) and Hydrogen/Anion source profiles. Recall that θ = 1 – θ.

The terms g incorporate any sources or chemical reactions that may impact the four concentrations. Hydrogen is affected by the buffering reaction with bicarbonate (which proceeds according to the law of mass action). Because the dissociation constant of hydrogen and bicarbonate is small (4.2 × 10−7 M [16]), we assume that this reaction is non-reversible. In [5], Schreiber et. al proposed a mechanism of hydrogen transport in which hydrogen ions are also produced by a source located at the edge of the mucus layer. Mathematically, we represent this as a source S(x) supported in the same region where the volume fraction transitions from mucus layer to lumen. This is to model the release of hydrogen due to enzymatic degradation of the mucus gel itself, in accordance with the mechanism proposed by Schreiber et. al [5]. Because the gel layer is fixed, the source of hydrogen is assumed to be independent of time, but localized spatially,

In [18], a hydrogen source was used that represents the “average” rate of hydrogen secretion in the gastric mucosa. In this work, we investigate the effects of changes in hydrogen secretion on model predictions, and therefore our source is given by , where is the “baseline” hydrogen source used in [18]. The constant S0 is a parameter which controls the magnitude of hydrogen production. The profile of is shown in Fig. 1. The concentration of bicarbonate is affected by the same hydrogen/bicarbonate buffering reaction, hence

Electro-neutrality implies that any hydrogen ion released by mucus degradation at the interface must be accompanied by a corresponding negatively charged particle. Therefore, the source term also affects the concentration of anions:

Finally, we assume that cations are not affected by any sources or reactions:

Equation (1) must be accompanied by boundary conditions at the wall (x = 0) and the lumen interior (x = L). At the wall (x = 0), we impose a flux relation that depends on the local ionic concentrations at the epithelial surface. We do not explicitly specify the flux of any ion through the epithelial surface. For each species, the total flux from the transport equation is

When evaluated at the mucosa (x = 0), this expression must be equal to the flux through the ion-exchange proteins expressed by the epithelial cells that make up the mucosal wall. These epithelial cells express two ion exchange proteins that we account for here: one which exchanges bicarbonate and chloride (in a 1-to-1 ratio), and another which exchanges hydrogen and sodium (again, 1-to-1) [27]. Mathematically, our boundary conditions at the left boundary are given by

The terms on the right hand side represent a simplified, linear model of flux due to anti-port ion exchange proteins [28]. Concentrations at the stomach interior (which we define as x = 0.2 cm) are given by Dirichlet boundary conditions.

These values may not be chosen independently. Electroneutrality in the stomach interior requires that

Several model parameters are either not well known, or may actively change in physiological situations. S0 represents the magnitude of hydrogen production by the gastric mucosa, which is known to change in response to physiological stimuli [14, 22, 29]. The parameters kHI and kAB quantify the rate at which hydrogen/sodium and bicarbonate/chloride exchange takes place. These values are presumably a function of the density of the anti-port proteins expressed by the gastric epithelium, as well as the specific thermodynamic properties of a single exchange event. These quantities are difficult to estimate, and may not remain constant in vivo. Similarly, δHI and δAB represent a “bias” in the respective ion exchange proteins, which can be related to the concentration of individual ionic species within the epithelial cells [18]. Existing estimates for these values are based on data for various cell types and may not necessarily be applicable to gastric epithelial cells. Finally, H and AL and represent ionic concentrations within the bulk of the stomach, which our model treats as prescribed, constant values, but which are known to vary temporally in vivo [13]. The one exception to this is B, which is the lumenal concentration of bicarbonate. This value is essentially zero in all situations, and which we take to be 10−14 M (simply for numerical purposes). Therefore, the model has 7 parameters that may reasonably vary in vivo, or are poorly estimated in the literature: S0, kHI, kAB, δHI, δAB, H , and I (the electro-neutrality constraint determines A). The main purpose of this work is to perform sensitivity analysis of the model predictions with respect to these 7 parameters. This analysis will help determine which parameters play a significant role in the model prediction of [18]. These “sensitive” parameters require accurate estimates. Conversely, “insensitive” parameters indicate values whose variations have little effect on model predictions.

To do this, we must first decide on a “Quanitity of Interest” (QoI). As our QoI, we will use the negative logarithm (base 10) of the hydrogen concentration at the left boundary

This corresponds to the pH at the mucosal surface of the stomach, and can be interpreted as a measurement of gastric health (a neutral stomach wall with pH approximately 6 or 7 would indicate a healthy stomach, while significantly lower pH would indicate an unhealthy stomach). Therefore, one may view sensitivity analysis of this QoI as quantifying the sensitivity of “predicted stomach health” to variations in physiological parameters.

Sensitivity Analysis

To analyze the sensitivity of the non-linear relationship between the pH at the mucosal surface of the stomach and the seven parameters (implicitly defined by Eqs. (1) to (13)), we use Sobol’ indices (SI). In [30] Sobol’ defined and developed the theory for SI, variance-based sensitivity measures that do not assume linearity and monotonicity in the mathematical model. This makes them attractive to a large class of models across different fields, and they have since been used in engineering [31], physical systems [32, 33], biological systems [34], and economic models [35]. Moreover, simultaneous variation in all of the parameters is considered over the entire parameter space to compute SIs, providing a comprehensive exploration of the sensitivity of the output to the inputs. Hence, SIs are classified as global sensitivity analysis indices. The numerical methods that we used to estimate total and first-order SIs in this work, are based on Monte-Carlo estimation of integrals over a volume in the full-dimensional parameter space. As a result, it can be computationally expensive to estimate SIs if the number of parameters, or the dimensionality of the parameter space, is large [36, 37]. However, as we have previously mentioned, less computationally expensive local methods like one-at-a-time sampling can miss important parameter interactions, giving misleading conclusions [19]. The reader is directed to [38-40] for a review of various sensitivity measures.

Mathematical theory and notation

Given a scalar quantity of interest, y, as a square-integrable function of parameters {p1, …, p}, where p ∈ [a],

Sobol’ showed in [30] that there exists a unique functional decomposition, where f0 is a constant, I ⊂ D {1, …, N}, and f is a function of parameters {p} where i ∈ I, and

Using Eq. (15), it can be shown that, if I, J are two distinct subsets of D, then ∫ fp = 0. Hence, f and f are orthogonal if I ≠ J. In other words, Eqs. (14) and (15) show that a square-integrable function f can be decomposed into lower dimensional orthogonal functions. Sobol’ [30] also gives an analytical method to evaluate f recursively. Furthermore, if we consider the parameters to be independent and uniformly distributed over their prescribed interval, Sobol’ showed that the variance in the distribution of the output, resulting from the variation in the parameters, can be decomposed as, where Y is the random variable corresponding to the output y = f(p1, …, p) and ∫ · dp represents integral over the parameters p1, …, p. This gives an ANOVA-like decomposition of the variance of the output by partitioning it into contributions from uncertainties in the parameters, considering not only the effects of uncertainty in an individual parameter but also the uncertainty propagation due to functional dependence of the output on any interaction between the parameters. For example, will be the contribution of the uncertainty in the parameter p1 to the uncertainty in the output due to the functional dependence of the output on the parameter p1, and will be the contribution of the uncertainty in the parameters p1 and p2 to the uncertainty in the output due to the functional dependence of the output on the parameters p1 and p2 as captured by . In this case, the contribution of is referred to as the main effect or the first-order effect of the parameter p1, the contribution of is referred to as the interaction effect of the parameters p1 and p2.

The Sobol’ indices are then defined as,

Intuitively, the ranking of will induce a ranking of the contribution of the corresponding Sobol’ function f to the variance of the output Y. We denote the random variables corresponding to the output y and a parameter p, by Y and P, respectively. Then, to consider the total contribution of the variance in P to the variance of Y, contributions from all of the Sobol’ functions of p should be considered. Following this line of reasoning, the total Sobol’ index for a parameter P, is defined by,

It follows from Eq. (18) that, if , then f = 0 for all I ⊂ D where i ⊂ I. This shows that the QoI does not functionally depend upon the input p in the chosen parameter space, , and we can conclude that the presumed variation in the parameter p is of no consequence to the QoI. Hence, in the sensitivity analysis using SIs, a parameter is deemed as unimportant if the corresponding total SI is negligible. Consequently, the choice of the value of the parameter p, in the analyzed parameter space, will not have any effect on the QoI.

Furthermore, a statistical interpretation for the total Sobol’ index can be seen from where P~ is the vector of random variables for all parameters except p [39]. V(Y|P~) is the conditional variance in the output if we vary only the parameter p and fix all of the remaining parameters, giving a measure of local sensitivity of the QoI to p. For a global sensitivity measure, variation in all of the parameters has to be considered simultaneously, and as a result, we can consider the statistical distribution of V(Y|P~) over the N — 1 dimensional space corresponding to P~. then measures the expected value of this distribution, i.e. the expected value of the conditional variation in the output with respect to p over the space of all of the remaining parameters (E [V(Y|P~)]), normalized by the variance in the output (V(Y)). A high therefore indicates that we expect the variance in the output to be high with respect to the variance in p, considering simultaneous variation in all of the remaining parameters as well, and we classify it as a significant parameter. On the other hand, a negligible will indicate that we do not expect any variance in the output with respect to p and it can be frozen at any value in the [a] without significantly affecting the QoI over the entire parameter space Ω.

The first-order or main Sobol’ indices are defined by,

We define the additional SI of a parameter p as the contributions of all of the interaction effects involving the parameter p, and denote it by . Hence,

Further analysis reveals the following relations between main Sobol’ indices and statistics of the QoI-distribution [39],

It follows from Eq. (22), that high indicates a high expectation of a reduction in the variance of the QoI if the value of P is known.

Numerical estimation of SIs

Since the integrals required to define SIs might not always be determinable analytically, numerical methods are often used to estimate these indices. The available methods can be classified into three categories: hybrid pick-and-freeze, extended Fourier amplitude sensitivity test (eFAST), and generalized polynomial chaos expansion (gPC). The concept for hybrid pick-and-freeze was provided by Sobol’ in his seminal paper for SIs [30], in which he used Monte-Carlo methods to estimate the integrals corresponding to the SIs. The Monte-Carlo scheme has since been improved, and reviews and comparisons of different hybrid pick-and-freeze strategies can be found in [39, 41].

In this work, we implement the hybrid pick-and-freeze of Jansen ([42]) to estimate total SIs. Jansen’s scheme has been shown to be more efficient than similar sampling-based schemes, especially when used in conjunction with quasi-random sequences [39, 41]. To estimate the first-order SI, we implemented the “Correlation 2” scheme proposed in [41]. This scheme uses a larger number of samples than many other first-order SI estimators, but we prefer this scheme because it has been shown to provide a better estimation in the case when the SIs for some of the parameters may be small. For more information on the estimation of Sobol’ indices, see Appendix D.

Results

As previously stated, our Quantity of Interest will be the steady-state pH at the left boundary (epithelial wall)

Intuitively, a large SI for a given parameter means that fluctuations in that quantity may cause large fluctuations in the pH experienced by the mucosal tissue of the stomach and potentially be detrimental to stomach health. Conversely, a small SI indicates that variations in that parameter do not appreciably contribute to variations in the pH at the stomach wall, which could be physiologically advantageous.

Full Sobol’ Analysis

One of the central results from [18] was the quantification of the mucosal pH as a function of kHI and kAB. In that study, all other parameters were fixed using values estimated from biological data. It was shown that the mucosal pH only varied appreciably in the regime of low kHI and/or low kAB. Here, we extend that study by including simultaneous variation in other model parameters (due to different experimental conditions, or simply natural variations between people) and use Sobol’ indices for sensitivity analysis. Table 1 shows all of the parameters of interest, their physical units, and the intervals on which we sample them. The final column also indicates which parameters were assumed to be log-uniformly distributed (as opposed to uniformly distributed).

Table 1:

This table lists the relevant parameters for our Sobol’ analysis, their units, and the intervals over which we allow them to range. The final column indicates if the parameter is assumed to be log-uniformly distributed (otherwise, parameters are assumed to be uniformly distributed).

Parameter	Description	Units	Interval	Log-unif.
k HI	Rate constant of Hydrogen/Sodium Exchange	cm/sec	[1 × 10−6, 1 × 10−1]	✓
k AB	Rate constant of Chloride/Bicarbonate Exchange	cm/sec	[1 × 10−6, 1 × 10−1]	✓
δ HI	Bias constant of Hydrogen/Sodium Exchange	[−]	[5 × 104, 15 × 104]
δ AB	Bias constant of Chloride/Bicarbonate Exchange	[−]	[5, 15]
HL	Lumenal concentration of Hydrogen	M	[1 × 10−4, 1 × 10−1]	✓
IL	Lumenal concentration of cations	M	[1 × 10−4, 1 × 10−1]	✓
S 0	Magnitude of Hydrogen secretion	[−]	[1 × 10−1, 1 × 101]	✓

We begin with varying all of the parameters except for S0, which is fixed at its nominal value of 1. This set of parameters captures those that were either varied, or roughly estimated in [18], and thus represents an approximate reproduction of the analysis from that work using the SI methodology. The estimated SIs are shown in Fig. 2a. The blue bars indicate estimated first-order SI (), while the orange bars indicate the additional SI (). The combined height of the orange and blue bars gives the total SI for each variable (). It is immediately clear that the two parameters with the largest total SI are the ion exchange rates (kHI and kAB and the bias constants δHI and δAB) have negligible total and first-order SIs. Furthermore, any variations in the lumenal concentrations (HL and IL account for strictly less than 20% of the variance in the QoI in this parameter regime, with their combined first-order effects accounting for less than 7% of the variance in the QoI. Therefore, the Sobol’ analysis shows that variations in kHI and kAB are necessary and sufficient to bring about a significant variation in the QoI, without any significant effects of variations in the other parameters. In this regard, Sobol’ analysis supports and extends the work in [18]. Together with the previous study, our analysis implies a remarkably robust mechanism for gastric pH maintenance.

Fig. 2:

Sobol’ indices for boundary pH for each parameter. First order () and additional () SI are depicted, while their sum indicates total SI (). (a) Source is fixed at 1, and all other parameters are varied in their region of interest, (b) All seven parameters are varied.

Next, we also include the effects of variation in the magnitude of the Hydrogen secretion, S0. The estimated SIs are shown Fig. 2b. Based on the convergence criterion that we define in Appendix D, 20000 evaluations of the QoI were required to estimate the total SIs and 48000 evaluations were required to estimate the first-order SIs. The mean of the evaluations used to estimate total SI is around 6.235 with a variance of 0.789. The comparison between total SIs of parameters again shows that the two rate constants are the most significant. The remaining parameters fall roughly into three categories. The bias constants δHI and δAB still have negligible SI (< 0.05). The lumenal boundary concentrations H and I have small SI (~ 0.1). However, the ion source magnitude S0 has significant SI (~ 0.2). This means that variation in S0 cannot be entirely ignored since it has non-negligible total SI along with a significant first-order effect – varying S0 by itself can account for 15% of the variance in the QoI.

To more precisely describe the impact each parameter has on the model predictions, we use violin/box-and-whisker plots of the QoI as each parameter is varied (see Fig. 3). Figure 3a illustrates the process of generating these plots for the parameter S0. The grey points indicate a simple scatter of the measured QoI as a function of the parameter S0. The dashed and solid purple lines indicate a running mean and mean ± standard deviation of the QoI. However, this format makes it difficult to visually infer trends in the distribution of the QoI as S0 is varied. Therefore, we divide the range of S0 into ten equal width sub-intervals and generate a violin for each. The white hash indicates the mean of the QoI for samples in that respective subinterval and the black hash indicates the median. The thicker black line indicates the region between the first and third quartiles. The black whiskers indicate the extent of the data (sans outliers), and the blue x-es indicate individual outliers defined as median ±1.5 times the inter-quartile range. Finally, the width of the red violins represent a smoothed histogram of the data within each subinterval. For visual clarity we present only the violin plots, and not the underlying scatter plots, for all remaining parameters.

Fig. 3:

Violin and box plots for log10 (H0 as a function of (a) ion source magnitude (S0), (b) hydrogen exchange rate (kHI), (c) bicarbonate exchange rate (kAB), (d) lumenal hydrogen concentration (HL), and (e) lumenal cation concentration (IL). White hashes indicate mean of QoI within a subinterval and black hashes indicate the median. Thick black lines indicate the range from first to third quartiles. Black whiskers indicate the extent of the data (sans outliers), and blue x-es indicate individual outliers (median ±1.5 times inter-quartile range). Panel (a) also shows the underlying data (with running mean and mean±std. indicated) for illustrative purposes.

Visual inspection of Fig. 3a shows that the mean value of wall pH changes by approximately 1 as S0 is varied over the given interval. From Eq. (23), we can interpret the ratio of this variation in the expected value of QoI to the total variation of QoI as the the first-order SI of the parameter under consideration (). Other changes in the distribution of the QoI (from violin to violin) are attributable to additional SI of that parameter (). For example, the distribution of QoI is bimodal when S0 is large but unimodal when S0 is small.

The observed change in mean QoI also has a very natural and unsurprising physiological interpretation: when the source of hydrogen due to gastric secretion is larger, the gastric wall (on average) is more acidic. However, we note here that variations in secretion rate (S0) over two orders of magnitude produces a variation in mean gastric pH of less than one. Thus, even though S0 exhibits a significant total SI, the mean mucosal pH is well within biologically acceptable limits (≈ 6 or 7). This seems physiologically advantageous, as one would expect in-vivo gastric secretion to vary throughout the day in response to various nervous and endocrine stimuli. In Section 4.2, we investigate more closely the behavior of mucosal pH in the limit of large S0 (i.e. secretion in response to a recent meal or histamine stimulus) and small S0 (i.e. secretion in response to weak/no stimulus).

Figures 3b and 3c show violin plots for QoI as a function of kHI and kAB. As expected, because these two parameters have the largest first order SI, we observe a notable variation in the mean of the QoI. As either exchange rate parameter is decreased, the measured pH decreases (on average), implying acidification of the gastric mucsoal wall and suggesting a failure to maintain homeostasis. This is in line with the results presented in [18], which indicated that as long as both exchange rates were sufficiently large (i.e. sufficiently rapid exchange of ions between the epithelium and gastric lumen), a healthy pH was robustly maintained. The present analysis also allows us to determine how the variance in the predicted wall pH changes as a function of ion exchange rates. In particular, we can see that when kHI is low, there is a significant increase in the variation of the QoI. As this behavior affects higher order moments of the QoI distribution, it may be related to the additional affects within the total SI of kHI. A similar (but distinct) result was presented in [18], where it was shown that mucosal pH could vary wildly when kHI = 0. In Section 4.3 we investigate more closely the behavior of mucosal pH in the limit of small kHI. Finally, we note that the parameter kAB does not obviously affect the spread of the distribution of QoI, and therefore exhibits less additional SI ().

The parameters H and I both have relatively small SIs. These parameters represent the concentration of ions in the interior of the stomach lumen, which could vary in the in vivo system of interest (in response to diet, for example). Small SI imply that the pH experienced by the stomach wall does not vary appreciably in response to these temporal fluctuations, implying a robust physiological control mechanism. In this regard, the data presented here represents a more quantitative (and global) generalization of the argument put forward in [18]. We note however, that these SIs are not completely negligible and are comprised of significant additional SI. For this reason we include the violin plots of both in Figs. 3d and 3e. Of particular note is that even though the mean QoI does not depend strongly on I ( is small), when I is large the variation in QoI is drastically decreased. In this regard, one could describe the model as extremely robust in the limit that I is large: regardless of other parameter values, the pH at the stomach wall will almost certainly be ≈ 6.

We do not include violin plots of δHI and δAB in this section, as both of those parameters have negligible SI. However, those graphs can be found in Appendix A for completeness. We also note here that there is no “accepted” value for these parameters in the literature. In [18] their values were roughly estimated based on data available from other epithelial cell types. The negligible SI for δHI and δAB imply that the results presented in [18] do not depend on the these estimated values and are generically representative of the model behavior.

Analysis for Fixed Source Magnitude

The results of Section 4.1 indicate that the ion exchange rates kHI and kAB exhibit the largest SI, which when interpreted in conjuction with [18], imply a robust physiological control of mucosal pH. However, they also indicate that the hydrogen source magnitude (or secretion rate) S0 may make a significant contribution to the sensitivity of the model predictions, which was not addressed in previous investigations. This is physiologically relevant, as the secretion rate of hydrogen in the human stomach is known to vary by a factor of 40 or more, in response to histamine levels or the recency of a meal [14, 22]. We analyze the sensitivity of the model predictions when secretion is at its extremal values, i.e. when the stomach is quiescent and not secreting much gastric juice, or after intense stimulation. In this section, we perform a second analysis where we fix the value of S0 and calculate the SI (primary and additional) for the remaining six parameters. This is done twice: with S0 = 0.1 and with S0 = 10.

Large S0:

Upon fixing S0 = 10 (which represents stimulated gastric secretion), we recalculate the SI of the other six parameters. This requires approximately 18000 evaluations of the QoI. The primary and additional SI are shown in Fig. 4. In this regime, kHI and kAB remain the dominant SI, and are predominantly characterized by first order SI. Remarkably, no other parameters exhibit significant SI. This can be interpreted to mean that in the regime of strong gastric secretion, the ion exchange rates nearly completely govern the mucosal pH.

Fig. 4:

Sobol’ indices for the other six parameters when hydrogen source is fixed and large (S0 = 10). First order () and additional () SI are depicted, while their sum indicates total SI ().

Figure 5 shows the violin plots for QoI as a function of kHI and kAB. Consistent with the SI shown in Fig. 4, when ion exchange rates are small the gastric wall is likely to be significantly more acidic. We note that there is a marked increase in the variation of QoI when kHI is small and this is the only parameter regime when the wall pH is strongly acidic (< 4). This result is also consistent with those presented in [18]: as long as ion exchange rates are “fast enough”, the mucosal wall is generically maintained at near-neutral pH. The violin plots for the remaining four parameters do not add much to our discussion here (and their SI are insignificant), but we include them in Appendix A for completeness.

Fig. 5:

Violin and box plots for log10 (H0) as a function of (a) hydrogen exchange rate (kHI) and (b) bicarbonate exchange rate (kAB), when hydrogen source is large (S0 = 10). White hashes indicate mean of Qol within a subinterval and black hashes indicate the median. Thick black lines indicate the range from first to third quartiles. Black whiskers indicate the extent of the data (sans outliers), and blue x-es indicate individual outliers (median ±1.5 times inter-quartile range).

Small S0:

After fixing S0 = 0.1 (which represents suppressed secretion of gastric juice), we again calculate the SI of the remaining parameters. Convergence required roughly 20000 evaluations of the QoI, and the calculated SI are shown in Fig. 6. The ion exchange rates still exhibit the largest SI; however in this regime we also observe significant SI for the lumenal ion concentrations H and I. This implies that the mucosal pH is more sensitive to the ionic composition of the gastric interior when secretion is suppressed. Furthermore, all four parameters exhibit appreciable additional SI, implying that there are significant interaction effects between parameters. As in all other analyses, the SI of both offset parameters is negligible.

Fig. 6:

Sobol’ indices for the other six parameters when hydrogen source is fixed and small (S0 = 0.1). First order () and additional () SI are depicted, while their sum indicates total SI ().

Figures 7a and 7b show the violin plots for kHI and kAB. As in all previous analyses, we observe that increased ion exchange generically leads to a less acidic mucosal wall. However, the mean QoI exhibits a relatively weak dependence on either exchange rate parameter. This is consistent with the relatively small first order SI shown in Fig. 6. We also note here that the mucosal wall is generically neutral (or near neutral), always with a mean pH ≳ 5.5. This is physiologically unsurprising, as these results are produced in the regime where less hydrogen is being secreted into the gastric lumen.

Fig. 7:

Violin and box plots for log10 (H0) as a function of (a) hydrogen exchange rate (kHI), (b) bicarbonate exchange rate (kAB, (c) lumenal hydrogen concentration (H), and (d) lumenal cation concentration (I), when hydrogen source is small (S0 = 0.1). White hashes indicate mean of QoI within a subinterval and black hashes indicate the median. Thick black lines indicate the range from first to third quartiles. Black whiskers indicate the extent of the data (sans outliers), and blue x-es indicate individual outliers (median ±1.5 times inter-quartile range).

Figures 7c and 7d show the violin plots for H and I. Figure 7c shows that an increase in H (which represents a more acidic gastric lumen) likely leads to a more acidic mucosal wall. However, it is also associated with an increase in the variance of the QoI. This means that when secretion is low, the mucosal pH is more susceptible to variations in lumenal acidity but may be able to maintain neutral (or near neutral) conditions depending on other factors. Conversely, Fig. 7d shows that the QoI exhibits a drastic decrease in variance associated with high values of I. This effect is even more pronounced than that observed in Section 4.1. The analysis shows that when gastric secretion is low, but sodium concentration in the stomach lumen is high, the mucosal pH will nearly always be ≈ 6.

Sobol’ Analysis for low kHI

A recurring theme in the above sections is that kHI exhibits a greater SI than the other parameters considered. Physically, this is somewhat surprising, as the analysis in [18] showed that under most conditions, the flux of hydrogen ions through hydrogen/sodium exchange is significantly smaller than the other sources/sinks of hydrogen in the system (secretion and bicarbonate buffering). However, the analysis of [18] also concluded that kHI is a critical parameter of the system: when it is identically zero (representing a lack of hydrogen/sodium exchange proteins), the system fails to robustly maintain a near-neutral mucosal pH. Our analysis here has corroborated this finding; in all cases we see that low values of kHI lead to drastically increased variance in QoI.

To further explore the behavior of the system in the limit of “vanishing” hydrogen/sodium exchange, we perform a series of calculations where kHI is fixed at a low values, and the SI of the remaining six parameter are calculated. We calculate the SI of the remaining parameters for five values of kHI ranging from 10−6 to 10−8. Calculating these SI required roughly 24000 to 26000 evaluations of the QoI, depending on the value of kHI. We do not explore the case when kHI = 0, as this caused the linear operators associated with our numerical scheme (see Appendices B and C) to become ill-conditioned.

Figure 8a shows the SI of the remaining six parameters when kHI = 10−6; Fig. 8b shows the same SI when kHI = 10−8. In both Figs. 8a and 8b, we see that though the numerical values have changed, the SI of the remaining parameters are qualitatively similar to the full Sobol’ analysis: kAB exhibits the dominant SI, followed by S0, I, and H with moderate SI, and finally the ion exchange offsets (δHI and δAB) have negligible SI.

Fig. 8:

First order and additional Sobol’ indices for the other six parameters when hydrogen exchange rate is fixed and small. Panels show (a) kHI = 10−6 and (b) kHI = 10−8. First order () and additional () SI are depicted, while their sum indicates total SI ().

Comparing Fig. 8a to Fig. 8b, we observe two trends. First, we note that decreasing kHI leads to significant increase in total SI for I and a moderate increase in the SI of H. This can be interpreted to mean that when there is a deficiency in (or shortage of) hydrogen/sodium exchange proteins in the epithelial wall, mucosal pH becomes much more sensitive to the ionic composition of the stomach lumen. To more clearly illustrate this phenomenon, Fig. 9a shows the total SI for the remaining size parameters as a function of kHI. It clearly illustrates the increase of as kHI → 0. One can also observe moderate increases in and , as well as a moderate decrease in . Generically, the SI for ion exchange offset parameters is negligible. These analyses all corroborate and extend the work presented in [18].

Fig. 9:

(a) Total Sobol’ indices of the other six parameters as a function of kHI. (b) Mean (blue diamonds) and variance (orange squares) of pH at the left boundary as a function of kHI.

The second trend of note in Fig. 8 is that decreasing kHI appears to decrease the first order SI and increase the additional SI for all parameters. This indicates that are significant interaction effects and may suggest variation in higher order moments of the QoI as parameter values change. For this reason we plot the mean and variance (over all evaluations used to calculate the SI) of the mucosal pH as a function of kHI. The results are shown in Fig. 9b. As can be seen, the mucosal pH becomes (on average) more acidic as kHI decreases. However, there is also a dramatic (> 3-fold) increase in the variance of mucosal pH; the coefficient of variation for the mucosal pH increases from approximately 0.23 (when kHI = 10−6) to approximately 0.62 (when kHI = 10−8). Taken together, these facts imply that as hydrogen/sodium exchange is down-regulated, the bicarbonate buffering mechanism becomes unable to robustly maintain a neutral environment adjacent to the mucosal wall.

Discussion

In this paper, we have presented a global sensitivity analysis of steady state epithelial pH predicted by an electrodiffusive model of hydrogen transport in the gastric mucus layer. Using Sobol’ Indices, we have provided a systematic extension of the (relatively local) sensitivity analysis of this model reported in previous work. Our results support the existing analysis of [18], while also providing further physiological insight by considering variation in biological parameters which were considered “fixed” in previous works.

The analysis presented here repeatedly shows that the rate constants of ion exchange through the epithelial surface (kHI and kAB) are the dominant parameters that impact the maintenance of the pH at the epithelial wall. These two parameters exhibited the largest sensitivities in all regimes considered, particularly when the source of hydrogen secretion is large (i.e. stimulated gastric secretion). In the regime of suppressed gastric acid secretion, the parameters representing the concentration of ions in the gastric lumen exhibited significant sensitivities. Taken together, these results indicate that when the stomach is not stimulated (and not producing large amounts of gastric juice), the contents of the stomach may impact the pH of the epithelium. However, when hydrochloric acid is actively introduced to the stomach lumen (i.e. secretion is stimulated by ingesting a meal, or treatment with histamine) the exchange of ions at the epithelial wall becomes the dominant process and exerts significant control over wall pH.

Furthermore, the large sensitivities of both exchange rate parameters (kHI and kAB) suggest that both hydrogen/sodium and bicarbonate/chloride exchange events are critical to the control of the hydrogen concentration adjacent to the gastric epithelium. This conclusion is in line with the work of [18], but to our knowledge has not been suggested anywhere in the gastric physiology literature, where the physiological purpose of the hydrogen/sodium exchange proteins is rarely, if ever, discussed.

The observation that epithelial pH is sensitive to hydrogen/sodium exchange even though this relationship is not appreciated in the physiology literature, together with the analysis of [18], inspired us to perform a sequence of sensitivity analyses where the SI of the remaining parameters were calculated for several (small and decreasing) values of kHI. Notably, the sensitivity measure of lumenal sodium concentration increases. This is perhaps unsurprising, as the electroneutrality constraint implies that the concentration of sodium and hydrogen at the epithelial wall are strongly dependent on one another, and impeding the regulatory function of hydrogen/sodium exchange at the epithelium would reasonably leave epithelial pH more sensitive to sodium diffusing from the lumen. Perhaps more surprising is that the sensitivity measure of hydrogen source magnitude decreases as kHI decreases. As the rate of hydrogen/sodium exchange decreases, one would expect the epithelium to become more sensitive to the amount of hydrogen secreted by the stomach. This illustrates a common misunderstanding in how SI are often interpreted in the broader community. Directly comparing SI values generated with two different parameter sets (kHI = 10−6 cm/sec vs. kHI = 10−6 cm/sec) may be misleading, as these values are expressed as a proportion of the total variation in the QoI, which may be changing (see below). What can be concluded is that as kHI decreases, the QoI becomes relatively more sensitive to I than S0, as when kHI 10−6, but when kHI = 10−8.

Finally and perhaps most importantly, our analysis shows that as the rate of hydrogen/sodium exchange decreases, the mean pH of the gastric epithelium also decreases, while the variance of this same quantity increases nearly five-fold. Both of these changes have drastic implications for healthy physiological function. A decrease in pH is generally detrimental to the health of the cells that make up the epithelial surface and is associated with numerous pathologies. Indeed, it is generally accepted that a major role of the gastric mucus layer is to protect the epithelium from the low pH of the lumen [11]. Therefore, any scenario which leads to a decrease in the lumenal pH would be physiologicaly disadvantageous. Furthermore, the increase in the variance as kHI decreases implies that the epithelial pH is generally more sensitive to perturbations in the other parameters, and thus one could expect larger variations in wall pH. Taken together, these results imply something of a “double whammy” for the maintenance of healthy gastric pH: as the rate of hydrogen/sodium exchange decreases towards zero, the system is not only more likely to maintain an “unhealthy” pH at the wall, but it is more likely to experience large swings in said pH.

In closing, we would like to note that several of these results represent a powerful use for Sobol’ Indices (and sensitivity analysis in general) which remains underappreciated in the literature. Sensitivity analysis is often used to simply quantify uncertainty in model predictions associated with uncertainty in parameter estimation, and no deeper implications are considered. However, there are specific biological insights that can be determined by judicious use of sensitivity analysis – this includes varying the QoI, parameter regimes, and exploring the QoI distributions. Furthermore, the physical and/or biological quantities that model parameters represent are often not fixed, but rather dynamic quantities subject to random or deterministic fluctuations. Sensitivity analysis can be used to make quantitative inferences about the robustness of the system of interest to these fluctuations, and such inferences can have profound implications for how biological systems function in the face of a dynamic environment.