Unidirectional Flux Balance of Monovalent Ions in Cells with Na/Na and Li/Na Exchange: Experimental and Computational Studies on Lymphoid U937 Cells

Monovalent ion traffic across the cell membrane occurs via various pathways. Evaluation of individual fluxes in whole cell is hampered by their strong interdependence. This difficulty can be overcome by computational analysis of the whole cell flux balance. However, the previous computational studies disregarded ion movement of the self-exchange type. We have taken this exchange into account. The developed software allows determination of unidirectional fluxes of all monovalent ions via the major pathways both under the balanced state and during transient processes. We show how the problem of finding the rate coefficients can be solved by measurement of monovalent ion concentrations and some of the fluxes. Interdependence of fluxes due to the mandatory conditions of electroneutrality and osmotic balance and due to specific effects can be discriminated, enabling one to identify specific changes in ion transfer machinery under varied conditions. To test the effectiveness of the developed approach we made use of the fact that Li/Na exchange is known to be an analogue of the coupled Na/Na exchange. Thus, we compared the predicted and experimental data obtained on U937 cells under varied Li+ concentrations and following inhibition of the sodium pump with ouabain. We found that the coupled Na/Na exchange in U937 cells comprises a significant portion of the entire Na+ turnover. The data showed that the loading of the sodium pump by Li/Na exchange involved in the secondary active Li+ transport at 1–10 mM external Li+ is small. This result may be extrapolated to similar Li+ and Na+ flux relationships in erythrocytes and other cells in patients treated with Li+ in therapeutic doses. The developed computational approach is applicable for studying various cells and can be useful in education for demonstrating the effects of individual transporters and channels on ion gradients, cell water content and membrane potential.


