Electrodiffusion-Mediated Swelling of a Two-Phase Gel Model of Gastric Mucus.

Gastric mucus gel is known to exhibit dramatic and unique swelling behaviors in response to the ionic composition of the hydrating solution. This swelling behavior is important in the maintenance of the mucus layer lining the stomach wall, as the layer is constantly digested by enzymes in the lumen, and must be replenished by new mucus that swells as it is secreted from the gastric wall. One hypothesis suggests that the condensed state of mucus at secretion is maintained by transient bonds with calcium that form crosslinks. These crosslinks are lost as monovalent cations from the environment displace divalent crosslinkers, leading to a dramatic change in the energy of the gel and inducing the swelling behavior. Previous modeling work has characterized the equilibrium behavior of polyelectrolyte gels that respond to calcium crosslinking. Here, we present an investigation of the dynamic swelling behavior of a polyelectrolytic gel model of mucus. In particular, we quantified the rate at which a globule of initially crosslinked gel swells when exposed to an ionic bath. The dependence of this swelling rate on several parameters was characterized. We observed that swelling rate has a non-monotone dependence on the molarity of the bath solution, with moderate concentrations of available sodium inducing the fastest swelling.

1. Introduction

Gastric mucus is a polyelectrolyte gel which serves an important biological function. It is generally accepted that the gastric mucus layer provides a protective barrier between the stomach wall and the stomach interior, shielding the wall from acid and digestive enzymes and preventing auto-digestion of the stomach epithelium. However, there are several outstanding questions regarding the maintenance of the gastric barrier. Digestive enzymes within the stomach are constantly degrading the mucus layer from the inside-out, cleaving the mucin glycoproteins and destroying the entangled network which makes up the mucus gel [1]. This degradation of the network must be balanced by secretion of fresh mucus from the stomach wall to maintain a stable layer.

Fresh mucus is secreted from epithelial cells on the stomach wall via exocytosis of microscopic vesicles filled with densely packed mucin polymers. This dense packing of negatively charged polymers is shielded by an abundance of divalent calcium ions within the vesicle [2]. When the contents of these vesicles are exposed to an aqueous environment, the mucus swells dramatically [3]. However, this swelling behavior is known to depend on environmental pH and ionic composition [4]. In particular, when dense mucus is secreted into a deionized environment, the dramatic swelling event does not occur [3]. Experiments have shown that swelling is accompanied by a massive (and rapid) transport of monovalent cations (such as sodium or hydrogen) into the densely packed mucus in exchange for the divalent calcium [5], removing the efficient charge shielding. It is possible that the divalent nature of calcium ions provides an explanation for the observed behavior. Due to having a charge of , calcium can form two bonds with the negatively charged mucin polymers, effectively forming a “crosslink”. Even though these crosslinks are much less stable than (for example) one potentially formed by covalent bonds [6,7,8,9,10], large numbers of crosslinks may allow for very dense network configurations. Conversely, replacing calcium with monovalent sodium removes the ability to form crosslinks, potentially causing rapid hydration of the gel.

The dynamics of swelling hydrogels have been studied for over seventy years, using theoretical frameworks of increasing sophistication [11]. However, the proposed mechanism must necessarily depend on electrodiffusive transport of ionic species, motion of the glycoprotein network and hydrating fluid, and chemical interactions between the network and dissolved ions. All of these processes are coupled and affect one another. In [12], the authors derived from first principles a model describing the transport of mono- and divalent ionic species which bind and unbind to a gel-like material described by a two-phase mixture model. The equilibrium behavior of the model is explored in [12], and the linear behavior near equilibria is explored in [13]. We present here the first general investigation of the dynamic swelling behavior of this model. In particular, we attempted to quantify the rate of swelling of a polyelectrolyte gel material as a function of the parameters which govern its chemical interaction with dissolved ions.

2. Mathematical Model

In this work, we choose to model a mucus gel using a two-phase framework. Similar frameworks have been used in a wide variety of biological contexts, including the modeling of blood clots, cellular mechanics, and bacterial biofilms [14,15,16,17]. In our framework, the local composition of the complex gel material is described by the volume fractions of two distinct phases which we call network (denoted with subscript n) and solvent (subscript s). The network represents the entangled mesh of polymeric proteins which give rise to the gel, while the solvent represents the interstitial hydrating fluid. Each material is allowed to move with its own distinct velocity. Dissolved ionic species are allowed to move through the solvent phase, and bind/unbind to/from the network phase.

