A quantitative analysis of electrolyte exchange in the salivary duct

Kate Patterson, Marcelo A. Catalán, James E. Melvin, David I. Yule, Edmund J. Crampin, James Sneyd


A healthy salivary gland secretes saliva in two stages. First, acinar cells generate primary saliva, a plasma-like, isotonic fluid high in Na+ and Cl. In the second stage, the ducts exchange Na+ and Cl for K+ and HCO3, producing a hypotonic final saliva with no apparent loss in volume. We have developed a tool that aims to understand how the ducts achieve this electrolyte exchange while maintaining the same volume. This tool is part of a larger multiscale model of the salivary gland and can be used at the duct or gland level to investigate the effects of genetic and chemical alterations. In this study, we construct a radially symmetric mathematical model of the mouse salivary gland duct, representing the lumen, the cell, and the interstitium. For a given flow and primary saliva composition, we predict the potential differences and the luminal and cytosolic concentrations along a duct. Our model accounts well for experimental data obtained in wild-type animals as well as knockouts and chemical inhibitors. Additionally, the luminal membrane potential of the duct cells is predicted to be very depolarized compared with acinar cells. We investigate the effects of an electrogenic vs. electroneutral anion exchanger in the luminal membrane on concentration and the potential difference across the luminal membrane as well as how impairing the cystic fibrosis transmembrane conductance regulator channel affects other ion transporting mechanisms. Our model suggests the electrogenicity of the anion exchanger has little effect in the submandibular duct.

  • salivary glands
  • tubular transport
  • membrane transport

saliva is composed of 99% water, and yet that other 1% helps to maintain oral health by providing an appropriate ecological balance in the mouth (45). The basic secretory unit of a salivary gland consists of a cluster of acinar cells and a segment of duct (Fig. 1). The acinar cells secrete an isotonic primary saliva high in Na+ and Cl but low in K+. These cells produce all of the water found in saliva as the duct is considered to be highly impermeable to water (21). As the fluid travels along the duct, the duct cells exchange Na+ and Cl for K+ and HCO3, producing a hypotonic final saliva. This electrolyte exchange creates a final saliva that acts as a buffer in the mouth, neutralizing acids and inhibiting caries progression. This hypotonicity also enhances the ability to taste salty foods and nutrient-rich foods (18). Primary saliva production has been modeled by Gin et al. (16) and Palk et al. (37); however, these models do not include any representation of the salivary duct. One aim in studying this system is to provide a quantitative understanding of how hypotonicity is achieved in the mouse submandibular gland duct with no apparent change in saliva volume. It is, to the best of our knowledge, the first quantitative description of a salivary duct and is validated against a variety of experimental data. This model is also a tool that can be used alone or incorporated into a multi-scale model of the salivary gland, spanning from molecular to the tissue level (31, 37, 47). Depending on the application, either tool can be used to study genetic and chemical alterations, including, but not limited to, salivary gland disease, knockout studies, and gene therapy through viral vectors.

Fig. 1.

Basic secretory unit. A basic secretory unit consists of a cluster of acinar cells and a section of duct. Each acinar cell could be modeled as in Palk et al. (37) or Gin et al. (16). In this study, we describe how to model the duct. Above and below are enlarged basolateral and luminal membranes, respectively. The ion pathways are described in further detail in the Introduction as well as methods.

The pancreas is known to have a similar structure to the salivary gland, and both glands have been used as biological models to investigate secretory epithelium. One important difference between pancreatic ducts and salivary ducts is that while the pancreatic ducts secrete near isotonic fluid, the salivary ducts are impermeable to water and the final saliva is hypotonic. We have combined ideas from both the pancreatic duct cell models (48, 49, 52, 53) and the tubule transport models (11, 22, 51) to constructed a modular, multicompartment mathematical model of the salivary duct. The ion transporters are displayed in Fig. 1. Cl uptake occurs via an anion exchanger (AE; Ref. 58) located in the basolateral membrane (BLM) and a luminal cystic fibrosis transmembrane conductance regulator (CFTR) channel (6, 19, 55, 57) while efflux occurs through a basolateral ion channel (43) and a luminal AE (46, 58). Na+ enters the cell via a sodium-bicarbonate cotransporter (pNBC1) in the BLM (1, 4, 17, 30, 39) and through a luminal epithelial sodium ion channel (ENaC; Refs. 6, 12, 24). Na+ leaves the cell through a luminal sodium-bicarbonate cotransporter (17, 28) and via the Na+/K+-ATPase exchanger. HCO3 uptake occurs through the luminal AE and the basolateral pNBC1. The HCO3 efflux mechanisms are a luminal NBC, a basolateral AE, and a basolateral HCO3 channel. HCO3 is also allowed to move through the CFTR channel, albeit at a lower conductance than Cl. Finally, K+ efflux occurs through ion channels (KCa1.1 in the LM) and uptake occurs via the Na+/K+-ATPase exchanger (refer to the methods for more details).

We validate our model against experimental data. Because of the modular nature of our model, we are able to investigate the effects of an electroneutral vs. an electrogenic AE in the luminal membrane (LM) and the tight junctional conductance for cations. Our model predictions are consistent with the CFTR deletion study by Catalán et al. (6) in which they find that the CFTR deletion also results in a reduction in functional ENaC sodium channels. Additionally, we predict 1) the LM potential of the duct must be very depolarized; 2) the CFTR channel is absorbing Cl from the lumen while the AE is secreting Cl (contrary to the currently proposed model; Ref. 5); and 3) the tight junction must have a low permeability for ions.


