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.


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

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.

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: ∂ ∂t θ n + ∇ · (θ n u n ) = 0, ∂ ∂t θ s + ∇ · (θ s u s ) = 0.
Here, θ n and θ s are the volume fractions of network and solvent, respectively, and u n and u s are their respective velocities. By definition, the volume fractions satisfy the constraint θ s ( x, t) + θ n ( x, t) = 1. As a consequence of this constraint, Equations (1) and (2) imply that the volume-averaged incompressibility condition ∇ · (θ n u n + θ s u s ) = 0, 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, ∇ · (θ n σ n ) − ξθ n θ s ( u n − u s ) − θ n ∇µ n − θ n ∇p = 0, Here, σ n and σ s 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; µ s , µ n , and µ j for j ≥ 3 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 ∇µ j 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), φ j 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 ν i is the viscosity of the respective phase.

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).

∂ ∂t
c Na + ∇ · (c Na u s ) = 1 θ s ∇ · D Na θ s (∇c Na + z Na c Na ∇Ψ) − k on Na θ n Mc Na + k off Na θ 2 s b Na , Here, c j is the concentration of species j measured in units of moles per liter solvent volume, D j is that ion's diffusion coefficient, z j its valence, and Ψ is the electric potential measured in terms of thermal voltage RT/F (R: ideal gas constant; T: absolute temperature; and F: Faraday constant). The parameter k on j is the binding rate of ion j to the gel network, while k off j is the corresponding unbinding rate. The ratio K j = k off j /k on j is the dissociation constant of that cation.

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, b Na 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 (b Ca ) and doubly bound (b C2 ). 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 wherez 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.

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 I in more detail below. Finally, the term labeled (III) represents the potential due to the charged nature of the mucus network.
Here, k B is the Boltzmann constant and N is the number of monomers in a typical polymer chain of mucin. The term σ I is the total ionic molality, given by where, again, φ j = c j /c water is the concentration of ion j divided by the standard molarity of water (c water = 55.5 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 z j is the valence of ion j. Equation (17) accounts for both entropic (I) and electric (II) potentials. We now return to the term I, 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 I is not a constant; it is given by where 1 and 2 are constants described below. The quantity α( x, t) is the ratio of the concentration of current crosslinks within the network (b C2 ( x, t)) to the maximum concentration of crosslinks possible (zθ n ( x, t))/2): Because α may vary spatially and temporally as the concentration of doubly bound calcium and network volume fraction do, so may I. 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 µ 0 n and µ 0 s (in Equations (14) and (15)) are the so-called "standard free energy" of the mucus phase and solvent phase, respectively. The solvent term µ 0 s 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 1 , 2 , 3 , and 4 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].

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 L = 25 µ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 θ n ,z 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 (c j and b j ) 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 x = 5 µ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. In all simulations, the initial volume fraction of network, the mechanical and energetic properties of the network (the parameters µ s , µ n , ξ, 1 , 2 , 3 , 4 and N), and the diffusion coefficients (D j ) of the ions were identical. The values of these parameters are listed in Table 1. In particular, the i 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 (α = 1), the network exhibits an interaction parameter I = 0, which does not drive swelling at all (this is analogous to a Flory interaction parameter of zero). However, when there is no crosslinking (α = 0), the network exhibits an interaction parameter I = −45, which is analogous to a large, negative Flory interaction parameter, and leads to rapid swelling. The parameter 3 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 4 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.
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 (z = 0.1) while a "dense" network has 1 mole of binding sites per liter of pure mucus (z = 1). We also investigated the case where a network-calcium bond is chemically preferred to a network-sodium bond (K Ca < K Na ), as well as the case where network-sodium bonds are preferred (K Ca > K Na ). 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 k on Ca , k off Ca , k on Na , and k off Na to be large and to produce dissociation constants which are roughly in line with literature values (10 −4 or 10 −3 ). 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 (α > 0.8). 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.   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 0.02 M with K Ca > K Na 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. To quantitatively measure the rate at which the network phase swells, we track the maximum network volume fraction (θ max n (t)) 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 θ ss n and note that, in these experiments, θ ss n ≈ 0.1. We may then define the quantity

Gel Swelling Experiments
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 (K Ca > K Na ). 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 0.002 M (green solid curve in Figure 3) compared to when it is 0.02 M (red dashed curve). However, at longer time scales, this behavior reverses, and the network immersed in a 0.02 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.

