Quantification of Water Flux in Vesicular Systems

Water transport across lipid membranes is fundamental to all forms of life and plays a major role in health and disease. However, not only typical water facilitators like aquaporins facilitate water flux, but also transporters, ion channels or receptors represent potent water pathways. The efforts directed towards a mechanistic understanding of water conductivity determinants in transmembrane proteins, the development of water flow inhibitors, and the creation of biomimetic membranes with incorporated membrane proteins or artificial water channels depend on reliable and accurate ways of quantifying water permeabilities Pf. A conventional method is to subject vesicles to an osmotic gradient in a stopped-flow device: Fast recordings of scattered light intensity are converted into the time course of vesicle volume change. Even though an analytical solution accurately acquiring Pf from scattered light intensities exists, approximations potentially misjudging Pf by orders of magnitude are used. By means of computational and experimental data we point out that erroneous results such as that the single channel water permeability pf depends on the osmotic gradient are direct results of such approximations. Finally, we propose an empirical solution of which calculated permeability values closely match those calculated with the analytical solution in the relevant range of parameters.

Alterations in the single channel water permeability (p f ) of aquaporin's (AQPs) are directly related to human diseases 1,2 . Alongside with AQPs, different classes of proteins like transporters 3 , ion channels 4 and receptors 5 represent highly efficient water pathways. To get an idea about the impact of each of these transmembrane proteins on the water permeability of the respective tissue, accurate ways of quantifying unitary water permeability values p f of wild-type and mutant proteins are of particular interest and importance. Increased computer performance currently enables micro-second long MD simulations of transmembrane proteins. Despite their indisputable value to uncover molecular details of the substrate transport process, in silico data lack credibility without experimental validation. Furthermore, accurate p f values promote further development of more precise in silico algorithms and models. Fresh water scarcity boosts the field of biomimetic membranes 6 which seeks to generate highly selective and efficient lipid or polymer based filter membranes. Selectivity and permeability can be tuned by the incorporation of aquaporins (AQPs) 7-10 and supra-11 or single-molecular 12 artificial water channels. Carbon nanotubes are potential candidates as simulations already predicted high rates of waterflux [13][14][15] , which is only awaiting experimental verification 16 . However the performance of these biological or synthetic water channels as well as the membrane matrix itself, which is critical for channel stability and functionality 6,17 , can only be tested by accurately determining membrane water permeabilities P f . Different methods 18 exist, which are capable of retrieving information on water flux through lipid or polymer based membranes and incorporated transmembrane proteins or artificial water channels. On the one hand, in vivo techniques are available that can be divided into single cell swelling experiments with e.g. oocytes [19][20][21] and transcellular flux through polarized cells grown on porous support 22 . While oocyte swelling experiments maybe limited by the unfolding of the oolemma 3 , measurements of transcellular flux through polarized cells offer the unique advantage of online inhibition and channel manipulation 22,23 . On the other hand, in vitro systems like planar lipid bilayers 24 or lipid vesicles 25 exist. Dilution or up-concentration of reporter ions in stagnant water layers close to planar lipid bilayers due to an osmotic water flux can be measured by ion sensitive microelectrodes [26][27][28] . This so called scanning electrochemical microscopy technique is time consuming and experimentally challenging, but it allows simultaneous counting of embedded electrically active channels [29][30][31] . However, measurements on transmembrane proteins have notably different success rates, since some proteins hardly incorporate into these planar lipid bilayers for unknown reasons. Therefore, experiments are most commonly performed with giant unilamellar vesicles (GUV) 32 that are tens of μm in diameter or large unilamellar vesicles (LUV) with a typical size of about 100-120 nm. In contrast to GUV's dimensions, LUV's size cannot be determined with conventional light microscopy. A suitable signal to assess time driven size changes of liposomes is fluorescence self-quenching 33 , absorption/turbidity 34 and scattering intensity 35 . Due to superior sensitivity and ease of use, the latter method became the standard method implemented in today's stopped-flow experiments. The scattering intensity depends on the vesicle's size and the refractive index. After a hyperosmotic shock, the concentration of entrapped osmolytes increases due to a decrease in vesicle volume. This increase in inner osmolarity increases the vesicle's refractive index which elevates their light scattering ability. At the same time, the loss of vesicle size, reducing its scattering ability, counteracts this effect. To relate vesicle volume to scattering intensity, empirical approximations ranged from double logarithmic 25 , over quadratic 36 to simple linear 37 relations. We found that the computed solution of the scattering intensity in dependence of vesicle volume is well described by a second degree polynomial function of the vesicular volume 38 .
In a stopped-flow device, vesicles in osmotic equilibrium are subjected to a hyperosmotic buffer with impermeable solutes of concentration c out at time zero. The vesicle shrinkage depends on the water permeability of the vesicular membrane P f : where V w , V 0 , A, c in,0 and c out are the molar volume of water, vesicle volume at time zero, constant surface area of the vesicle, the initial osmolyte concentration inside the vesicles, and osmolyte concentration in the external solution, respectively. We assume that the osmolarities change indirectly proportional to the vesicle volume. This is justified by an almost linear relation between concentrations and osmolarities for various substances in the relevant concentration range 39 . Several decades ago these differential-algebraic equations were solved numerically 36,40 to calculate V(t) for various P f values. Time constants τ of mono-exponential fits to the permeability dependent V(t) were matched to experimentally fitted τ values to obtain P f . However, in the meantime different approximations 25,41,42 were proposed to calculate P f directly from τ. These models only differ in the way the osmotic conditions in the vesicles interior and exterior are considered ( Fig. 1b-d). An analytical solution exists for Eq. (1) at hyperosmotic conditions 38 (Fig. 1a), but still approximations are applied. Even though an exponential function (Fig. 2, black line) is not an exact solution to Eq. (1), it resembles the analytical solution. However, these approximations based on exponential fits to the scattering data lead to systematic errors and false positive effects as described below. Therefore we question all published approximations of estimating P f from light scattering traces of vesicle shrinkage and search for a new relation based on the time constant of exponential fits to the scattering data. To do so, we compare the discrepancies between the different methods by combined experimental and computational approaches. Permeability values obtained from different approximations are analyzed in respect to the absolute concentration of osmolyte, the osmotic gradient, the actual permeability and the vesicle radius. Finally we propose a new relation based on τ values which closely resembles P f,analyt calculated with our analytical solution.