We employ a radially symmetric advection-diffusion model to a fixed radius (RA) lumen with length L, closed at the acinus end. Surrounding the lumen is a ring of epithelia with radial height, h. We use the following notation where x is the position along the duct and t is time; c is one of the following solutes: Na+, Cl, K+, or HCO3; k is a placeholder for i, l, or e representing the cellular, luminal, or interstitial compartment; † is A or B representing the LM or BLM, respectively; [c]k(x, t) is the concentration of the solute c in compartment k; v(t) is the fluid velocity, positive in the direction of increasing x; N†,c(x) is the solute flux per area, across membrane †; and q†,c(x) is the osmotic flux per area, across membrane †. We follow the convention that a current or flux is positive when positive ions leave the cell.

Using standard conservation arguments, we derive differential equations for concentration, cell volume and luminal fluid velocity. [c]lt=2NA,cRA+D2[c]lx2v[c]lx,1 d([c]iw)dt=2π[(RA+h)NB,c+RANA,c],2 dwdt=2π[(RA+h)qBRAqA],and3 vx=0.4 For clarity, the variable dependence, (x, t), has been suppressed. Assuming a cylindrical lumen, [c]l depends on the solute transport across the LM, and the diffusion and fluid flow along the duct. We assume a continuum epithelium model for the cells. RA + h is the distance from the center of the lumen to the BLM. We solve for h using the formula h = w/π+RA2RA, allowing the cell height to change as water flows across the BLM. The cellular concentration, [c]i, and the cross-sectional area, w = π(2hRA + h2), are then calculated from the radial ionic (N†,c) and osmotic (q) fluxes, respectively. This model does not allow diffusion or transport between cells.

We assume that the LM is rigid, is essentially impermeable to water (qA = 0), and has a fixed radius and that saliva is an incompressible fluid. As a result, the fluid velocity in the lumen will be equal to the primary saliva flow rate, v0, at every point along the duct (uniform in x).

Boundary conditions.

We set the concentrations at x = 0 to the primary saliva composition determined experimentally by Mangos et al. (34). We do not have a primary saliva concentration for HCO3 but will assume that saliva is electroneutral and that the charge is carried by Na+, Cl, K+, and HCO3. We then use electroneutrality to determine the primary saliva bicarbonate concentration: [HCO3]l = [Na+]l + [K+]l − [Cl]l. Based on the observations of Martin and Young (35) and Mangos et al. (34), this is a reasonable assumption.

For our second boundary condition, we assume that for a sufficiently long duct, x = L, the concentration in the lumen, at the end of the duct, will attain some steady state. This steady assumption is supported by experimental measurements (see Fig. 5 in Ref. 34). The saliva composition is then calculated on a more realistic length, x = Lf. We found L = 2.5 mm to satisfy the steady-state assumption and estimated Lf = 1 mm from Flint (15). By solving the system in this manner, the concentration profile on the region of interest is not artificially affected by the steady assumption at the end of the duct.

Numerical simulations.

To solve this system of differential equations, we employ a second order, backward differentiation formula with a first order, upwind scheme. For the very first two time steps in a simulation, we calculate one based on the initial condition and the second using a single step of a first order backward Euler, central difference method. We use the steady-state values obtained under a 1 μm/min flow and unstimulated primary saliva concentrations as the initial condition along the duct. In Table 1, we list the primary saliva concentration.

View this table:
Table 1.

Initial conditions and their sources

For these numerical methods, the duct is partitioned spatially into segments of length Δx so that [c]l,0 is the concentration at x0 (this is primary saliva) and [c]l,n is the nth partition along the duct. The cytosolic concentrations are computed using finite differences for each cell of width Δx at each time step. These computations have been programmed and computed in MATLAB.

Luminal, cellular, and interstitial transport.

The flux of ions traveling from the cell via ion channels, exchangers and cotransporters is represented by N†,c. We represent the current through an ion channel by I†,c and the flux from an exchanger or a cotransporter is represented by J†,*(σ) we distinguish an electroneutral from an electrogenic flux by σ = 1 or σ = 2, respectively, and * is a place holder for the ions involved in the flux. The various ion transport mechanisms are shown in Fig. 1. The individual currents and fluxes are described in more detail below.

Thus N†,c represents the number of moles of ions traveling from the cell per area and has units of mol·s−1·m−2. NA,[Na]=IA,NazNaF+JA,NaHCO3(1),NA,[Cl]=IA,ClzClF+JA,ClHCO3(σ),NA,[K]=IA,KzKF,NA,[HCO3]=IA,HCO3zHCO3FσJA,ClHCO3(σ)+JA,NaHCO3(1),NB,[Na]=IB,NazNaF+3JB,NaK+JB,NaHCO3(2),NB,[Cl]=IB,ClzClF+JB,ClHCO3(1),NB,[K]=IB,KzKF2JB,NaK,NB,[HCO3]=IB,HCO3zHCO3FJB,ClHCO3(1)+2JB,NaHCO3(2).

The BLM has permeability LB and we represent the osmotic flux across the BLM as: qB=RTLB{[Na+]e+[Cl]e+[K+]e+[HCO3]e([Na+]i+[Cl]i+[K+]i+[HCO3]i+χw)}, where χ/w is the concentration of impermeable ions and RT is the gas constant and temperature. Because the duct is nearly impermeable to water (45), we set the LM osmotic flux to zero, qA = 0.

Membrane potential difference.