2.1. Gel Evolution

As is standard in a two-phase model, there are four equations for the composite material composed of network and solvent. The first two equations describe conservation of mass of the solvent and network phases, respectively: Here, and are the volume fractions of network and solvent, respectively, and and are their respective velocities. By definition, the volume fractions satisfy the constraint . As a consequence of this constraint, Equations (1) and (2) imply that the volume-averaged incompressibility condition must hold. The velocity fields of each phase are determined by conservation of momentum. Because the spatial scale of interest is small (microns) and the velocities of the solvent and network are low, it is reasonable to ignore inertial forces. Hence, momentum balance reduces to force balance equations for each phase, Here, and are stress tensors which encapsulate the internal stresses within the network and solvent, respectively; is the coefficient of the drag force that arises whenever there is relative motion between the two materials; , , and for are chemical potentials that act on the solvent, network, and ionic species, respectively; and p is the hydrodynamic pressure. The final term in Equation (5) replaces the drag forces between the ions and the solvent. In effect, the forces from the ionic potentials are transferred to the solvent because the only forces that act on the ions of type j are that from and the drag between the solvent and these ions, and these forces must sum to zero. The pressure p is determined by the incompressibility constraint in Equation (3). In Equation (5), is the ratio of the number density of ion j to the number density of water particles. This ionic contribution is valid in the limit that ions are dilute in the solvent.

The internal stresses within each phase must be given by a constitutive law. In this work, we assume that both the solvent and network exhibit a viscous response (though the modeling framework is adaptable to account for elastic or viscoelastic constitutive laws as well). In one spatial dimension, the stress tensor for each phase reduces to where is the viscosity of the respective phase.

2.2. Dissolved Ion Evolution

The ionic species dissolved within the solvent phase are each subject to a Nernst–Planck type equation, modified to account for diffusion through a fluid of non-uniform volume fraction. Cations may bind to the network via mass-action type reactions, while anions may not. We are specifically interested in the behavior of the model when we have three distinct types of ions present: a monovalent cation (which we call sodium), a divalent cation (calcium), and a monovalent anion (chloride). Here, is the concentration of species j measured in units of moles per liter solvent volume, is that ion’s diffusion coefficient, its valence, and is the electric potential measured in terms of thermal voltage (R: ideal gas constant; T: absolute temperature; and F: Faraday constant). The parameter is the binding rate of ion j to the gel network, while is the corresponding unbinding rate. The ratio is the dissociation constant of that cation.

2.3. Bound Ion Evolution

Cations that are bound to the mucus network advect with the network velocity. The concentration of bound ions may also change through binding and unbinding reactions described by Here, is the concentration of sodium bound to the network, measured in moles per total volume. Because calcium is divalent, it may bind twice, and therefore exists in one of two bound states: singly bound () and doubly bound (). The quantity M is the concentration of negative charge on the mucus network which is not bound to a cation; it is determined by the relationship where is the density (per network volume) of negative sites on a unit of network. We are assuming that negative charge on the network and binding sites which cations may occupy exist in a one-to-one ratio.

2.4. Driving Potentials

The potentials which act on the volume-occupying species (network and solvent) account for entropic, long-range electrostatic, and short-range internal energy effects. The potentials acting directly on the solvent are given by and the potentials acting directly on the network are given by Terms in Equations (14) and (15) labeled (I) represent entropic effects, and capture osmotic pressure acting on the solvent. Terms labeled (II) represent potentials due to short range interactions and generalize the standard Flory–Huggins mixture theory. We discuss the factor in more detail below. Finally, the term labeled (III) represents the potential due to the charged nature of the mucus network. Here, is the Boltzmann constant and N is the number of monomers in a typical polymer chain of mucin. The term is the total ionic molality, given by where, again, is the concentration of ion j divided by the standard molarity of water ( M).