Introduction
The concept of the pump-leak flux balance as the basis of monovalent ion gradients at the animal cell membrane is universally accepted. A number of various transporters and channels are involved in continuous ion traffic across the membrane and many of them are capable of passing ions both inward and outward. However, discrimination between fluxes via specific ways is not a trivial problem because any macroscopic ion transfer is accompanied by disturbance of cell water and electrical balance. Fluxes of different ions and via different routes appear to be interdependent due to the mandatory conditions of electroneutrality and osmotic balance. In addition, some transporters operate as a co-or counter-transporters. Calculation of the overall flux balance and prediction of its dependence on specific properties of transporters and channels can be performed by the computational solution of a set of nonlinear differential equations [1][2][3][4][5][6][7][8][9]. However, there are no sufficiently simple computational tools for solution of real cell physiology problems. Most experimentalists continue to neglect computational approaches because many parameters are required for modeling, whose evaluation is difficult and unreliable. Not all types of the monovalent ion movement across the cell membrane are accounted for in the available models. Ion traffic of the self-exchange type that comprises a significant portion of Na + and Cl − fluxes across the membrane remained beyond the scope of previous models. We aimed to develop relatively simple software for analyzing the effects of various transporters and channels on cell water-volume, membrane potential and related cell properties under various conditions suitable for researchers with limited programming expertise.
Our approach was developed initially for studying Li + transport. Li + is the closest physiological analogue of Na + and the Li/Na exchange is the closest analogue of the Na/Na exchange. Li + is a poor substrate for the Na/K-ATPase pump but it passes through the same channels as Na + , and their gradients on the cell membrane are comparable. For example, the ratio of balanced intracellular to extracellular concentrations in U937 cells is 0.82-0.96 for Li + and 0.28-0.30 for Na + , whereas for K + it is 30-32 [10]. It is the Li/Na exchange that mediates secondary active Li + transport out of cells [10][11][12][13]. The mechanism of Li + transport and of Li/Na exchange, in particular, is important for a number of practical reasons: alteration of Li/Na exchange in erythrocytes accompanies widespread human pathologies (hypertension, diabetes, nephropathy etc.); Li + is used as a medication for treatment of neuropsychiatric disorders and testing renal clearance [10,[14][15][16][17][18].

Materials and Methods
U937 human histiocytic lymphoma cells were obtained from the Russian Cell Culture Collection (cat. number 160B2). The cells were cultured in RPMI 1640 medium (Biolot, Russia) with 10% fetal calf serum (HyClone, USA). Ouabain and dimethylamiloride (DMA) were purchased from Sigma-Aldrich (Germany), Percoll was from Pharmacia (Sweden) and the salts (all of analytical grade) were from Reachem (Russia). Intracellular cation content was determined by flame emission on a Perkin-Elmer AA 306 spectrophotometer, Cl − content by distribution of 36 Cl − (Isotope, Russia) and chemical external Cl − assay, cell water was evaluated by cell buoyant density in Percoll density gradient, and protein was measured by the Lowry method. The experimental methods used in this work have been described in detail earlier [10,[19][20][21][22][23][24]. Na + and Cl − fluxes were estimated using 22 Na + and 36 Cl -. For quantifying the pump-mediated Rb + influx cells were incubated in the medium with 2.5 mM RbCl with and without 0.1 mM ouabain for 10 min. To study ion efflux, cells were pre-equilibrated in the presence of ions for either 1.5 h (with 5 mM LiCl, 22 Na + or 36 Cl -) or 20 h (with 1 mM RbCl). Next, the cells were rapidly washed with 96 mM MgCl 2 and placed into fresh RPMI medium. The radioactivity of 36 Cl was analyzed using a Beckman LS 6500 liquid scintillation counter, and radioactivity of 22 Na + was counted using a Na-Iodide detector. The equilibration rate coefficients k were calculated according to y t = y 1 (1-exp (-kt)) for influx and y t = y 0 exp (-kt) for efflux, where y t is the content of ion-label at time t (5, 10, 20 min) and y 1 and y 0 are the final and initial content of the corresponding ion. Li + balanced distribution in U937 cells was studied at 1-10 mM LiCl added to the culture medium. To replace intracellular K + and Na + with Li + , cells were incubated in a medium containing 140 mM LiCl, 5.5 mM KCl, 0.42 mM CaCl 2 , 0.41 mM MgCl 2 , 10 mM Hepes, pH 7.4 was adjusted with LiOH. To minimize errors in determination of the Li/Na discrimination coefficients, the medium was analyzed simultaneously with the cell samples. The standard value of 140 mM for Na + concentration in RPMI medium was used in other calculations.
The results were analyzed using Student's t-test. Differences were considered significant at P<0.05.

Overview of the movement and distribution of monovalent ions in U937 cells
The kinetics of gain and loss of monovalent ion tracers is shown in Fig 1. The rate of equilibration is the highest for Cl − and Na + , intermediate for Li + , and the lowest for Rb + , which is a known analogue of K + . By contrast, with Na + , the rate coefficient for Li + exit is not reduced in the presence of ouabain (Fig 2). This agrees with the notion that the Na/K ATPase pump is not involved in Li + efflux out of cell. The efflux of Li + is reduced by about two-fold by DMA, an inhibitor of the Na/H exchanger. Intracellular Na + and Li + concentrations in U937 cells equilibrated with 1-10 mM external Li + are presented below. The rates of ion tracer equilibration, as in Fig 1, correspond to uptake and loss of about 15% of cell Na + and Cl -, 5% of Li + , and 1% of Rb + (K + ) content per min. The corresponding data in absolute units for the U937 subline 160B2 cells [10,19,22,24] were 30-40 μmol -1 min -1 per g protein for Na + and Cl − and 6-7 μmol -1 min -1 for K + .
Computational dissection of ion traffic via specific routes and determination of parameters characterizing the cell as a dynamic electrochemical system to the Na/H [6]. In this approach, the entire multiplicity of ion transfer systems is replaced by a reduced number of ion pathways defined thermodynamically. All the major pathways are subdivided into five subtypes by ion-driving force: ion channels, where the driving force is the transmembrane electrochemical potential difference for a single ion species; NKCC, NC and KC cotransporters, where the force is the sum of the electrochemical potential differences for all partners; and the Na/K ATPase pump, where ion movement against electrochemical gradient is energized by ATP hydrolysis. This allows characterization of intrinsic properties of each pathway by a single rate coefficient. The following abbreviations are used to designate ion transporters: NKCC (LKCC) indicates the known cotransporters of the SLC12 family carrying monovalent ions with stoichiometry 1Na + (Li + ):1K + :2Cland KC and NC (LC) stand for cotransporters with stoichiometry 1K + :1Cl-or 1Na + (Li + ):1Cl -. The latter can be represented by a single protein, the thiazide-sensitive Na-Cl cotransporter (SLC12 family), or by coordinated operation of the exchangers Na/H (SLC9) and Cl/HCO3 (SLC26). Monovalent ions are the main factor in maintaining electrical and water balance in animal cells. Because concentration of free Ca ++ and Mg ++ in the cytoplasm is low; their impact on electrical and osmotic balance, as well as on the transmembrane flux balance, is small, and these ions were not included in the calculations. Numerical integration of the flux equations was carried out as described in the Appendix. The variables and parameters used in the software are listed in Table 1. Some of the data necessary for the initial file DATAP.txt (Table 2) were obtained experimentally; these are the intra-and extracellular monovalent ion concentrations, extracellular concentration of membrane-impermeant osmolytes (B0), the ratio of ouabain-sensitive to ouabain-resistant Rb + (K + ) influx (OSOR) and the Na/K pump rate coefficient β, which is calculated from the measured [Na] i and ouabain-sensitive Rb + influx. The process of finding the rate coefficients characterizing the ion pathways-pna, pk, pl, pcl, inc, ikc, ilc, inkcc, ilkcc is discussed below.
The executable file LIP.exe generates a file RESP.txt which shows the evolution of the cell toward the balanced state, when the influx and efflux of every ion species become equal ( Table 3). The fluxes via the major pathways are given for the last state. In some cases the balanced state cannot be reached, e.g. when the conditions of macroscopic electroneutrality and osmotic balance are not fulfilled or the membrane potential falls outside the range +5 mV to -173 mV (computer reports RANGE LIMIT). One should keep in mind that the intracellular concentrations [Na]    The parameters na0, k0, l0, cl0 and B0 are extracellular ion concentrations (mM); na, k, l and cl are intracellular concentrations; kv is the ratio of the external to initial internal osmolarity; alpha and beta are the rate coefficients (min -1 ) for Li/K and Na/K pumps, respectively; gamma is the pump Na/K stoichiometric coefficient; pna, pk, pl, pcl are the channel permeability coefficients (min -1 ); inc, ikc, ilc, inkcc, ilkcc are the rate coefficients for the NC, KC, LC (mlÁμmol -1 Ámin -1 ) and for NKCC, LKCC (ml 3 Áμmol -3 Ámin -1 ) cotransporters; kp is the rate coefficient for Li/Na countertransport. The parameter hp is the number of integration steps between data output and corresponds to the real time of the transient process in min. If one wishes to follow the initial stages of the process, hp should be made small , and if the process is slow and one is interested in the balanced state, hp can be increased to 1000-5000. The data from a typical experiment with U937 cells are used in this example.
doi:10.1371/journal.pone.0153284.t002 compared with experimental data as a test that the parameters have been chosen correctly. The inverse problem (finding the parameters from experimental data) is much harder. There are mathematical methods [26] and software packages [27] for finding optimal parameters in similar cases. However, the use of automated procedures requires more mathematical sophistication than a simple trial and error approach, as well as consideration of specific properties of the system. Most importantly, the system of equations for describing cell ion balance has only a unique solution for a balanced state within the acceptable range of variables as can be verified by using our soft. It should be noted that the accuracy of parameter estimation depends on the dispersion of the experimental data and not on the numerical solution procedure, because the calculation process continues until the difference between the calculated and experimental data becomes less than the experimental error. In the case of U937 cells the parameters ikc, inkcc, ilkcc can be excluded for the following reasons. The inside Cl − electrochemical potential under the normal balanced state significantly differs from that outside (mucl = +33.7 mV, Table 3). Therefore, a secondary active transport of Cl − into the cells should be present. The driving force for Cl − ions involved in cotransport is determined by the sum of the electrochemical potential differences for all the partners. Under normal conditions, these sums for NC, KC and NKCC pathways in U937 cells are -48.2, +73.9 and +25.7 mV, respectively (Table 3, mun+mucl, muk +mucl, mun+muk+2mucl). Therefore, the net Cl − flux should be directed outward through KC and NKCC and inward through NC. As the integral Cl − transport, which creates the +33.7 mV net electrochemical potential difference for Cl -, is directed into the cell, the net fluxes through the KC and NKCC routes have to be negligibly small in comparison with the NC cotransport. The small NKCC and KC impact is confirmed by the fact that the effects of bumetanide and DIOA, used as inhibitors of NKCC and KC pathways, are small in the U937 cells compared to epithelial cells [20]. It is reasonable to accept that permeability of electroconductive channels for Li + and for Na + are similar, i.e. pl % pna [28]. The channel permeability coefficients pna, pk, pl, pcl were found by trial and error by seeking the same balanced state as was observed in the experiment. This is not a particularly time-consuming procedure if one takes into account that [K] i and [Na] i are affected mostly by the parameter ratios pna/β and pk/β, while [Cl] i and cell water content per unit of impermeant intracellular osmolyte A are mostly affected by inc/β and pcl/β [9]. Analysis of the present model shows that [Li] i depends mostly on ilc/β, pl/β and kp.
Computation shows that the balanced intracellular Li + concentration of 4.37 mM, as well as the other variables observed in U937 cells upon addition of 5 mM LiCl, could be obtained at many combinations of ilc and kp. To specify ilc and kp, kinetic data on Li + equilibration were used. Only a unique pair of ilc and kp (ilc = 0.00018, kp = 0.0002) corresponds to the real kinetics of Li + gain and exit. Curves 1 in Fig 3A and 3B show the results obtained at values of ilc and kp that fit the observed balanced state but not the transition kinetics, whereas curves 2 fit both.  Table 2.
doi:10.1371/journal.pone.0153284.g003 Table 3 imitates the table displayed when the calculation is finished. In addition to the initial parameters already contained in Table 2, the parameters and variables in Table 3 include the following: current intracellular ion concentrations na, k, l, cl (mM); membrane potential U (mV); cell water content V/A x 1000 (ml per mmol of A); transmembrane electrochemical potential differences mun, muk, mul, mucl (mV) for Na + , K + , Li + and Cl -, respectively, calculated as mun = 26.7Áln(na/na0)+U, muk = 26.7Áln(k/k0)+U, mul = 26.7Áln(l/l0)+U, mucl = 26.7Áln(cl/cl0)-U; the derivatives prna, prk, prl, prcl, characterizing the rate of change in [Na] i , [K] i , [Li] i and [Cl] i ; net fluxes of Na + , K + , Li + and Cl − via separate ion pathways (μmolÁml -1 Ámin -1 ); A/V (mM) and the mean valence of membrane-impermeant intracellular osmolytes, z. The displayed values of fluxes, as well as OSOR, the ratio of ouabain-sensitive to ouabain-resistant K + influx, are calculated only in the bottom line, i.e. for the latest time point. The fluxes at any time can be found by varying hp. The values of fluxes per g of cell protein can be found by multiplying the displayed values by cell water content per g of cell protein, which is close to 6.5 mlÁg -1 in "standard" U937 cells under normal conditions. Transmembrane electrochemical potential differences mun, muk, mul, mucl expressed in mV can be converted to J/ mol as 1 mV = 10 -3 F -1 (J/mol).

Validation of the model by studying Li + balanced distribution at 1-10 mM LiCl
We showed earlier that the presence of 3-5 mM Li + in the medium for 12-24 h has no effect on cell proliferation, cell morphology or the balance of monovalent ions in U937 cells [10]. Table 4 presents experimental data on Li + balanced distribution in U937 cells incubated in RPMI medium with addition of 1-10 mM LiCl. The ratio of balanced intracellular to extracellular concentrations were 0.82-0.96 for Li + , 0.28-0.30 for Na + , and 30-32 for K + . These results confirm that the distributions of Li + and Na + are qualitatively similar and, at the same time, quite different from the distribution of K + . Quantitatively, there is still a three-fold difference between the distribution coefficients for Li + and Na + when the Na/K ATPase pump is active. The concentration of Li + in the medium slightly influences the [Li] i /[Li] o and [Na] i /[Na] o ratios. The coefficient of Li / Na discrimination, c d , better characterizes the differences between the distributions of Li + and Na + because it is not affected by possible errors in determination of the cell water content; c d increases slightly but significantly, when [Li] o increases from 1 to 10 mM. The value of c d % 3 in U937 cells is within the range found for human erythrocytes at very low (0.1-1.0 μM) Li + concentrations in the plasma [29], for patients regularly taking Li + when its concentration in the plasma of about 0.5-1 mM [14] and for the muscle tissue of  [30][31][32]. Evidently, the Li + -Na + discrimination is determined by very conservative factors.  Table 6B. The general conclusion is that the model based on constant parameters predicts the general trends in Li + distribution reasonably well; nevertheless, further improvements achieved by fine-tuning of the parameters suggest that they vary slightly with [Li] o . We also conclude that the maintenance of the balanced [Li] i % 0.9[Li] o , as observed in the experiments, can indeed be provided by an outward net Li + flux due to Li/Na coupled transfer, LN. If LC were primarily responsible for Li + maintenance then the net Li + flux would have the opposite sign, i.e. directed inward, down the Li + electrochemical gradient. The Li + and Na + net fluxes via the LN pathway are smaller than the active Na + flux through the pump by a factor of 18 (Table 3). Thus, only a minor fraction of the Na/K ATPase energy consumption in this example falls on the secondary active transport of Li + . The inward Li + net flux of 0.0839 down the gradient splits up between the LC pathway (62%) and the channels (38%). The dynamics of Li + gain and K + exit following a complete external Na + for Li + substitution resembles the case when the pump is inhibited by ouabain (Fig 4A and 4B). As soon as intracellular Na + is lost, the Na/K-ATPase pump stops (Fig 4C). Therefore, in cells with low internal Na + , 80% of K + is replaced with Li + by 4 h. In spite of this, the cells maintain their size and water content at nearly the same level as in the control [10]. The normal K/Na ratio can be restored by returning cells to RPMI medium ( Fig 4F). The data on transients caused by substitution of Na + for Li + , as well as by application of ouabain, support the validity of the model in general. A good match between the observed and predicted dynamics of ion disturbance following ouabain treatment has been reported previously [9]; the more detailed analysis presented here shows that further improvement can be achieved by assuming a small (2.5%) residual activity of the pump (the dotted line in Fig 4A). The match between the theory and the Li + experiment is less perfect (Fig 4B); nevertheless, the reciprocal equivalent changes in intracellular Li + and K + concentration are predicted sufficiently well. The fit is improved if one assumes a Li/K pump activity with the rate coefficient 0.004, i.e. at Whole-Cell Unidirectional Ion Fluxes 10% of the normal Na/K pump activity (Fig 4D). The other way to improve the fit is to decrease pl by a factor 2.3 and to increase ilc by a factor 1.3 (Fig 4E). Varying the other parameters does not improve the agreement between the measured and computed data. The best fit of the transient kinetics following reverse replacement of external Li + for Na + is achieved by assuming activation of the Li/K and Na/K pumps (alpha = 0.04 and beta = 0.08, Fig 4F). Comparison of the calculated and experimental data is the way to find out how ion-transporting mechanisms are altered under variable conditions. Note that the rate of the transient process associated with the macroscopic changes in cell ion content may be much lower than the rate of the tracer 22 Na + or 36 Cl − equilibration alone due to rapid Na/ 22 Na or Cl/ 36  Li/Na exchange with and without specific Li/Na countertransporter The Li/Na exchange is usually evaluated by comparison of Li + efflux from cells loaded with Li + into the media with or without Na + . Computational modeling allows separating the effect of the specific Li/Na countertransport and the effect of the mandatory electroneutrality of ion transfer. Fig 5 shows the changes in Li + , Na + , K + and Cl − concentrations in U937 cells preloaded with Li + during recovery of the ion balance in Li + -free media with and without Na + . Experimental results are compared to the calculations based on the models with or without Li/Na countertransporter. When cells are placed in RPMI half of the Li + exit during the first 15 min is balanced by Na + gain, whereas the other half is balanced by K + gain. Later a decrease in both in Li + and Na + concentrations takes place, which is balanced by accumulation of K + . Both experiment and computation show that Li + exit into a Na + -containing medium during the first stage is balanced to a large extent by Na + gain even in the absence of Li/Na countertransport (kp = 0) while later it is balanced by K + gain. By contrast, Li + exit into a Na + -free medium is accompanied by K + exit and the loss of both cations is compensated by an equivalent Cl − loss. The difference in the recovery kinetics in these cases is associated with different changes in the membrane potential and water balance (Fig 5, graphs on the right). Li + loss and Na + gain are enhanced in the model with a Li/Na countertransporter (kp = 0.0002) but only during the first stage, when the outward Li + net flux is counterbalanced by an inward Na + net flux. These complicated relationships between ion fluxes should be taken into account when Li/Na substitution is used to test for alteration in a specific Li/Na countertransporter.

Full list of major unidirectional fluxes across the plasma membrane in U937 cells under the balanced state. Coupled Na/Na and Cl/Cl exchange
Once validated, the model can help estimate every unidirectional and "electroneutral" coupled flux. Table 7 shows unidirectional fluxes through the pump, the channels and the NC, LC, and LN transporters in U937 cells balanced in RPMI+5 mM LiCl. Summation of all unidirectional effluxes or influxes gives the overall flux representing continuous traffic of ions into and out of cell, the "turnover" flux. This flux determines the rate of equilibration of ions-tracers. The computed Li + entire flux of 0.1849 μmolÁmin -1 Á (ml cell water) -1 at [Li] i = 4.4 mM corresponds sufficiently well with the observed equilibration rate coefficient of 0.042 min -1 at [Li] o = 5 mM (Fig 3A and 3B). The sum of the computed unidirectional Na + influxes through the channels, NC and LN (or the sum of effluxes through the pump, channels, NC, and LN) is significantly less than the observed turnover flux, which is about 5.8 μmolÁmin -1 Á(ml cell water) -1 (Fig 3C). Therefore, there has to be one more unaccounted for Na + flux. It is the equivalent Na/Na exchange, which has no impact on ion concentrations or on membrane potential but adds significantly to the turnover flux.
In order to find the rate coefficient for the Na/Na coupled exchange, the 22 Na gain or exit was computed by using Eq (9) for Li + with parameters similar to those of the Na + carrier. A Li + -free media was assumed in this case and the Li/Na coupled exchange had the meaning of 22 Na/Na exchange, while the parameter kpNa fitting the experimental 22 Na curve represented the rate coefficient for the 22 Na/Na exchange. The curves 1 and 2 in Fig 3C show the computed  kinetics of 22 Na exit for kpNa = 0 (no Na/Na exchange) and kpNa = 0.0008, respectively. As the cells cannot distinguish between 22 Na and Na, the rate coefficients kp for both isotopes should be equal. Therefore, the macroscopic flux related to the Na/Na coupled exchange in the normal Li + -free media according to Eq (10) is J NN = kpNaÁ[Na] o Á[Na] i = 0.0008Á140Á37 = 4.14 μmolÁmin -1 Á(ml cell water) -1 . Interestingly, the rate coefficient kpNa is four times larger than the coefficient kp for the Li/Na exchange. The Na/Na coupled exchange appears to be very large in comparison with the fluxes via other pathways (Fig 6). It exceeds the pump-mediated Na + efflux by a factor of 2.8, the channel Na + influx by a factor of 4, and accounts for 71% of the total Na + flux. Therefore Na + influx involved in secondary active Li + efflux relates to active Na + efflux by the Na/K-ATPase pump as 0.1225 to 1.5051 (at 5 mM external Li + ). We may conclude that Li + active transport in other human cells at therapeutic Li + concentration does not load significantly the sodium pump. A major mechanism for Li + efflux appears to be the Li/Na exchange whereas for Li + influx it is the LC cotransport.
It should be noted that the current concept of ion transfer via channels and transporters assumes that there are forward and backward fluxes in every pathway, including electroconductive channels. Evidently, we should distinguish the genuine coupled exchange and "phenomenological" self-exchange which is determined by reversibility of transport process irrespective of its molecular mechanism. In the used formalism, the difference between direct and reverse fluxes is much greater for channels than for transporters. There is much evidence that Li/Na exchange may occur via various pathways where sodium-linked transport takes place, e.g. NHE [11,33,34], NKCC [13], NCC [35], Na + coupled phosphate transporter [12,36,37] etc. Similar transporters may be responsible also for specific Na/Na exchange.
The K + flux balance in the minimal cell model is relatively simple, because the impact of fluxes via NKCC and KC cotransporters in U937 cells should be small due to small forces driving the ions by these mechanisms. The sum of K + fluxes via the pump and channels fully corresponds to the measured total K + turnover flux and there is no need to hypothesize a K/K equivalent exchange. The coupled Cl/Cl exchange has been described in various cell types [38][39][40][41]. The approach developed in the present study allows computation of fluxes due to Cl/Cl equivalent exchange along with other unidirectional Cl − fluxes although. However, the executable file should be slightly modified. That would go beyond the scope of the present paper; nevertheless, we have included in Fig 6 our unpublished experimental and computational results on Cl/Cl exchange in U937 cells. It accounts for 91% of the total Cl − flux in U937 cells under considered conditions. Parameters used in the computation as in Table 3; the fluxes in μmolÁmin -1 Á (ml cell water) -1 doi:10.1371/journal.pone.0153284.t007 Whole-Cell Unidirectional Ion Fluxes

Conclusion
Monovalent ion fluxes across the plasma membrane are not independent, but linked through the factors that can be classified into three groups. First, fluxes are constrained by the conditions of osmotic balance and electroneutrality. Second, they are linked through coupled ion transfer, such as occurs in the Na/K-ATPase pump and various co-and counter-transporters. Third, transporters can share common dependence on the same factors, such as pH (as is the case with Na/H and Cl/HCO 3 countertransporters) or on mediators of intracellular signaling: such dependence can manifest itself in variable rate coefficients. The factors of the first two groups can be accounted for by mathematical modeling and numerical integration of the flux balance equations. This has been demonstrated in a series of fundamental works from several research groups [3,5,6,9,25,[42][43][44]. Unfortunately, most experimentalists continue to neglect computational approaches because modeling includes many parameters whose evaluation is difficult and unreliable and because of the scarcity of software available for researchers with limited programming expertise. We thus attempted to develop such a computational tool. Following the previous authors [24][25][26]28] we used a mathematical model of ion traffic across the cell membrane that differentiates ion transfer systems by a thermodynamic characteristic, i.e. by the driving force acting on ions. Each pathway can be characterized by a single rate coefficient that can be rigorously determined in the experiment and used as an intrinsic physiological parameter of a resting cell. These resting values can then be applied to a system where the balance has been disturbed to see if the model continues to accurately describe cell behavior. If it does, that can be taken as an indication that the parameters have not changed. We applied such an approach to an experiment where all the external sodium was replaced with lithium. Our study shows that the analysis of the balanced state allows a rather accurate prediction of the following: (1) changes in the internal Li + concentration under variation of external Li + within 1-10 mM; (2) kinetics of the transients caused by blocking the pump with ouabain. The latter result exceeded expectations because ouabain treatment causes profound changes in the cell ion content, and one could expect that the properties of transporters and channels would be affected as well. Nevertheless, the parameters remained unchanged. Not so good was the match between the computed and measured kinetics after full replacement of external Na + for Li + . However, the general time course of changes in ion and water content and the reciprocal changes in intracellular K + and Li + were as predicted. In interpreting the results of numerical modeling, one should keep in mind that the main limiting factor in solving physiological problems is the accuracy of experimental data, and not the fine details of the model.
The experimental data necessary for rigorous determination of the parameters for U937 cells under the balanced state include monovalent ion concentrations, cell water content, ouabain-sensitive and -resistant Rb + influx and the rate coefficients of equilibration of 22 Na, 36 Cl and Li + at low concentrations. The NKCC and KC transporters were not included in the model because of their low expression in U937 cells. However, our software allows their inclusion in the model if desired. For determination of the NKCC and KC rate coefficients, additional information would be needed. That information could be obtained, for instance, from measurements of the bumetanide-and DIOA-sensitive portions of the Rb + influx. Thus, the developed computational approach can be applied to any cell types.
The developed software enables one to estimate quantities that are difficult to measure directly, i.e. unidirectional fluxes, one-to-one coupled ion exchange, transmembrane Computed net and unidirectional Na + , K + , Li + , and Cl − fluxes via major pathways in cells equilibrated in RPMI+5 mM LiCl. The parameters used in the computation are as in Table 3; the flux units are defined in Table 1.
doi:10.1371/journal.pone.0153284.g006 electrochemical and electric potential differences, and the rate coefficients for ion transfer via specific pathways. Furthermore, it makes it possible to distinguish between the relatively simple effects caused by linking of fluxes due to the basic conditions of osmotic balance and electroneutrality and more complex effects caused by alteration of ion pathways per se, i.e. the rate coefficients of ion transfer. We believe that this software can be useful both in research (for example, for finding conditions where the effects of alteration of specific pathways on measurable characteristics are most pronounced) and in education, for demonstration of regulation of cell water, membrane potential and electrochemical gradients by individual channels and transporters.

Appendix Equations
The mandatory conditions of macroscopic electroneutrality and osmotic balance can be written when Li + is presented in the media or in cell as: The subscripts i and o indicate the intracellular and extracellular concentrations in mM, respectively; V is the volume occupied by intracellular water (ml), which assumed approximately equal to the cell volume; A is the total amount (in mmol) of membrane-impermeant osmolytes with mean valence z; [B] o is the external concentration of the membrane-impermeant solute (see Table 1 for symbols and definitions).
From Eqs 1 and 2 it follows that Dependences of ion fluxes via electroconductive channels and NKCC, KC, NC cotransporters and the pump on ion concentrations were formulated by the simplest way ignoring saturation phenomena [3][4][5][6]25]. The equation for Li + fluxes was similar to that for Na + , but without pump component. The Li/Na countertransport added which was formulated like the Na/H countertransport in [6].
To express the fluxes, the following equations used: dð½Li i VÞ dt ¼ Vfp L uð½Li i expðuÞ À ½Li o Þ=g À a½Li i þ J LC þ J LKCC À J LN g ð 9Þ Here u is the dimensionless membrane potential related to the absolute values of membrane potential U (mV) as U = uRT/F = 26.7u for 37°C and g = 1 − exp(u).
The left-hand sides of Eqs 6-9 represent the rates of change in the cell ion content. The right-hand sides express net fluxes via channels, the Na/K pump (β[Na] i , β[Na] i /γ), the Li/K pump (α[Li] i , α[Li] i /γ), cotransporters NC, KC, LC, NKCC, LKCC (J NC , J KC , J LC , J NKCC , J NKCC ), and countertransport Li/Na (J LN ). The rate coefficients p Na , p K , p L , p Cl characterizing channel ion transfer are similar to Goldman's coefficients. Fluxes J NC , J KC, J KC, J NKCC , J LKCC and J LN are related to intracellular and external ion concentrations as: Here i NC , i KC , i LC , i NKCC , i LKCC , and k P are the rate coefficients. By multiplying expression 1 by V and taking a time derivative, we obtain the following relationship: From Eqs 6-9 and considering Eq 11, a transcendental equation for potential u follows: FðuÞ ¼ expðuÞ À p Na ½Na o þ p K ½K o þ p L ½Li o þ p Cl ½Cl i þ b½Na i ð1 À 1=gÞ=u þ a½Li i ð1 À 1=gÞ=u p Na ½Na i þ p K ½K i þ p L ½Li i þ p Cl ½Cl o þ b½Na i ð1 À 1=gÞ=u þ a½Li i ð1 À 1=gÞ=u Algorithm Modeling of the flux balance in a cell requires a numerical solution of a Cauchy problem for the set of Eqs 6-9 [9,45]. The values of z and V compatible with the conditions of electroneutrality and osmotic balance have to be computed at the outset. Cell volume V and potential u are found at every step of numerical integration. The transcendental Eq 12 is solved in two stages. First, one finds the values of u that ensure that F(u) crosses zero at physically meaningful values of u. Next, the process of root refining ZEROIN from the calculation package FMM [46] is applied. Numerical integration of 6-9 can be accomplished by the explicit Euler method.
To monitor the approach to stationary state, time derivatives of sodium, potassium and chloride concentrations are displayed.

Executable file
This system of equations does not have a physically meaningful solution for every value of parameters or initial concentrations. This was taken into account in program development, and the following diagnostic tools were provided to assist the user: 1. In studying the dependence of the process on external ion concentrations, one enters new concentrations along with parameter kv = S oN / S oO , where the second subscripts stand for the old (O) and new (N) values.
If the parameter kv or the initial concentrations do not meet the requirements of electric or osmotic equilibrium, a message "BAD INITIAL DATA" is generated.
2. The program finds potential u in the voltage range +5 mV to -173 mV as sufficient to the possible physiological values and then stops. The negative limit may be changed if necessary. If u is not found because of unrealistic parameters or initial concentrations, the program reports "RANGE LIMIT".
3. If, in the process of computation, the internal sodium concentration drops below 0.1, a message "LOW SODIUM" appears, and the computation stops.
The executable file LIP, as well as the files: How to use LIP.doc, DATAP.txt, LIP.txt, and RESP control.txt is included in the zipped file LIP.zip available at www.cytspb.rssi.ru/lab_ vereninov/LIP.zip. The source code LIP-pas.zip can be obtained upon request from vereni-no@gmail.com, via1938@mail.ru, v.yurinskaya@mail.ru. It can be freely reproduced and distributed with explicit reference to or acknowledgment of the material author.