The potential difference as measured from the exterior of the cell to the cytoplasm is represented as follows: CmdVAdt=IA,NaIA,KIA,ClIA,HCO3F(σ1)JA,ClHCO3(σ),5 CmdVBdt=IB,NaIB,KIB,ClIB,HCO3F[JB,NaK+JB,NaHCO3(2)]6 for the LM and BLM, respectively. Cm is the cell capacitance per membrane area in F μm−2. As the cell capacitance for a duct cell is quite small [4–12 pF Zeng et al. (57) and personal communication with the authors of Catalán et al. (6)], we will solve Eqs. 5 and 6 using a quasi-steady-state approximation. This approximation assumes that the potential difference adjusts to changes in concentration immediately, while the ion fluxes occur more slowly. As such, we set dVA/dt = 0 and dVB/dt = 0 and solve for VA and VB algebraically.

Ion channels.

We model transmembrane ion currents, I†,c, as: I,c=g,c(VV,c)whereV,c=RTzcFlog([c]k[c]i).7 For the LM, † = A and k = l; for the BLM, † = B and k = e. zc is the valence of the ion c, g†,c is the conductance per membrane area, and V†,c is the Nernst potential for the channel. As is standard, R, T, and F are the universal gas constant, the temperature, and Faraday's constant, respectively.

We use the whole cell, ENaC- and CFTR-mediated current-voltage data from Catalán et al. (6) to estimate the conductances for the sodium and chloride channels (ENaC and CFTR, respectively) given in Table 2. The CFTR channel is known to be permeable to other anions such as HCO3 (29, 40, 41); we use the permeability ratio from Poulsen et al. (40) to determine HCO3 conductance through the CFTR channel.

View this table:
Table 2.

Conductance parameters

The only potassium channel reported in the LM is the outward rectifying, maxi-K channel (36). To our knowledge, whole cell, KCa1.1-mediated current-voltage data have not been measured for a duct cell. We approximated the potassium channel conductance using the whole cell cation I-V plot for a granular duct cell published by Nakamoto et al. (36) (Table 2).

The BLM contains an ion channel for each ion. These conductances have been chosen to approximate the potential difference and cytosolic concentrations as listed in Table 1.

Exchangers and cotransporters.

We derive a general three-state model similar to Palk et al. (see Ref. 37, Appendix E). We assume the exchanger/cotransporter has two distinct conformations, I and O, where I is the exchanger/cotransporter with ions bound and O is unbound. We assume that the ions bind and release simultaneously. With the use of C1 and C2 to represent the ions, the three-state model can be described as O+C1k1k1+Ik2k2+O+C2.8 From this model we derive a system of differential equations and determine the turnover rate at steady state to be: J,*(σ)=g,*(k1k2C1k2+k1+C2)k2C1+k2++k1+C2+k1.9 The constant σ represents the stoichiometry of the exchanger/cotransporter. If it is electroneutral, σ = 1, otherwise σ = 2. This general formulation of the turnover rate will be applied below to an electroneutral AE, an electrogenic AE, and the luminal NBC.

For JA,ClHCO3(1), the electroneutral AE in the LM we use Eq. 9 with C1 = [Cl]i[HCO3]l, C2 = [Cl]l[HCO3]i, k1+ = 0.02 s−1·mM−2, k1 = 0.1 s−1, k2+= 0.0067 s−1, and k2 = 0.018 s−1·mM−2.

For the basolateral AE JB,ClHCO3(1), we assume the same form and rate constants as in JA,ClHCO3(1), replacing luminal concentrations with interstitial concentrations, [Cl]e, [HCO3]e and the density with gB,ClHCO3.

The electrogenic AE (stoichiometry of 1Cl:2HCO3) is represented by Eq. 9 with C1 = [Cl]i2[HCO3]l and C2 = [Cl]l2[HCO3]i. To keep the focus on the electrogenic contribution, we held the rate constants and densities the same as the electroneutral exchanger, only introducing the electrogenic step: k1+ = 0.02 exp (FVART) s−1·mM−3, k1 = 0.1 exp (2FVART) s−1, k2+= 0.0067 s−1, and k2 = 0.018 s−1·mM−3.

For the NBC in the LM we use Eq. 9 with C1 = [Na+]l[HCO3]l, C2 = [Na+]i[HCO3]i, k1+ = 0.005 s−1·mM−2, k1= 0.03 s−1, k2+ = 0.002 s−1, and k2= 0.08 s−1·mM−2. The rate constants and densities were fit to attain mouse SMG final saliva.

While we could design a similar three state electrogenic model for the pNBC found in the BLM, we instead implement the carefully derived, previously tested, pNBC turnover rate from Sohma et al. (49) for the basolateral pNBC, JB,ClHCO3(1). We use their parameters, and assume the stoichiometry is 1Na+:2HCO3. The density, gB,NaHCO3, is fit to this model.

To represent Na+/K+-ATPase exchanger, we implement a simplified model [see Appendix F in Palk et al. (37)] for further details. The density, gB,NaK, is fit to this model.

Saliva composition.

Stimulation of acinar cells, by applying an agonist such as carbachol (CCh), increases fluid production and slightly changes the primary saliva composition (34). We have made the assumption that CCh concentration positively correlates to the fluid flow. To model stimulation, we increase v0 from 1 to 8 μl/min for 0.3 μM CCh + 5 μM isoproterenol (IPR) and to 20 μl/min for 0.5 μM CCh. In addition, we alter the concentration of the primary saliva, or equivalently, the boundary condition at the beginning of the duct, as shown in Table 1.

We convert the flow rate from μl/min to μm/s assuming the average SMG has a mass of 67 mg, that there are 2.26 × 107 acinar cells per 100 mg (37), the duct has a radius of 3 μm, and that an acinus has six cells (Table 3).

View this table:
Table 3.

Duct parameters