Because the ionic species are regarded as massless, the forces on each species are in balance. Thus, the chemical potential force and the solvent drag force on each ionic species sum to zero. Since the drag force the ions exert on the solvent is the opposite of that the solvent exerts on the ions, the net effect is that the chemical potential forces appear to act on the solvent itself (as incorporated in the last term of Equation (5)). The potentials which act on dissolved ions are given by where is the valence of ion j. Equation (17) accounts for both entropic (I) and electric (II) potentials.

We now return to the term , which appears in Term (II) in Equations (14) and (15). We refer to this as the “interaction parameter”, and note that it is somewhat analogous to a typical Flory interaction parameter [11]. However, here is not a constant; it is given by where and are constants described below. The quantity is the ratio of the concentration of current crosslinks within the network () to the maximum concentration of crosslinks possible (): Because may vary spatially and temporally as the concentration of doubly bound calcium and network volume fraction do, so may . Thus, our modeling framework is capable of describing a gel network with a “Flory interaction parameter” that varies spatiotemporally in response to the local ionic concentrations.

Finally, terms and (in Equations (14) and (15)) are the so-called “standard free energy” of the mucus phase and solvent phase, respectively. The solvent term is a constant, and therefore does not affect the behavior of the system (as only gradients of potentials appear in the force balance). The network standard free energy is given by Note that the crosslinking fraction appears in this expression. The parameters , , , and are referred to as the interaction energies, and arise in a standard mean-field calculation of the internal energy associated with a given mucus/solvent mixture. A detailed derivation of this modeling framework, especially the potentials given in Equations (14), (15) and (17), may be found in [12].

3. Results

3.1. Numerical Experiments and Initial Conditions

Our goal was to investigate the effect of the ambient ionic concentration on the dynamic swelling of high density network. To this end, we constructed initial profiles of volume fraction and ionic concentrations which represent a dense, highly crosslinked amount of network in a small region immersed in a fluid of known ionic composition. We refer to the region with appreciable volume fractions of network as the “inner” region, and initially placed it at the left end of the one-dimensional domain of width µm. The rest of the domain is referred to as the “bath” and initially has essentially no network and a prescribed ionic composition. To construct the initial profiles, we prescribed , and the total amount of calcium, sodium, and chloride. Under the assumption that all species are spatially uniform, Equations (7) to (12) may be solved to determine the individual concentrations ( and ) at which the binding and unbinding chemistry is at local equilibrium. This process was carried out separately for the inner region and the bath. We then used a hyperbolic tangent function (centered at µm) to construct spatial profiles that transition smoothly but sharply between the inner and bath concentrations. These profiles were used as the initial conditions for Equations (7) to (12). An example of the initial conditions can be seen in Figure 1. We simulated Equations (7) to (12) using no-flux boundary conditions for all species and zero Dirichlet conditions for both velocities (to represent a closed container) to investigate the dynamic swelling behavior of the network. Please see Appendix A.1 for details of the numerical techniques used.

Figure 1

Example initial profiles of: volume fraction (left); bound ionic concentrations (middle); and dissolved ionic concentrations (right).

In all simulations, the initial volume fraction of network, the mechanical and energetic properties of the network (the parameters , , , , , , and N), and the diffusion coefficients () of the ions were identical. The values of these parameters are listed in Table 1. In particular, the were chosen to produce a network which has a very weak propensity to swell when fully cross-linked, but swells quickly when no calcium crosslinking is present. When fully crosslinked (), the network exhibits an interaction parameter , which does not drive swelling at all (this is analogous to a Flory interaction parameter of zero). However, when there is no crosslinking (), the network exhibits an interaction parameter , which is analogous to a large, negative Flory interaction parameter, and leads to rapid swelling. The parameter was chosen to be small and negative, as large negative values of this parameter can lead to gels which tend to phase separate, which is not a behavior we wanted to investigate. The final interaction energy does not impact the system dynamics, as it is simply an additive constant that produces no potential gradient. In all simulations, the bath concentration of sodium was large, while the calcium concentration was small. The chloride concentration was chosen to maintain electroneutrality. This choice was made to allow monovalent sodium from the bath to displace calcium crosslinks within the network.

Table 1

Parameters that do not vary between sets of experiments.