Front Propagation
For expositional purposes, suppose that a region of network of size L * (t) and constant volume fraction θ(t) swells to a new, steady state size L and volume fraction θ ss n . Conservation of mass stipulates that, at all points in time, L * (t)θ(t) = Lθ ss n . 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 0.01 θ n L * (t) = 0.01.
We refer to the location L * 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 L = 25 µm, as it fills the entire domain. Figure 4 shows an example time course of ∆ L −1 . The data shown were generated from the same simulation illustrated in Figure 2, where the bath concentration of sodium was 0.02 M, K Ca > K Na , and binding sites were sparse. Figure 4 shows several behaviors. At extremely short times, prior to approximately 0.2 s, ∆ L −1 decreases rapidly, indicating a rapid propagation of front location. At longer times, after approximately 1.4 s, we again see a rapid decrease in ∆ L −1 , corresponding to the front reaching the right end of the domain (this occurs in finite time, causing ∆ L −1 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 ∆ L −1 (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 quantitỹ L = ln(∆ L −1 ).
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 γ = 2.39 s −1 , 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 (K Ca < K Na ), while blue curves indicate networks that chemically prefer sodium binding (K Na < K Ca ). Dashed lines indicate networks with sparse binding sites (z = 0.1) and solid lines indicate networks with dense binding sites (z = 1). 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).
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 4 × 10 −3 and 2 × 10 −1 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 4 × 10 −2 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 4 × 10 −3 and 4 × 10 −1 M for sparse binding sites and 2 × 10 −1 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.

Swelling rates measured by front propogation
K Ca > K Na , Dense binding K Ca > K Na , Sparse binding K Ca < K Na , Dense binding K Ca < K Na , Sparse binding 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.

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 k ) 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. Here, we outline the numerical scheme used to simulate our system of of equations. To solve the momentum equations while respecting the volume-averaged incompressibility constraint (Equations (3) to (5)), we utilize a relatively standard second order finite difference method, in conjunction with the "ε-regularization" outlined in [25]. To evolve the volume fractions of network and solvent, we utilize a variation of second-order finite volume method (with van Leer limiter) [26]. To evolve the equations which govern ion transport and chemistry, we use a semi-implicit variation of the electrodiffusive solution scheme described in [27].
We begin by discretizing space using a so-called "staggered grid". Two collections of spatial points are defined, and various quantities more naturally "live" at each. The first collection of points we refer to as "cell centers," and they are defined by The second collection of points are referred to as "cell edges" and are defined by where ∆x = L/J is the spatial resolution of the grid. There are in total J − 1 cell edges and J cell centers. A schematic of the spatial discretization is shown in Figure A1. Finally, we have a single temporal discretization t k = k ∆t, k = 0, 1, . . . (A3) x 3/2 = ∆x We approximate ionic concentrations (bound and dissolved), as well as volume fractions and pressure at cell centers. Where necessary, we use a second subscript j to denote the spatial location where an approximation takes place, while a superscript k denotes the temporal location.
Because concentrations are defined at cell centers, particle fractions are as well.
Conversely, we approximate velocities, potential gradients, and the electric potential gradient at cell edges.
For brevity, we introduce the quantity V i to represent the potential gradient associated with species i, and Φ to represent the electric potential gradient.
The overall time integration scheme can broadly be described in the following steps carried out at each timestep: 1. Given ionic concentrations (C k i, j and B k i, j ), electric potential gradient (Φ k j ), and both volume fractions (Θ k i, j ), evaluate the potential gradients which act on the solvent and network phase (V k i, j ). 2. Given the potential gradients on each phase (V k i, j ), solve the momentum equations to determine the network and solvent velocities (U k i, j ). 3. Given their respective velocity fields, update the solvent and network volume fraction (Θ k+1 i, j ).
4. Given the two transport velocities and volume fractions (U k i, j and Θ k+1 i, j ), evolve the ionic concentrations while simultaneously solving for a new electric potential gradient (C k+1 i, j , B k+1 i, j and Φ k+1 j ).

Momentum Solve
To solve the momentum equations for the two velocity fields, we use the ε-regularization described in [25]. This method is designed to avoid complications where the linear operators in Equations (4) and (5) become ill defined in regions where θ n ≈ 0. Simply put, we define "regularized" volume fractions which avoids such an issuê We then use these volume fractions in discrete versions of the momentum and volume-averaged incompressibility constraint to avoid any poorly scaled linear operators. These "regularized" volume fractions are not used later when we update the actual volume fractions, but are simply used to create a regularization of the elliptic operators that define the velocity fields. Using standard second order finite difference schemes and moving all potentials (aside from pressure) to the right hand side, we can write Equations (4) and (5)  (A14) At cell edges which are adjacent to the domain boundary (j = 3/2 and j = J − 1/2), these equations contain terms which are not represented on our numerical grid. However, the closed domain and corresponding no-flux conditions on solvent and network imply that U k n, 1/2 = U k n, J+1/2 = 0. Similarly, Again, we stipulate that U k s, 1/2 = U k s, J+1/2 = 0. The discrete form of the incompressibility constraint is enforced "at" the cell centers. Therefore, for j = 1 . . . J, we have (for j = 3/2 . . . J − 1/2) and P k j (for j = 1 . . . J). Solving these equations allows us to determine the velocity fields and pressure at the k-th time step.