To compare our results to Nakamoto et al. (36) and Catalán et al. (6), we compute the luminal solute concentrations on a duct of length L = 2.5 mm for 10 min under unstimulated conditions to allow the system to reach a steady state. We then stimulate the gland for 10 min and begin collecting the final saliva. After 10 min, we stop stimulating the gland. For the Nakamoto et al. (36) experiments, the saliva collection period is complete; for the Catalán et al. (6) experiments, we continue collecting unstimulated saliva for a further 10 min.

Concentrations in the collected final saliva are estimated by computing the average at x = Lf over the collection period. [c]f=0T[c]l(Lf,t)vπRA2dt0TvπRA2dt,10 where [c]f represents the final saliva concentration for solute c; t is from 0 to T, and T = 10 or T = 20 min depending on the experimental collection time; [c]l(Lf, t) is the concentration of the solute at x = Lf; v is the fluid velocity; and RA is the radius of the duct. The choice of Lf (determined in Boundary conditions) will affect the final saliva composition; the smaller the value of Lf, the less duct available for electrolyte exchange, resulting in higher NaCl and lower K in the final saliva. Table 4 lists the WT final saliva composition at Lf = 1.0 mm and at Lf ± 10%.

View this table:
Table 4.

Wild-type final saliva composition when Lf is 0.9, 1.0, and 1.1 mm

View this table:
Table 5.

Steady-state solution at x =1.0 mm, v0 =1 μl/min


Wild-type model.

In an early model of the salivary duct, we modeled the AE as the primary mechanism for absorbing chloride ions while the CFTR channel recycled Cl back into the lumen. This is the traditional model of chloride transport in the duct [see for example Boron and Boulpaep (5)]. In the wild-type model, we could obtain the appropriate final saliva concentrations if we included an active potassium flux. However, we were unable to capture the observed increase in chloride concentration associated with a CFTR channel deletion or block. Once the CFTR channel was compromised, modeled by a severe decrease in conductance, the model predicted a decrease in final saliva chloride, contrary to the data (see Fig. 3B). Due to this result, we determined that the CFTR channel must be absorbing chloride. Therefore, the potential difference across the LM must be depolarized enough to allow for chloride absorption until final saliva chloride concentrations are attained. Further, for the CFTR channel and the AE to reach a steady state, the AE must be moving chloride from the cell to the lumen at low luminal chloride concentrations.

With a luminal potential difference of about −20 mV in the wild-type model, a potassium channel is sufficient to attain final saliva potassium composition. Further, any CFTR channel dysfunction results in increased final saliva chloride, as desired.

For our parameter values, we incorporate channel conductances from data whenever possible and fit the remaining parameters (Table 2). This model predicts final saliva composition under two different experimental conditions: 0.3 μM CCh + 5 μM IPR (6) (Figs. 2 and 3) and 0.5 μM CCh (36) (Fig. 4). Additionally, the model predicts that the cytosolic concentrations are essentially constant along the length of the duct. At steady state (v = 1 μl/min), the maximal changes along the duct are 0.1, 0.3, 0.3, and 0.2 mM, for [Na+]i, [Cl]i, [K+]i and [HCO3]i, respectively.

Fig. 2.

Final saliva comparisons between model simulations and data for wild-type and chemically blocked channels. A: application of amiloride blocks the ENaC sodium channel. We represent this as a reduction in the sodium channel conductance (gblock Na+ = 0.2 × gA,Na+). B: cystic fibrosis transmembrane conductance regulator (CFTR) channel chemically blocked (gblock Cl= 0.15 gA,Cl).

Fig. 3.

Model variations against data from the wild-type and CFTRΔF/ΔF mouse final saliva. CFTRΔF/ΔF mouse model has the CFTR channel deleted. A: CFTR channel conductance is reduced (gΔCl = 0.15 × gA,Cl). Notice how the model sodium concentration for the deletion is unable to capture behavior as measured in the CFTRΔF/ΔF mouse. B: in addition to the CFTR channel reduction performed in A, the ENaC channel conductance is also reduced (gΔNa+ = 0.4 × gA,Na+).

Fig. 4.

A: comparison between our model and the wild-type and KCa1.1 channel deletion final saliva (36). The KCa1.1 channel deletion is represented by a reduction in channel conductance, gΔK+ = 0.01 × gA,K+. We also reduce the CFTR conductance due to the absence of IPR (36), gNo IPR Cl = 0.5 × gA,Cl. A further reduction in the channel conductance will raise the final Cl concentration. B: AE activity is reduced and compared with the CFTRΔF/ΔF mouse. In addition to the parameter modifications in Fig. 3B, the AE activity is also reduced to 60%, gΔClHCO3 = gA,ClHCO3× 0.6 (in black); to 40%, gΔClHCO3 = gA,ClHCO3 × 0.4 (in gray); and to 20%, gΔClHCO3 = gA,ClHCO3× 0.2 (in white). Slight decreases in Cl and increases in HCO3 occur as the activity is reduced.

We have tested the sensitivity of our parameter choices by comparing the steady-state final saliva resulting from changing one parameter by ±10% to the wild-type steady-state saliva composition. The parameters tested were the densities and conductances in Table 2 as well as each of the rate constants in the three-state AE and NBC (see methods). The ENaC conductance proved to be the most sensitive resulting in a maximal change of ∼11%. All other parameters resulted in changes ≤7%.

Channel blocks and deletions.

To validate our model, we compute the final saliva composition and compare to the results of Catalán et al. (6) and Nakamoto et al. (36). In these studies, they investigate how various channel deletions and blockers affect the mouse SMG duct by measuring the final saliva composition of the gland. A channel deletion or block can be modeled by a reduction in the conductance. Since we do not know how strongly each deletion or chemical block affects the channel, we determine the reduction in conductance by fitting our model to the data. The conductance parameters we used to fit the model to the data are listed in Table 2.