Parameter	Value
ϵ_1	−45
ϵ_2	25
ϵ_3	−0.5
ϵ_4	0
N	6
inner θ_n	0.5
bath θ_n	1×10^−9
Total bath Ca	1 mM
Total bath Na	0.2–200 mM
Total bath Cl	0.2–202 mM
D_Na	2.5×10^−5 cm^2/s
D_Ca	2.5×10^−5 cm^2/s
D_Cl	2.5×10^−5 cm^2/s

We performed four major sets of experiments to investigate how two basic aspects of the network chemistry affect swelling behavior. Specifically, we investigated the regime where there is an extremely high density of binding sites on the network (referred to as “dense binding”), as well as the regime where this is a more moderate density of binding sites (referred to as “sparse binding”). Estimates for the density of cation binding sites on salivary mucus have been reported to be less than one molar [18]. However, the experiments in [19] suggest that gastric mucus may bind hydrogen ions in concentrations above one molar. Therefore, we assumed that a “sparse” network exhibits 0.1 mole of binding sites per liter of pure mucus () while a “dense” network has 1 mole of binding sites per liter of pure mucus (). We also investigated the case where a network–calcium bond is chemically preferred to a network–sodium bond (), as well as the case where network–sodium bonds are preferred (). We refer to these cases as “calcium preferred” and “sodium preferred”, respectively. We investigated both cases, as estimates for the dissociation constants for both sodium and calcium vary over orders of magnitude depending on the experimental setup and type of mucus used [18,20]. Furthermore, we know of no reliable estimates for the individual kinetic rate constants for either ion. The steady state analysis in [12] assumed that the timescale of chemical reactions was fast compared to other timescales in the problem. Following in this spirit, we chose , , , and to be large and to produce dissociation constants which are roughly in line with literature values ( or ). We note that the phrase “calcium preferred” only describes the relative size of the two dissociation constants, and does not necessarily imply that there is actually more calcium bound to the network. The actual concentration of bound ions of each type depends on the chemistry, as well as the total amounts of each ion available for binding. To construct a network which is initially in a highly crosslinked state, we chose inner sodium and calcium concentrations that, for each case, yield a binding/unbinding equilibrium where most of the network is crosslinked (). Chloride concentrations were then specified to maintain electroneutrality. In each of these four regimes, we carried out swelling experiments for a variety of bath concentrations of sodium. Again, chloride concentrations were adjusted to maintain electroneutrality. The inner and bath concentrations, as well as the binding/unbinding rates and binding site density used in each of these four cases are listed in Table 2.

Table 2

Parameters that vary between experiments.

Parameter	K_Ca>K_Na,	K_Ca>K_Na,	K_Ca<K_Na,	K_Ca<K_Na,
	Dense Binding	Sparse Binding	Dense Binding	Sparse Binding
Total inner Ca	250 mM	25.3 mM	250 mM	25.3 mM
Total inner Na	5 mM	1 mM	5 mM	1 mM
Total inner Cl	5 mM	1.6 mM	5 mM	1.6 mM
kCaon	1×10^6 M^−1s^−1	1×10^6 M^−1s^−1	1×10^7 M^−1s^−1	1×10^7 M^−1s^−1
kCaoff	1×10^3 s^−1	1×10^3 s^−1	1×10^3 s^−1	1×10^3 s^−1
kNaon	5×10^6 M^−1s^−1	5×10^6 M^−1s^−1	5×10^5 M^−1s^−1	5×10^5 M^−1s^−1
kNaoff	5×10^2 s^−1	5×10^2 s^−1	5×10^2 s^−1	5×10^2 s^−1
z˜	1 M	0.1 M	1 M	0.1 M

3.2. Gel Swelling Experiments

Figure 2 shows an example time evolution of the network and solvent volume fraction, as well as dissolved and bound ion concentrations. The shown data were generated using a bath sodium concentration of M with and sparse binding sites. The initial condition was the same as that shown in Figure 1. The images show that, as time progresses, sodium from the bath diffuses into the inner region, reducing the bath concentration of sodium. As this occurs, the concentration of bound sodium drastically increases and propagates left, towards the interior region of dense network. This necessarily drives a decrease in the concentration of bound and doubly-bound calcium, and this release of bound calcium leads to an increase in the concentration of calcium dissolved in the solvent. As time progresses, the ionic concentrations (both bound and dissolved) approach uniform profiles. As the ions diffuse into and out of the region of network volume fraction, there is a rearrangement of the network phase. In particular, volume fraction of network flows from left (high density) to right. This decrease in maximum network density combined with an increase in the spatial region where network volume fraction exists is interpreted as swelling of the network phase. Eventually, over the course of several seconds, the network phase spreads across the domain and approaches a uniform steady state.