Results
The rate of vesicle shrinkage depends on P f (Eq. 1). This differential equation was analytically solved by Horner et al. 38 (Eq. 7) and is depicted in Fig. (2) (black line). It may be fitted by an exponential function (Fig. 2, orange line) even though the latter is not an exact solution to Eq. (1). Three different approximations ( Fig. 1b-d) connect P f with the exponential time constant τ (Eq. 6) of vesicle shrinkage. They only differ in the way c in,0 and c out are considered. To test these models for their deviation from the analytical solution depending on P f,analyt , we first overexpressed, purified and reconstituted the aqua-glyceroporin of E.Coli, GlpF, into large unilamellar vesicles (LUVs). We subjected these proteoliposomes to a hyperosmotic solution and followed the decrease in vesicle volume by light scattering (Fig. 3). The raw data are either fitted by the analytical solution (Eqs 7 and 8) considering two vesicle populations 38 or by a double exponential function where the slow time constant corresponded to bare lipid vesicles and the fast time constant to GlpF containing vesicles. By adjusting the protein to lipid ratio during reconstitution, we experimentally varied the water permeability calculated with the analytical solution, P f,analyt , between 5 and 175 µm/s. The deviation of P f,exp from P f,analyt depends crucially on the approximation used (Fig. 4). However, the relative deviation is independent of P f,analyt . Comparing these approximations (Fig. 1b-d) for the  relation between τ and P f,exp revealed that the arithmetic mean of the two approximations (b, c) closely resembles P f,analyt in the investigated range of parameters with a maximal error of 5% (Fig. 1, new): To computationally verify this result we computed the analytical solution V(t) for a series of P f values between 2 and 200 µm/s. The time constants of the exponential fits to the data and the three published approximations (Fig. 1b-d) as well as our newly found approximation (Fig. 1, new) were used to calculate the corresponding membrane permeabilities P f . These permeability values are consistent with the experimentally found discrepancy between the approximations based on exponential time constants and the analytical solution (Fig. 4).
Next, we studied the deviation of P f,exp for different inner and outer osmolarities. Therefore, we varied c in,0 (97 to 919 mOsm) and c out (133 to 1778mOsm) of control vesicles (bare lipid vesicles) resulting in a gradient factor G = c out /c in,0 between 1.2 and 8.3. These experiments ( Supplementary Fig. S1) revealed that the deviation of the approximate models from the analytical solution do not depend on the absolute inner osmolarity c in,0 , but only on the gradient factor G (Fig. 5). Computational analysis confirmed these results ( Fig. 5 and Supplementary Fig. S2). The deviation of permeabilies P f,exp obtained from exponential fits to the analytical solution P f,analyt for the osmotic shrinkage of vesicles with radius r 0 = 60 nm and a water permeability of 20 µm/s depends on the relative inner and outer osmolarities but not their absolute values (Fig. 5 and Supplementary Fig. S2). To ensure comparability with experimental data we normalized all permeabilities P f,exp with P f,analyt . By calculating G for the proteoliposomes shown in Figs (3) and (4) and plotting them into Fig. (5), we strengthen the argument that the relative deviations between P f,exp and P f,analyt are independent of P f,analyt . The normalized permeability values perfectly match our control measurements. Hence Fig. (5) illustrates that the accuracy of the different models (Fig. 1b-d) only depends on G and in the range of investigated ratios, P f,exp may be decreased or increased compared to P f,analyt by a factor of six or several orders of magnitude, respectively. Our new approximation (Eq. 3, Fig. 1, new) gave similar values to P f,analyt , independent of the experimental conditions. The simple dependence of the error on the osmotic conditions allows recalculating an erroneous P f obtained by one of the old models (Fig. 1b-d) with the corresponding correction factor depicted in Fig. (1).