In the bar charts depicting final saliva ion composition (Figs. 24), the solid bars represent the model while the striped bars are data from Catalán et al. (6) and Nakamoto et al. (36). We represent the following in order from left to right: sodium, potassium, chloride, bicarbonate, and osmolality. The wild-type or control experiment is represented by a white bar for the model, or a sparsely striped bar for the data, while the blocked or deleted channel experiment is represented by a shaded bar for the model and a densely striped bar for the data. When the bicarbonate or osmolarity data were not published, as in Catalán et al. (6), we have approximated the bicarbonate value under the assumption that [HCO3]l = [Na+]l + [K+]l − [Cl]l, i.e., that the luminal composition is electroneutral. However, our approximation may be slightly off as other ions are likely involved.

Catalán et al. (6) focus on the deletion and block of the CFTR channel and the block of the ENaC channel. These are the primary chloride and sodium channels found in the LM of the salivary duct cells. We show their ex vivo results in Figs. 2, 3, and 4B.

ENaC channel block.

Amiloride is known to block the ENaC sodium channel, resulting in an increased NaCl concentration in final saliva. We model this chemical block by reducing the value of gNa+, the ENaC channel conductance, by 80%. Our model accurately captures the experimental findings, Fig. 2A.

CFTR channel block.

We next consider the chemical block of the CFTR channel. In the experiment, Catalán et al. (6) use a highly selective CFTR channel inhibitor, CFTRinh-172. We find a reduction in the conductance to 15% of the wild-type CFTR channel conductance provides a good fit the CFTR blocked channel data (Fig. 2B).

Our model, consistent with the data, finds an increase in luminal chloride, and a slight increase in sodium. Our model also predicts an increase in potassium. This potassium increase is due to the LM becoming even more depolarized, an event which does not seem to happen in the mouse model (compare potassium concentrations from the model and data in Fig. 2B).

CFTR mutation.

For both blocked channels, our model is in good agreement with the data. The final dataset we consider from the study of Catalán et al. (6) is a CFTR mutation. Interestingly, the final sodium concentration in the data is much higher in the CFTR mutation (Fig. 3) when compared with the CFTR channel block (see Fig. 2B).

In Figs. 3 and 4B, we compare several variations of our model to the CFTR ΔF508 dataset from Catalán et al. (6). We first present the model with the CFTR channel conductance greatly reduced to represent the cystic fibrosis (CF) mouse model carrying the ΔF508 mutation in the CFTR channel. As Fig. 3A shows, we are unable to accurately capture the sodium concentration for the CF mouse model. Note that for ease of comparison, we reduce the conductance by the same value we used in the CFTR channel block in Fig. 2B, where we were able to represent the data accurately.

Catalán et al. (6) also note the difference in sodium concentration between the CF mouse model and the CFTR channel block. They find that the whole cell ENaC-mediated sodium inward current decreases and the amount of α-ENaC protein, a protein necessary for ENaC functionality (3), is reduced in the mouse model carrying the ΔF508 mutation. We reduce the ENaC channel density by 60%, in addition to the CFTR reduction, and find we are able to accurately capture the sodium concentrations reported in the CF mouse model (6) (Fig. 3B).

KCa1:1 channel deletion.

Nakamoto et al. (36) find the KCa1.1 channel is the primary pathway involved in K+ efflux. In Fig. 4A we compare our model to the data published in Nakamoto et al. (36). This data comes from a gland stimulated with 0.5 μM CCh (36), whereas the Catalán et al. (6) experiments use 0.3 μM CCh and IPR. The higher level of agonist suggests a faster flow rate, resulting in less NaCl absorption. However, as Romanenko et al. (43) showed, when IPR and CCh are applied, the final saliva chloride concentrations are considerably lower than when IPR is not included, presumably because CFTR is activated by IPR. To obtain the experimental results, the KCa1.1 deletion is represented here as a reduction to 1% of the wild-type KCa1.1 channel conductance. We have also reduced the CFTR conductance by 50% to account for the different agonist response.

Cl/HCO3 exchanger.

Zhao et al. (58) demonstrated the existence of a chloride-bicarbonate exchanger in the luminal and BLMs of duct cells. Shcheynikov et al. (46) found the electrogenic chloride-bicarbonate exchanger SLC26A6 to be responsible for the majority of bicarbonate transport in the LM of parotid duct cells. However, Shcheynikov et al. (46) also show immunolocalization of SLC26A4, an AE thought to be electroneutral (36a), to the LM of submandibular ducts.

We compare the effect of an electrogenic (1Cl:2HCO3) vs. an electroneutral AE. Then, we investigate the dependence of the AE activity on the CFTR channel. In Fig. 5A, we compare the luminal concentration profiles along the first 1.0 mm of the duct (at steady state) produced by the model with an electrogenic vs. an electroneutral AE in the LM. In Fig. 5, left, the electrogenic AE, JA,ClHCO3(2), is active while Fig. 5, right, has the same parameter values except the AE is electroneutral, JA,ClHCO3(1). We have not refit any parameters to obtain this data. Some differences exist along the duct, but nothing that clearly suggests which stoichiometry should be used in the duct. At 1.0 mm, the electrogenic exchanger results a decrease of 1.3 mM Na+, 6.5 mM K+, and 13.7 mM HCO3 in the lumen, while the Cl increases 5.9 mM. The potential difference across the LM is also decreased 3.5 mV. We have adopted the electroneutral AE for the figures in this study.

Fig. 5.