Figure 2

Example time evolution of: volume fraction (first column); bound ionic concentrations (second column); and dissolved ionic concentrations (third column).

To quantitatively measure the rate at which the network phase swells, we track the maximum network volume fraction () as a function of time. Due to conservation of mass, the volume fraction of network at steady state (when the gel is completely uniform) may be determined directly from initial conditions. We denote this quantity as and note that, in these experiments, . We may then define the quantity as a measure of the network swelling behavior. As the network swells towards its final configuration, approaches zero. Figure 3 shows the time evolution of for several bath concentrations of sodium. All of these data were produced in the case when binding is sparse and sodium binding is preferred (). It is immediately apparent that the bath concentration of sodium influences the rate of swelling, although the nature of this influence is not necessarily simple. The rate at which decreases does not appear to be monotonic with respect to bath sodium concentrations (we return to this below). Furthermore, the effect of bath sodium on swelling rate varies depending on the time one considers. For example, at times prior to roughly 1 s, the network swells more quickly when the bath sodium concentration is M (green solid curve in Figure 3) compared to when it is M (red dashed curve). However, at longer time scales, this behavior reverses, and the network immersed in a M sodium bath swells more rapidly. Finally, we note that, at longer times (greater than approximately 1 s), the behavior of appears exponential (linear on the semilog plot in Figure 3). Exponential swelling behavior has been predicted by numerous other studies. However, on time scales less than a second, the behavior of is distinct. This “short time” behavior is a novel prediction of the model, and something that we investigate further below. The transition between this “short time” behavior and the more obvious exponential decay is difficult to precisely quantify, but generically occurs sometime prior to 1 s. To focus on the short time swelling behavior, we turn to a separate quantification of swelling rate.

Figure 3

Decay of maximum network volume fraction to steady state for several concentrations of external sodium. Vertical dash-dot line indicates (roughly) the transition from early-time expansion behavior to long-time exponential behavior.

3.3. Front Propagation

For expositional purposes, suppose that a region of network of size and constant volume fraction swells to a new, steady state size L and volume fraction . Conservation of mass stipulates that, at all points in time, . If we then assume that volume fraction approaches its steady state value in an exponential manner then we may deduce that From Figure 2, it is clear that the network in our experiments does not have a spatially uniform profile, and thus does not meet the assumptions of the above calculation. However, this argument provides a scaling law for the size of the network globule which we investigate. Now, to define the width of the network globule in our experiments, we identify the location where network volume fraction is equal to We refer to the location as the “front” of the network, as there is 99% solvent (by volume) in the space to the right of this point. We know that the steady state width of the network is µm, as it fills the entire domain. Figure 4 shows an example time course of . The data shown were generated from the same simulation illustrated in Figure 2, where the bath concentration of sodium was M, , and binding sites were sparse.

Figure 4

Time evolution of network front towards right wall (at ). Dashed line indicates exponential fit from to s. The slope of this line ( s) represents the early-time expansion rate of the network.

Figure 4 shows several behaviors. At extremely short times, prior to approximately 0.2 s, decreases rapidly, indicating a rapid propagation of front location. At longer times, after approximately 1.4 s, we again see a rapid decrease in , corresponding to the front reaching the right end of the domain (this occurs in finite time, causing to reach zero, which cannot be illustrated on the semilog plot). It is worth noting that the front approaching the right hand end of the domain approximately corresponds to the transition between the “short time” and “long time” behaviors seen in Figure 3 (see dashed blue line). However, at intermediate times, we see an approximately exponential decrease in (again, indicated by an approximately linear behavior on the semilog plot). To identify the rate of this exponential decrease, we perform a linear fit to the quantity The slope of this linear fit can be interpreted as the exponential decay rate (see Equation (23)). To eliminate the contribution of the extremely early time start-up behavior, as well as the effects of the front reaching the right end of the domain, we only use data points between 0.2 s and 0.8 s. For reference, Figure 4 shows a dashed line indicating a purely exponential behavior with the decay rate s, which is the value calculated from the numerical simulation.