Discussion
Direct p f measurements are subject to large inherent technical difficulties. They include effects of stagnant water layers in membrane vicinity and uncertainties in the actual channel density 43 . Additional errors are introduced by using approximations (Fig. 1b-d) to calculate P f from the time constant τ of a single exponential fit to the scattering data. For example model d depicted in Fig. (1) produces an error that depends on the osmotic gradient 44,45 (see Fig. 5d), which can neither be explained by a structural resistance of LUVs to volume changes 33 , nor by an intra-vesicular unstirred layer 35,38 . Alternative methods, like scanning electrochemical microscopy 29 and fluorescence self-quenching experiments combined with a numerical solution 33 , confirm that P f is independent of the osmotic gradient. Still this model has widely been used to report for example water permeabilities for KcsA 4 , several AQPs 8,17,42,46-52 and artificial water channels 11,12 . Models b and c shown in Fig. (1) exhibit a smaller relative error when changing the ratio c out to c in,0 . Figure (5) depicts that approximation b and c represent a systematic over-or underestimation of P f,analyt with an increasing relative osmotic gradient. To conclude, varying G from 1 to 10, models b, c and d produce a systematic relative error up to a factor of 2, 6 and several orders of magnitude, respectively. We further demonstrate by means of computational and experimental data that the deviations of P f,exp from P f,analyt do not depend on P f,analyt or c in,0 but solely on G, the ratio c out /c in,0 . This is evident as models b to d only differ in the way they consider c out and c in,0 .
The above discussed approximations tend to systematically deviate from P f,analyt , due to oversimplifications in the respective derivations. For example, model b linearly approximates V(t) based on the initial change in volume 25 . The elaborate calculation yields a similar result as a heuristic approach that assumes an exponential volume change with the final volume depending on the ratio c in,0 /c out : Comparing the derivative of Eq. (4) at time zero ,0 with the initial change in volume (Eq. 1 at t = 0), one obtains approximation b. Yet the best-fit global time constant does not necessarily give the best fit to the initial rate of volume change which ill-poses the assumption of a strict exponential shrinkage. But, the flatter the kinetics, the better this approximation is. An example is provided by experiments with small relative osmotic gradients (Fig. 5). Approximation d can also be obtained heuristically, but with the assumption that the volume drops exponentially to zero instead to the ratio c in,0 /c out . Especially at very low osmotic gradients or gradient factors close to 1 model d leads to a divergent solution since only for large gradients c in,0 /c out approaches zero. Approximation c results in rather reliable values for small relative osmotic gradients as its derivation relies on small volume changes 41 (Fig. 5). Generally P f is independent of the initial vesicle radius r 0 38 . Any vesicle population with a unimodal radius distribution can be fitted by considering only one mean vesicle diameter. Still under-or overestimation of r 0 leads to an error in P f , which linearly depends on the ratio of assumed to real vesicle radius (Supplementary Fig. S3). Hence, r 0 must be known with even higher accuracy to resolve small changes in P f . It is important to note that for a scattering technique it is mandatory to estimate the mean radius of the intensity distribution and not the mean radius of the volume or particle distribution. Finally, we established a new relation to calculate accurate P f values from the time constant of an exponential fit to scattered light of vesicle shrinkage. Our new model (Fig. 1, new) yields an osmotic gradient independent permeability value consistent with the analytical solution of Eq. (1), which was already published in 2015. However, the analytical solution was only applied in our work group so far to calculate p f values of AQP1, GlpF, AqpZ, KcsA 38 , AQP4 53 and hSGLT1 3 . We attribute this to the fact that the Lambert function as part of the analytical solution is not implemented in common tools for stopped-flow data analysis. This survey enables the community to obtain correct permeability values from exponential fits or recalculate previous erroneous results.