A: luminal concentration and potential difference along the duct for the electrogenic and the electroneutral AE. The luminal concentration profile from our model has the same behavior as in Fig. 5 from Mangos et al. (34). Here we show the saliva composition spatially along the first 1.0 mm of the wild-type salivary duct at steady state under low flow (1 μm/min), unstimulated conditions. Notice how the bicarbonate profile is nearly constant along the duct. B: concentration of ions vs. flow rate. Circles and squares represent ion concentration at 1.0 mm along the duct for the given flow rate. Units are scaled so that a flow rate of one (dashed line) is equal to the flow rate we used when modeling the Catalán et al. (6) data. As the flow rate increases (along the x-axis), it is clear that the Na+ and Cl absorption decreases, resulting more Na+ and Cl in final saliva. Similarly, K+ secretion decreases, resulting in less K+ in final saliva. This is comparable to Fig. 1 in Mangos et al. (34).

The CFTR channel activity has also been tied to the activity of the AE (46, 50). In Fig. 4B, we have reduced the electroneutral AE density, gClHCO3 to 60% (black bars) to 40% (gray bars) and to 20% (white filled bars) of the wild-type value; all other parameters are as in Fig. 3B. As the AE activity is reduced, the final saliva Cl and HCO3 are affected but not greatly.

Duct LM is depolarized.

When the CFTR channel is mutated or blocked, the chloride concentration found in final saliva is greatly increased (6, 19). To investigate the requirements this places on the LM potential, consider the luminal Cl current in Eq. 7. To obtain the desired final saliva concentrations, the CFTR channel current must be positive until [Cl]l ≈ 38 mM (34, 6). With the use of Eq. 7 and [Cl]l = 38, and [Cl]i = 25 mM (56) with the restriction that the current must be positive, we estimate VA ≥ −11 mV. This seems very depolarized, but by reducing [Cl]i to 18 mM, we can decrease this bound to VA > −20 mV. We are unaware of any experimental results reporting a salivary duct luminal potential difference.

Tight junction has low permeability to cations.

We explored the possibility of allowing cations to move across the tight junction, by implementing the following current: It,c=gt[(VAVB)RTzKFlog[c]l[c]e],11 where gt is the tight junctional conductance and c is either Na+ or K+. We find we can attain the final saliva measured by Catalán et al. (6) and Nakamoto et al. (36) when the tight junctional conductance is gt = 0.1 × gA,K+. This value of gt correlates to a resistance of 70,000–250,000 Ω/cm2, consistent with other tight epithelia (2). However, if the tight junctional conductance is increased to gt = gA,K+, the resulting final saliva is high in Na+ (99–115 mM) and low in K+ (24–25 mM), when simulating the Catalán et al. (6) and Nakamoto et al. (36) experiments, respectively.

For the simulations presented in this study, we have not included any tight junctional current. We have observed through simulations that we can include a tight junctional current, but it must be at least an order of magnitude smaller than the KCa1.1 cannel conductance. In addition, the duct is considered to be a very tight epithelium, also suggesting there is essentially no tight junctional flow.


There are three major salivary glands in most animals, the parotid, the sublingual, and the submandibular (or submaxillary) glands. The parotid gland produces a serous, watery secretion and is the largest of the glands. The other two glands produce a serous and mucous solution. Whenever possible, we used data from the SMG from the mouse. We choose the SMG gland because the majority of the experimental data on the salivary ducts seem to come from this gland. One possible reason for this preference may be because the NaCl absorption is more extreme compared with the parotid gland (34, 33, 32, 54).

We chose the mouse model for a similar reason. Due to the high density of the CFTR channel in the salivary ducts, researchers often use salivary ducts to study cystic fibrosis. Because of the mouse cystic fibrosis model system (55), many of the recent salivary duct experiments have been performed on the mouse model (6, 7, 20, 26, 36, 56). However, mouse SMG duct is not necessarily a good model for all salivary gland ducts in all animal models. Young and Schneyer (54) reviewed the various animal models and their differences. At maximal rates of pilocarpine stimulated secretion, rat SMG saliva is reported to have 5 mM Na+ and 35 mM K+, while the rat parotid saliva has 115 mM Na+ and 10 mM K+, and, finally, monkey SMG saliva has 90 mM Na+ and 40 mM K+. Just considering these three, the variation from animal to animal in the same gland or gland to gland in the same animal is so large, one model cannot hope to describe all possible scenarios. However, with the model that we have developed, it is straightforward to change conductances and densities or to introduce/remove an ion transporting mechanism. It is in this way that we can use this model as a tool and apply it to other projects. For example, one avenue of particular interest is studying the effects of genes introduced using viral vectors, such as ductal aquaporins, (10). In general, the structure remains the same, so the only changes necessary are the primary saliva compositions and the ion transporting mechanisms found along the duct.

Comparisons to other epithelial models.

To the best of our knowledge, there does not exist another model of a salivary duct or salivary duct cell. However, there are many other models of tubules with ion transport across the membrane; for some examples, see Weinstein (51), Dias et al. (11), and Keener and Sneyd (22). There are also other models of exocrine gland duct cells, some examples of which are Sohma et al. (48, 49), Whitcomb and Ermentrout (52), and Yamaguchi et al. (53).

Our model, in construction, is similar to other tubule transport models in that we have some flow along a tube and ion flux across the membrane. In particular, the ascending loop of Henle in the kidney is water impermeable and ion absorbing. However, the set of ion transporting mechanisms known to exist in duct cells differs from any of the tubule models we are aware of.

The exocrine duct cells models (48, 49, 53, 52) are all of pancreatic cells. They represent ion transport mechanisms in the LM and BLM and as a result have a system of differential equations for the cytosolic and luminal concentrations. The major differences between the pancreatic models and our model are 1) a water permeable vs. impermeable membrane, 2) that the final pancreatic juice is near isotonic and high in NaHCO3 whereas final saliva is hypotonic and low in NaCl, and 3) that they model a single cell while we model a section of duct.