Volume Fraction Advection
Once the velocities and pressure field have been determined, we use the advection equations (Equations (1) and (2)) to update each volume fraction. To do this, we use a modification of a standard high-order advection scheme. Given solvent volume fractions at the cell centers (Θ k s, j , j = 1 . . . J) and velocities at cell edges (U k s, j , j = 3/2 . . . J − 1/2), we define fluxes at the cell edges using a van Leer limited high order method [26] (F k s, j , j = 3/2 . . . J − 1/2). We then update the volume fraction via the following relation: To enforce the no-flux boundary conditions consistent with a closed domain, we assume that F k s, 1/2 = F k s, J+1/2 = 0. Finally, we update the volume fraction of network via At this point, the concentration of dissolved ions must be numerically corrected. To illustrate, consider the concentration of chloride in the domain. The evolution equation (Equation (9)) (and the derivation of our model in [12]) conserve the total amount of chloride in the domain This is particularly important in problems of electrodiffusion, as a non-conservative numerical scheme may lead to a local imbalance of charge within the domain, which will result in extremely large electric potential gradients that may prove numerically untenable. However, even though Equation (A17) is a conservative scheme for the quantity Θ s , the update may have altered the discrete integral in Equation (A20). To correct this, we first calculate the flux quantitiesF k i, j (at cell edges) for each ion i by using the same van Leer limited high order method on the quantity Θ k s, j C k i, j (again, using solvent velocities). We then define the corrected concentration This corrected concentration ensures the discrete conservation of each ionic species At this point, having corrected any ionic imbalances that the volume fraction advection step may have created, we are able to take a time step of the electrodiffusive equations.

Ion Update
To time-step the evolution equations for ionic species, we use a methodology similar to that employed in [27]. The goal is to use an implicit-explicit (IMEX) time integration technique on all ionic species simultaneously, to avoid any further time-splitting errors. Unfortunately, the electric flux term and several reaction terms involve non-linear products of ionic concentrations and/or the electric potential gradient. Therefore, we treat these terms semi-implicitly in time. To do so requires extrapolating "forward" in time. We use first order extrapolation whenever possible, and zeroth order whenever we must (i.e., at the very first time step). Doing so allows us to define the quantities as well asB We then define F k+1 i, j as the fluxes of C * i, j due to solvent velocity (using the same van Leer limited method) and G k+1 i, j as the fluxes of B k i, j due to network velocity. We are now able to discretize Equations (7) to (9) for the evolution of each of the three dissolved ionic species. At cell centers which are not adjacent to the boundary (j = 2, 3, . . . J − 1), we have Na, j−1 At the first cell center (j = 1), the no flux boundary condition yields Na, 1 At the last cell center (j = J), the no flux boundary condition yields Na, J−1 Similarly, the equations for the dissolved calcium are (for j = 2 . . . J − 1) Ca, j−1 At the first cell center (j = 1), the no flux boundary condition yields  Because chloride is not allowed to bind to the network, the equations are relatively simple, as there are no reaction terms.

C k+1
Cl, j −C k Because the equations which govern bound ionic species are not diffusive, they are quite a bit simpler. For sodium, we have (j = 1 . . . J)
Equations (A25) to (A37) constitute 7J − 1 linear equations for the 7J − 1 unknown electrodiffusive quantities (J each for the C k+1 i, j s and B k+1 i, j s, and J − 1 for the Φ k+1 j s). Therefore, solving these equations simultaneously evolves the electrodiffusive model one time step.
Notice that we do not enforce electroneutrality at the right most cell center (j = J). This is because the quantities Φ k+1 j may be regarded as Lagrange multipliers on Equations (A25) to (A36) which serve to enforce the constraint in Equation (A37). As we only have J − 1 Lagrange multipliers, we are only able to enforce electroneutrality at J − 1 cells. However, as our numerical scheme is conservative, the total amount of each ion within the domain is conserved. Therefore, enforcing electroneutrality at all but one cell is guaranteed to enforce electroneutrality at the remaining cell.