Protein Overexpression, Purification and Reconstitution. The sequence of GlpF (Escherichia coli)
is inserted in a pTrc plasmid, which is transformed into C43 (DE3) cells. The cells are grown in LB overnight, diluted 40-fold and grown until OD 0.6. Expression is induced by 1 mM IPTG for 3 hours. Cells are harvested and pellets frozen at −80 °C. Purification and Atto488 labeling is performed as we described elsewhere 38,54,55 . In brief, cell pellets are lysed and the pelleted cell fraction is solubilized in a detergent buffer. After the removal of the insoluble material by ultracentrifugation, the supernatant is further purified using affinity (Ni 2+ -column). GlpF is reconstituted into proteoliposomes (PLs) with small modifications as previously described 38,54,55 . E.coli polar lipids (PLE, Avanti Polar Lipids) doped with 0.004 m% Atto633PPE are dried on a rotary evaporator. The thin lipid film is rehydrated in Reko buffer (100 mM NaCl, 20 mM MOPS, 1.4% OG, pH 7.4) to attain a final lipid concentration of 20 mg/ml. Subsequent to bath sonication, the clear suspension is incubated with equal amounts of protein diluted in Reko buffer at room temperature for an hour. With stepwise addition of Biobeads SM-2 (Bio-Rad), we remove the detergent within 36 hours. PLs are harvested by ultracentrifugation. The resuspended vesicles are centrifuged to remove aggregates and put through 21 extrusion cycles stacked with two polycarbonate filters with 100-nm pore sizes using a mini-extruder from Avanti Polar Lipids. This results in a unimodal radius distribution as seen from dynamic light scattering measurements. In addition protein reconstitution results in a fraction of bare lipid vesicles and a fraction of proteoliposomes with a heterogeneous but very sharp distribution of proteins between the lipid vesicles 38 . Control vesicles are treated similarly. All samples are assayed without delay.
Bare Lipid Vesicle Preparation. Large unilamellar vesicles (LUVs) were prepared from an E.coli polar lipid extract (PLE) mixture in chloroform as described elsewhere 56 . In brief, PLE was dried on a rotary evaporator, hydrated in a solution containing 100 mM NaCl and 20 mM MOPS buffered at pH 7.4, and extruded through 100 nm polycarbonate filters as described above. The final stock solution subsequently contained 10 mg/ml lipids. For experiments with different c in,0 the NaCl concentration was varied accordingly.

Stopped-Flow & Data Analysis.
PLs and LUVs are subjected to a hyperosmotic solution in a stopped-flow apparatus (SFM-300, Bio-Logic, Claix, France) at 4 °C. As previously described 3,4,38,53 , we monitor the intensity of scattered light at 90° at a wavelength of 546 nm. To calculate water permeability values from light scattering we use our recently found analytical solution 38 and three common approximations 25,41,42 as well as our new approximation all based on single-exponential functions with a time constant τ where depending on different models Π is equal to c out −1 , (c out -c in,0 ) −1 , c in,0 ·c out −2 or (c in,0 + c out )/(2· c out 2 ). The analytical solution for Eq. (1) at hyperosmotic conditions can be written as 38 :  2 However, in the relevant range of r 0 (~30 to ~100 nm) Eq. (6) can be approximated with a linear dependence of V(t) on I(t) with acceptable accuracy (Fig. S1 38 ). In combination with our new approximation (Eq. 3, Fig. 1, new) this results in P f values deviating only a few percent from P f,analyt (Fig. 4, cyan dots). In addition we take into account the fraction of bare lipid vesicles in a proteoliposome preparation as described elsewhere 38,53 . Osmolarities of c in,0 and c out were routinely checked using a vapor pressure osmometer 5500 from (Wescor, EliTechGroup, Utah). numerical V(t) values of the analytical solution (Eq. 7) to these arguments and returns the time constant τ of an exponential fit to the numerical values. An exponential fit to the analytical solution with initial and final volume as fit parameters violates the initial volume change (Eq. 1) as the fit overestimates V(t) at time zero as seen from the initial slopes in Fig. (2) However, we are only interested in the kinetic parameter τ, which is subsequently used in conjunction with Eq. (6) to derive approximate P f,exp values. These purely computationally derived P f,exp values are subsequently compared to P f,exp values obtained by the time constants of exponential fits to the scattering raw data, which comprises a linear dependence of V(t) on I(t). If not varied for systematic comparison, P f,analyt and r 0 are set to 20 µm/s and 60 nm. If not stated differently c in,0 is set to 200 mOsm and c out to 300 mOsm.

Computational Analysis.
Availability of materials and data. The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.