CFTR channel and membrane potential.

Although the CFTR channel has been modeled as a Cl secreting mechanism in the pancreas, we find the CFTR channel in the duct must be Cl absorbing. In particular, if Cl ions are being secreted by the CFTR channel, when the channel is blocked or deleted, the concentration in the lumen will decrease. However, that is not what has been observed in the SMG duct (6, 19). In further support of this, Ishibashi et al. (19) injected either CFTR-siRNA or CFTRinh-172 into the duct lumen and reported an increase in final saliva Cl in both experiments.

Because it has been previously proposed that the CFTR channel and the AE are reversed, where Cl is secreted through the CFTR and absorbed through the AE (5), we investigated our model under these conditions. For the CFTR channel to secrete Cl, the membrane potential difference must be larger. With this change several problems arose. First, we were unable to find an ion with enough charge difference across the LM to set the potential difference. Furthermore, any impairment of the CFTR channel resulted in a decrease, rather than increase as seen in the block and deletion data (6), in the final saliva Cl concentration. Even if we set the potential difference artificially, we were unable to find a set of valid parameters in which the impaired CFTR channel resulted in an increase in final saliva chloride.

In Ko et al. (23), using a transepithelial voltage and the BLM potential measurements from the literature, they estimate the LM potential for a pancreatic cell to be −20 to −30 mV. However, we are not aware of similar measurements in the salivary duct. With the final saliva being low in sodium and chloride while high in potassium, the membrane potential could vary greatly from a pancreatic cell. In addition, the transepithelial voltage reported for the salivary duct ranges greatly, from −71 mV in rat to −13 mV in rabbit (14), and, to the best of our knowledge, no BLM potential has been reported. Given the available evidence, we predict that the LM of the duct will be depolarized.

Another variation on the behavior of the CFTR channel has been suggested by Park et al. (38). They have observed that the ratio of HCO3 permeability to Cl permeability, through the CFTR channel, can change. They propose that in the proximal duct, the CFTR channel is primarily a Cl channel but smoothly changes to primarily a HCO3 channel by the distal duct. In our model, we are able to attain reasonable values for bicarbonate concentration in the final saliva with constant CFTR conductance along the duct and hence do not find it necessary to consider variation of the current carried by CFTR in this manner.


We have investigated either an electroneutral or electrogenic (1Cl:2HCO3) AE. To the best of our knowledge, the electrogenic SLC26a6 has only been found in the parotid ducts, not the SMG ducts. Here, we aim to model the SMG duct, but we include the electrogenic exchanger to understand the affect, if any, the electrogenic exchanger may have. In particular, we wanted to see how the exchanger affected the membrane potential and if the direction of Cl through the CFTR channel could be reversed (Cl secreting) while still achieving an increased Cl concentration in the case of a blocked or deleted CFTR channel. Having explored all of the deletions and blocks, using the exact same parameters and conditions except for the AE representation, we found the differences in the final saliva concentration to be so near the experimental error that the two models are essentially the same. It is possible that in the parotid ducts the electrogenic AE has a more obvious contribution to the final saliva, but for the mouse SMG, the electrogenicity of the AE has little effect on the overall behavior of the duct.

We have fit the rate constants as well as the densities of the ion exchangers and cotransporters to the data. For each parameter (rate constants, conductances, and densities), we have tested the sensitivity by calculating the final saliva obtained when the parameter is altered by ±10%. The most sensitive parameter is the ENaC conductance, which results in a change as large as 11% in final saliva sodium. All other parameters result in ≤7% change.

If we set the membrane potential artificially so that the CFTR channel is secreting, we can find parameter values that generate realistic SMG saliva concentrations. However, we are unable to capture the experimental results in the case where the CFTR channel is blocked or deleted. Even when we reduce the AE activity by 40, 60, or 80%, if the CFTR channel is secreting, the CFTR channel deletion causes a decrease in final saliva Cl, opposite to the experimental observations. For the parameter values listed in Table 2, the CFTR channel absorbs Cl into the cell, and the AE secretes Cl into the lumen at steady state.

We investigated reducing the AE flux in the LM for the CFTR ΔF508 mutation but not for the CFTR channel block. Wang et al. (50) found that the chemical block for CFTR, CFTRinh-172, did not inhibit Cl/HCO3 exchange activity in wild-type pancreatic ducts. However, Wang et al. (50) and Shcheynikov et al. (46) both report siRNA against CFTR inhibits/alters Cl/HCO3 exchange activity. Our model suggests that the AE does not appear to contribute significantly to alterations in final saliva due to CFTR channel dysfunction.


Bicarbonate is an essential buffer for pH in saliva. However, it is difficult to measure directly. As such, concentrations of bicarbonate along the duct are unknown. However, because the salivary gland is thought of as a bicarbonate secreting epithelium, we felt the inclusion of bicarbonate as a variable in the final saliva composition was useful. Here we have provided one scenario based on the available data.

We have made an assumption that final saliva is electroneutral. Fortuitously, our model [HCO3]l is relatively constant as flow rate is increased, which looks quite similar to the mouse SMG displayed in Young and Schneyer (54) (Fig. 5B). In addition, we have found the luminal bicarbonate to be rather robust as we changed the AE from electroneutral, as described by Lee et al. (27), to electrogenic, as found in the parotid (23, 46).

Cellular and luminal pH.

We realize the our pH assumption may result in a coarse approximation of the duct behavior. pH is obviously an important component that cannot be disregarded entirely. In particular, understanding the transporters that control pH are likely key to understanding the mechanisms that affect caries progression. Unfortunately, we are not aware of sufficient data being available to constrain these fluxes in the model.