Figure 5 shows the calculated decay rate across all experiments. Red curves indicate cases where calcium binding is preferred (), while blue curves indicate networks that chemically prefer sodium binding (). Dashed lines indicate networks with sparse binding sites () and solid lines indicate networks with dense binding sites (). For extreme concentrations of sodium in the bath, the behavior of gel swelling is relatively easy to characterize. At very high bath concentrations of sodium, all networks appear to swell at the same rate. In this limit, bath sodium is so prevalent that all calcium crosslinks on the network are rapidly replaced with sodium binding (regardless of the relative size of dissociation constants, or number of binding sites), and the networks all swell in an identical manner. In the limit of low sodium concentration, we observe two clear patterns: a network with dense binding sites swells more slowly than one with sparse binding, and a network that chemically prefers calcium bonds swells more slowly than one that prefers sodium bonds. To understand this behavior, note that, if the bath concentration of sodium is very low, then the total number of sodium ions in the domain is severely limited. This means that a network with a high density of binding sites cannot break a large proportion of calcium crosslinks in favor of sodium bonds, and thus will not swell very rapidly (i.e., solid curves indicate slower swelling than their dashed counterparts in Figure 5). By a similar logic, if the number of sodium ions is small, then a network that prefers calcium bonds will not break crosslinks in favor of sodium binding, and will not swell rapidly (i.e., red curves indicate slower swelling than their blue counterparts in Figure 5).

Figure 5

Computed swelling rate as the bath concentration of sodium is varied. Data are shown for sparse (dashed) and dense (solid) network binding in the cases where sodium (blue circles) and calcium (red crosses) binding is preferred.

However, at moderate concentrations of sodium, we see regimes where these patterns are reversed. For networks that prefer sodium bonds (blue curves), for sodium concentrations between approximately and M, the network with dense binding sites swells more rapidly than the network with sparse binding sites. A similar behavior is observed for networks that prefer calcium bonds (red curves), though at higher bath concentrations of sodium (approximately to 2 M), as more sodium ions are required to displace calcium bonds. Similarly, for certain regimes of moderate bath sodium concentrations, we see that networks which prefer calcium bonds (red curves) actually swell faster than networks which prefer sodium bonds (blue curves). This behavior occurs between roughly and M for sparse binding sites and and 2 M for dense binding sites. This result is counter-intuitive, as one would suppose that chemically preferring calcium bonds would always retard the breaking of calcium crosslinks, and thus result in slower swelling behavior.

4. Discussion

In this paper, we have presented the first investigation of the dynamic swelling behavior of the model of polyelectrolyte gels first derived in [12]. We characterized the early-time swelling rate of the gel network as a function of the chemical parameters that govern the network’s ability and propensity to form crosslinks via calcium bonds, as well as the bath concentration of sodium. Our primary focus was on how the loss of calcium crosslinks alters the energetic landscape of the gel network, thereby driving a swelling event. The fact that swelling rate depends in a non-monotonic way on the concentration of sodium in the bath is a major result of this work. We do not currently have a simple explanation for this counter-intuitive behavior, and future investigations will address this issue. Experimental validation of this prediction would be exciting, but we know of no directly comparable studies. Experimental data would also provide an opportunity to estimate several parameters in the model (specifically, the interaction energies ) for which we have no reliable estimates.

Finally, we note that, for the sake of simplicity, we chose to treat the rheology of the mucus network as constant and viscous. However, it is known that mucus networks may exhibit a wide range of viscoelastic behaviors depending on the ionic milieu and other factors [21,22,23], including the presence of permanent covalent crosslinking within the gel network [24]. It may be that the formation and breakage of calcium crosslinks alters the rheology of mucus gel in a manner which we do not account for here. If so, the dynamics of swelling are likely governed by a complex interplay between the changing the internal energy and the evolving rheology of the gel mixture as its crosslinking structure varies in space and time. Such considerations are beyond the scope of this work, but may be included in future modeling attempts and potentially fully characterize the ionic regulation of mucus swelling in the human stomach.