Increased apical Na+ permeability in cystic fibrosis is supported by a quantitative model of epithelial ion transport

Cystic fibrosis (CF) is caused by mutations in the cystic fibrosis transmembrane conductance regulator (CFTR) gene, which encodes an anion channel. In the human lung CFTR loss causes abnormal ion transport across airway epithelial cells. As a result CF individuals produce thick mucus, suffer persistent bacterial infections and have a much reduced life expectancy. Trans-epithelial potential difference (Vt) measurements are routinely carried out on nasal epithelia of CF patients in the clinic. CF epithelia exhibit a hyperpolarised basal Vt and a larger Vt change in response to amiloride (a blocker of the epithelial Na+ channel, ENaC). Are these altered bioelectric properties solely a result of electrical coupling between the ENaC and CFTR currents, or are they due to an increased ENaC permeability associated with CFTR loss? To examine these issues we have developed a quantitative mathematical model of human nasal epithelial ion transport. We find that while the loss of CFTR permeability hyperpolarises Vt and also increases amiloride-sensitive Vt, these effects are too small to account for the magnitude of change observed in CF epithelia. Instead, a parallel increase in ENaC permeability is required to adequately fit observed experimental data. Our study provides quantitative predictions for the complex relationships between ionic permeabilities and nasal Vt, giving insights into the physiology of CF disease that have important implications for CF therapy.


GHK flux equation
is the permeability of membrane x per unit area to ion n, is the electric PD across membrane x, F is the Faraday constant, R the universal gas constant, T is the temperature in Kelvin, and is the valence of the ion in question.
is the thermodynamic activity of n in the cell ( is extracellular activity), which is related to its chemical concentration via , where is the activity coefficient. The current through ion channels in the apical and basolateral membranes can be determined quantitatively by substituting the appropriate values of , , and into equation S1-S4. Paracellular flux of Na + ( ), K + ( ), Cl -( ) and gluconate ( ) can also be modelled using the GHK formulism. For serosal to luminal paracellular current (flow of positive ions), the appropriate electrical driving force is -, and the relevant concentrations are and , which are used in equation S5-S8.

NKCC co-transporter model
Full description of kinetic model and parameter values used is given in (Benjamin & Johnson, 1997).
We describe it briefly here.
The basolateral flux due to co-transport is given by where (a free parameter in our model) is the number of co-transporters per unit area of the basolateral membrane, and is the turnover rate of a single co-transport protein, given by , and they are also taken from reference (Benjamin & Johnson).

Na + -K + pump model
Full description of sodium potassium pump model are available in reference (Smith & Crampin, 2004), we use their model description and parameter values. A brief description of the pump model is outlined here.
Flux from the sodium potassium pump, , is given by where is the turnover rate of an individual pump protein, and is the number of pump proteins per unit area of the basolateral membrane (free parameter in our model).
The turnover rate is given by where are the forward and backward rate constants of the reduced 4 stage pump cycle: Normalised concentrations: Voltage dependent dissociation constants: Free inorganic phosphate is given by: The term in the denominator of is given by a sum of permutations of the rate constants

Determining paracellular permeability from shunt resistance
In our model we assume that under resting conditions the paracellular pathway is only permeable to Na + , K + and Clions, therefore the total paracellular or "shunt" current will be given by Assuming the paracellular currents are accurately described by the GHK formulism (equations S5-S7), and that the cation/anion permeability ratio is , then in the limit of identical bathing solutions this expression reduces to where is the permeability of the pathway to cations.
Shunt resistance is related to shunt current via therefore Mean shunt resistance was measured (Willumsen & Boucher, 1989) to be and for non-CF and CF epithelia respectively, with solution composition , and . For a non-selective pathway , we find that If we assume a cation selective paracellular pathway with (Levin et al., 2006;Flynn et al., 2009) the values we find are All parameter estimation and Monte Carlo filtering analyses were carried out with the following paracellular transport configurations: (a) unchanged in CF and non-selective, (b) reduced in CF relative to non-CF and non-selective, (c) unchanged in CF and cation-selective, (d) reduced in CF relative to non-CF and cation selective.
For simulations where we assumed a cation selective , we also assumed that .
We estimated the ratio of 1.4 based on recreating the experiments of Coakley et al. who found a permeability ratio (Coakley et al., 2003).