As a first attempt at a salivary duct model, we have assumed that the changes in cellular and luminal pH are not large enough to have a significant affect on the mechanisms represented in this model. This simplification helps us to reduce the complexity of the salivary duct and create a reasonable first approximation. In both compartments, pH is, in reality, affected by sodium proton exchangers (NHE; not represented in the model), which in turn may alter both the bicarbonate and the sodium concentrations in the model. Furthermore, it has been shown that the parallel activity of the apical NHE and the AE can mediate NaCl absorption in some epithelia (for example, Ref. 25). However, Catalán et al. (6) found that blocking the ENaC channel with lower concentrations of amiloride nearly abolished the NaCl reabsorption, suggesting that this channel is the primary Na+ influx mechanism. In addition, they found that the Nhe2−/− and Nhe3−/− mice did not have impaired NaCl reabsorption. This may indicate that NHE is not necessary for Na+ absorption or that the NHE is inhibited when the gland is stimulated, or it could be that other mechanisms are compensating for the loss of NHE in the Catalán et al. (6) knockout experiments. Obviously, this is an area that will require further investigation and modeling. Due to the success of our simulations capturing the wild-type behavior as well as the observed results associated with several impaired or knocked out transporters, we are optimistic that this simplified form will provide insight.

Flow rate.

We have used data from Nakamoto et al. (36) and Catalán et al. (6) to validate our model. Both report on ex vivo experiments where glands are stimulated by CCh, but the concentrations of CCh differ between the studies. Surprisingly, Catalán et al. (6) report a higher flow rate with 0.3 μM CCh than Nakamoto et al. (36) with 0.5 μM CCh. The flow rate differences may be explained by the changes in protocol since both the CCh concentration and the temperature during the experiments were changed; The Nakamoto et al. (36) data was collected at room temperature while the Catalán et al. (6) data was collected at 37° C. In addition, Catalán et al. (6) also apply IPR during stimulation. However, since it is well known that increasing CCh increases the flow rate, we have assumed a maximal flow rate of 8 μl/min for the Catalán et al. (6) data and a maximal flow rate of 20 μl/min for Nakamoto et al. (36) data. To convert from a whole gland flow rate to an acinar flow rate (the input flow, v0, for this model), the rate must be converted from μl/min to μl/s. As such, in addition to the glandular flow rate, we must have an estimate for the number of acinar cells per gland, the number of acinar cells per acinus, and the radius of the duct. For simplicity, we have assumed fixed values for all parameters except the flow rate and represent all stimulation by changes in the primary saliva and flow rate.

The final saliva composition is also dependent on the flow rate. At a low flow rate, the [K+]l is high and the [Na+]l and [Cl]l are low. As the flow rate is increased, the concentrations adjust: [K+]l decreases while [Na+]l and [Cl]l increase (Fig. 5B) before reaching a steady state similar to primary saliva. This is the same behavior reported by Mangos et al. (34) in the mouse parotid gland.

Strangely they do not observe the same behavior in the mouse SMG. Instead all three ion concentrations, [Na+]l, [Cl]l, and [K+]l, decrease before reaching a plateau. This difference could result from the choice of agonist, pilocarpine, used by Mangos et al. (34). Romanenko et al. (44) showed that the SMG and the parotid gland have significantly different responses to pilocarpine.

Future directions.

The duct has been represented here as a radially symmetric tube of length L. We do not represent any ductal branching or change in radius. It would be straightforward to include a changing radius, though Schneyer et al. (45) states that the luminal diameter of any part of the ductal network, excluding the excretory duct, averages 4 μm. Our aim has been to represent a basic secretory unit. In the future, we would like to continue this model to include branching, but for simplicity we have restricted our initial model to an unbranched duct.

In our model, we have assumed that pH does not affect these ion transporters. An important next step would be to explore the effects of modeling the pH along the duct. We expect pH to have an effect on bicarbonate due to buffering, but pH may also affect the membrane potential and the sodium concentration through the ENaC channel (42).


This work was supported by the National Institute of Dental Research Grant R01-DE-19245.


No conflicts of interest, financial or otherwise, are declared by the author(s).


Author contributions: K.P. designed model, carried out simulations, and drafted manuscript; M.A.C. designed model; J.E.M. performed analysis and interpretation of data; D.I.Y. peformed analysis and interpretation of data; E.J.C. designed model and revised draft for important intellectual content; J.S. designed model and revised draft for important intellectual content.


K. Patterson thanks Ted B. Begenisich at the University of Rochester and Ivo Siekmann at the University of Auckland for comments on an earlier version of this article.


  1. 1.
  2. 2.
  3. 3.
  4. 4.
  5. 5.
  6. 6.
  7. 7.
  8. 8.
  9. 9.
  10. 10.
  11. 11.
  12. 12.
  13. 14.
  14. 15.
  15. 16.
  16. 17.
  17. 18.
  18. 19.
  19. 20.
  20. 21.
  21. 22.
  22. 23.
  23. 24.
  24. 25.
  25. 26.
  26. 27.
  27. 28.
  28. 29.
  29. 30.
  30. 31.
  31. 32.
  32. 33.
  33. 34.
  34. 35.
  35. 36.
  36. 36a.
  37. 37.
  38. 38.
  39. 39.
  40. 40.
  41. 41.
  42. 42.
  43. 43.
  44. 44.
  45. 45.
  46. 46.
  47. 47.
  48. 48.
  49. 49.
  50. 50.
  51. 51.
  52. 52.
  53. 53.
  54. 54.
  55. 55.
  56. 56.
  57. 57.
  58. 58.
View Abstract