S2 Model validation
We carried out model validation to assess the stability of the dynamical system to small perturbations of model variables. The perturbations were introduced for numerical validation only, and not with the purpose of simulating a particular physiological phenomenon. The effect of this perturbation on cellular variables such as concentrations and membrane potentials can be seen in Figure S1. The system remains stable during the forcing, and relaxes to the same steady state as it was initially in, after the input term is switched off.
(b) Sinusoidal input: we implemented a sinusoidal input into the ODE for moles of Na + , with See Figure S2 for plots of physiological quantities over time for this perturbation.

S3 Numerical parameter estimation
The problem of estimating model parameters from observed experimental data can be formulated formally as the following optimisation problem: (S19) subject to the conditions: is the uncertainty in the experimental data measured at .
The first parameter estimation problem we solved was to find which minimised the residuals between observed data from non-CF HNE cells and model output, for a combined data set from two experiments, "+amiloride " and " ". These in silico experiments can be simulated by setting up the model equations and initial conditions appropriately. Amiloride block of ENaC channels is simulated assuming in model equations but with initial conditions found by solving where . The assumption that this inhibitor only affects the apical Na+ transport pathway is commonly made in modelling studies (Horisberger, 2003;Falkenberg & Jakobsson, 2010;Garcia et al., 2013) and also is the basis for interpretation of amiloride sensitive measurements in the standard nasal PD test (Knowles et al., 1995).
For modelling how reducing luminal affects system kinetics, initial conditions are found which solve with , and model equations then are numerically integrated from this initial condition with now at 3mM and = 117mM (note we assume that a gluconate concentration is introduced to replace the osmolarity of the Clions which are removed in this experiment). In both cases the change in model variables due to the "+amiloride" or " " perturbation can be found by taking the difference between the initial conditions and the new steady state values, The data used along with the resulting model output is shown in Table S3.
The second, separate parameter estimation problem was to find which would minimise residuals between observed data from CF HNE cells and model predicted variable values, again for data from a combination of "+amiloride" or " ". The data used in this problem is listed in Table S4.
In order to implement this formal parameter estimation we created a multistart optimisation problem in Matlab. A sample of randomly generated sets of parameter values were chosen from a uniform distribution (where is the set of baseline parameter values described in Table 1 of the main text). These were used as start points for multiple runs of the minimisation algorithm. The estimation problem (equation S19) was then solved, using fmincon with the interior point algorithm, for each randomly chosen start point. The solution which gave the lowest objective function value, out of resulting sample of minima found by the solver, was considered to be the global minimum in the region of parameter space searched and hence the solution to our problem.
The multi-start optimisation problem was run in four different scenarios, to assess the influence of paracellular permeability ( ) and selectivity ratios ( and ) on the estimated values of and . Initially, we used our baseline estimate for and assumed the paracellular pathway was equally selective (results in Table S5(a)), then we assessed the influence of an increased shunt resistance in CF, again without selectivity of the paracellular pathway (Table   S5(b)). The effect of a selective paracellular permeability was then assessed with the baseline value of in both CF and non-CF simulations (Table S6(a)), and the effect of selectivity combined with different CF and non-CF (  Table S6 (b), we can see the CF estimate of is greater than the non-CF estimate, but this is much greater than the likely error in both individual estimates, which is in both cases.
Hence this increase is significant given the different sets of HNE cell data. Figure S3 shows a plot of the model kinetics given by the optimal parameter values we determined in the baseline case, estimated from CF HNE cell data. The results of solving S19 for data from non-CF HNE cells, can be seen in figure 2 of the main text for comparison.

S4 Determining feasible non-CF and CF parameter distributions Monte Carlo sampling of parameter values from baseline estimates
To generate a population of individual parameter sets (a unique set of parameter values ) representative of the region of parameter space around the baseline parameter set , individual were randomly chosen from a uniform distribution . Each element of was chosen from its own uniform distribution (e.g. from etc), and each element was chosen independently of all the other elements.
The factor of 5 is arbitrary, this range is an attempt to represent the variability in transport parameters likely to be found in different cultures and cells, by not overlooking any physiologically relevant area of parameter space.

Physiologically realistic steady state behaviour
Equations S13-S18 describe a set of coupled ordinary differential equations (ODE's) of the form , where is the vector of transport parameters as described previously, and We chose the upper and lower bounds on allowed cellular variables to be the 90th and 10th percentile respectively, of the distribution measured for each variable, wherever an appropriate distribution had been published. A distribution of steady state cellular variables has been published for in normal and CF HNE cells (Willumsen et al., 1989a(Willumsen et al., , 1989bWillumsen & Boucher, 1991b, 1991a. The upper and lower bounds on the allowed change in cellular variables due to a "+amiloride" or " " experiment was determined differently, as distributions of these quantities was not presented for measurements. Here data for initial and final variable values (say before and after amiloride addition) was published as the mean value of the quantity in question, , plus or minus the standard error . In order to calculate the relevant upper and lower bounds on , we used since for the variance in , when , and are the variances associated with measurements of respectively.

Physiologically realistic transport kinetics
The second stage of model verification was to find parameter sets which predicted realistic transient behaviour of physiological variables, as well as realistic homeostatic behaviour. To do this we simulated a number of biologically relevant in silico experiments, for which experimental data on membrane potential and intracellular concentration kinetics was available.
The two types of simulated experiment carried out were (i) pharmacological block of ENaC channels and (ii) changing luminal Clconcentration, and these were implemented by appropriate choice of initial conditions and model equations, as follows: i. Amiloride block of ENaC channels.
Initial conditions are chosen which give a steady state with a non zero ENaC permeability (i.e. ). With these initial conditions, the system is numerically integrated assuming the ENaC current is completely blocked by amiloride, that is we set in the model equations.
ii. Reducing [Cl -] in the luminal solution.
Initial conditions are chosen which give a steady state when .
Model equations are numerically integrated, from this initial condition, with now set to 3mM and a gluconate concentration replacing the lost in the luminal compartment . We assume gluconate cannot permeate through the cell membrane, but will diffuse along the paracellular pathway, creating a current .
The net effect of (i) is to block the apical Na + current (i.e. ), and the net effect of (ii) is to increase the driving force for Clsecretion across the apical membrane (luminal osmolarity remains the same.
For each parameter set which remained after the initial model verification, the two in silico experiments as described above were simulated by numerically integrating the set of ODE's S13-S18, with the steady state value for that particular parameterisation used as the initial condition.
Numerical integration of the system was carried out in Matlab using the ode15s function, for a time interval much greater than the relaxation period of the system (t=3600s), to allow it to reach a new steady state. The changes in physiological variables between old and new steady states, can be recorded and compared with experiment.
which predict that are not physiologically realistic can be discarded. The expected changes in steady state variables used as filters are listed in Table 2 in the main text.
Parameter sets which predict steady state and kinetic model behaviour in quantitative agreement with observed data from non-CF HNE cells are shown in Figure S4, and those which can satisfy the CF constraints in Figure S5.

Influence of paracellular permeability & selectivity on feasible parameter distributions
As well as performing the MC filtering analysis for our baseline estimate value of , we also repeated the analysis to assess whether increased shunt resistance in CF or selective paracellular permeability would significantly alter the outcomes. In line with our parameter estimation work, we separately investigated the effect of (i) reduced in CF HNE cells, (ii) differential ionic selectivity in and (iii) reduced in CF combined with selective paracellular transport. The results of these further MC filtering analysis can be seen in figures Figure S6 - Figure S8. Qualitatively the relative difference between CF and non-CF & distributions is not altered in either of these scenarios, although the estimates of values are decreased if we assume .

S5 Model sensitivity analysis
To quantitatively assess the relative influence of transport parameters on physiological variables of interest, we carried out a sensitivity analysis using the data resulting from the model verification process. The approach we take is similar to that of Taylor, Goaillard and Marder (Taylor et al., 2009) where the authors are determining how different conductances affect electrophysiological properties of an multi-compartmental cellular model of a neuron, and Sobie (Sobie, 2009) Once the population of parameter sets is z-scored, the regression model was fit using Matlab function regstats with option 'quadratic'. The regression coefficients determine the strength of the linear correlation between parameter and output . By comparing the magnitude of the linear regression coefficients we have a way of objectively determining the relative influence that transport parameters have on the physiological variable of interest. In our study we carry out sensitivity analyses to determine the relationships between transport parameters and basal , , and . The results of these analyses are listed in Table S9.    Table S3: (a) Experimental data used for estimation of non-CF transport parameters. Concentrations are given in mM, membrane potentials in mV, and time in seconds. Equivalent data predicted by optimal parameter set in the case of baseline and nonselective paracellular transport are shown in (b).

Sources of observed data:
for n=1, 2 : Data found by digitizing plots in Figure 8 of reference (Willumsen & Boucher, 1991a), using the initial steady state and final data points from each time course plot. for n = 3, 4: Table 4 in (Willumsen et al., 1989a)* for n = 5, 6: Table 3 in (Willumsen et al., 1989a) Table S4: (a) Experimental data used for estimation of CF transport parameters. (b) Equivalent data predicted by optimal parameter set in the case of baseline and non-selective paracellular transport.
Sources of observed data are: for n = 1, 2: Table 3 in (Willumsen et al., 1989b)* for n = 3, 4: Data found by digitizing plots in Figure 7 of reference (Willumsen & Boucher, 1991b), using the initial steady state and final data points from each time course plot.
for n = 5, 6: Table 2 in (Willumsen et al., 1989b)* * Assuming system is at steady state when measurements are made at .
Note since measurements of and from two separate "+amiloride" experiments are included in the objective function, we weight the residual error from each of these data points by ½, in order to give the same weight to these measurements as to and , for which only 1 set of measurements is available.  Table S5: Results of parameter estimation via minimisation of residual error between model predictions and experimental data given in tables S3 and S4. (a) Estimates of transport parameters found when is assumed to be fixed at baseline estimate in both CF and non-CF simulations, and is assumed to be the same for all ions. (b) Parameter estimates when paracellular permeability is assumed to be decreased in CF (Willumsen & Boucher, 1989) Table S9: Linear regression coefficients ( ) found by fitting multiple regression between z-scored transport parameters ( ) and physiological variables ( ), for use in sensitivity analysis. Figure S1: Model validation with constant input. The rate of change of cell volume is altered so that there is a constant forcing input between and . Other cellular quantities change accordingly, and the system reaches a new steady state , before relaxes back to the initial steady state once the constant input is switched off. Note reversal potentials are normalised to their initial steady state value ( : +44.0 mV, : -23.3 mV, : -85.5 mV) Figure S2: Model validation with sinusoidal input. The ODE for rate of change of Na + moles in the cell is perturbed with a sinusoidally varying input term between and . All other cellular variables are affected to some extent by this forcing of and vary with the same periodicity as the input, before returning to the same steady state as they were initially in, after the forcing term is switched off. Note reversal potentials are normalised as in figure S1.  Table S5(a)).  ), for the distribution of non-CF values found. For a given low Cl response, hyperpolarisation or depolarisation of basal is possible, but there appears to be a limit on the hyperpolarisation possible (suggested by dashed line), which is significantly less than the magnitude observed in CF.