Protons in small spaces: Discrete simulations of vesicle acidification

The lumenal pH of an organelle is one of its defining characteristics and central to its biological function. Experiments have elucidated many of the key pH regulatory elements and how they vary from compartment-to-compartment, and continuum mathematical models have played an important role in understanding how these elements (proton pumps, counter-ion fluxes, membrane potential, buffering capacity, etc.) work together to achieve specific pH setpoints. While continuum models have proven successful in describing ion regulation at the cellular length scale, it is unknown if they are valid at the subcellular level where volumes are small, ion numbers may fluctuate wildly, and biochemical heterogeneity is large. Here, we create a discrete, stochastic (DS) model of vesicular acidification to answer this question. We used this simplified model to analyze pH measurements of isolated vesicles containing single proton pumps and compared these results to solutions from a continuum, ordinary differential equations (ODE)-based model. Both models predict similar parameter estimates for the mean proton pumping rate, membrane permeability, etc., but, as expected, the ODE model fails to report on the fluctuations in the system. The stochastic model predicts that pH fluctuations decrease during acidification, but noise analysis of single-vesicle data confirms our finding that the experimental noise is dominated by the fluorescent dye, and it reveals no insight into the true noise in the proton fluctuations. Finally, we again use the reduced DS model explore the acidification of large, lysosome-like vesicles to determine how stochastic elements, such as variations in proton-pump copy number and cycling between on and off states, impact the pH setpoint and fluctuations around this setpoint.

Introduction Acidification of intracellular organelles, such as lysosomes, is achieved by the V-ATPase proton pump. That said, how organelles set and maintain specific pH environments remains poorly understood. In addition to H + pumping, counter-ion movement is needed to oppose the buildup of a positive membrane potential [1]. The chloride-proton antiporters and chloride channels of the ClC family are thought to play this role in many organelles [2][3][4][5]; however, this role is controversial in the lysosomal where cation channels have been implicated [6]. Another determinant of lumenal pH is the proton leak across the membrane, which may be membrane mediated, occur through voltage-gated proton channels such as H V 1 [7,8], or through other unknown proteins. Lumenal pH is also impacted by internal buffer molecules and the numerous H + -dependent chemical reactions that occur in cellular compartments.
Many pH regulation studies have focused on assembling a parts list of the components involved in acidification, and in doing so, they have produced valuable information regarding the impact of specific proteins/molecules on pH. However, there are still gaps in these lists for most organelles with few studies producing a comprehensive dissection of the protein makeup as was done for the synaptic vesicle by the Jahn lab [9]. Additionally, most cell based studies have reported only ensemble averaged pH measurements for different organelles, potentially masking important compartment-to-compartment variations. Exceptions include a handful of studies on endosomes [10], synaptic vesicles [11], and reconstituted proteoliposomes [12]. In particular, the Krishnan lab nicely shows spatial and temporal variation in the distribution of pH values of maturing endosomes with final lysosomal values of 5.0 ± 0.16 for 80 individual compartments [10]. Nonetheless, the fluctuations in pH of individual compartments in time has not been presented in a cell based assay to our knowledge. In a similar vein, the Grinstein lab showed that lysosomes adjacent to the nucleus are more acidic than those that are distantly positioned [13], and such heterogeneity may have functional consequences related to nutrient response and organelle position [14]. Nonetheless, without knowing how the pH in isolated organelles vary we have an incomplete picture of organelle function. For instance, are all lysosomes near pH 5, or is there a wide distribution of values in the cell ranging from mildly acidic to very acidic? If this spread is too large, it is possible that some organelles may fail to achieve the chemical environment required for proper function. On a related note, how does the pH of an individual organelle change in time as enzymes and pumps activate and deactivate, compartments fuse together changing the chemical make up of the membrane proteins and lumenal contents, and the organelle ages?
Previous models based on ordinary differential equations (ODEs) have proven beneficial in teasing out the contribution of specific components to this rather complex processes. These models include components such as ATP-driven proton pumps, counter-ion fluxes, membrane potential effects, and buffering capacity [15], and they have highlighted the importance of the balance between V-ATPase proton pump and proton leak in setting pH along the secretory pathway [16], helped tease apart glutamatergic and GABAergic effects in synaptic vesicles [11], elucidated the role of ClC-7 antiporters in aiding lysosomal acidification [17], and generated hypotheses regarding the complex orchestration of plasma membrane, cytoplasmic, and ruffled border elements during osteoclast-mediated bone resorption [18]. Nonetheless, ODE models only provide an average description of a population of compartments, which may limit their usefulness for studying isolated compartments of a heterogeneous nature. In particular, the continuum approximation may fail for small compartments containing a handful of chemical species. A synaptic vesicle, for example, is typically 20 nm in radius and has to maintain a lumenal pH of approximately 5.5 to properly package neurotransmitters. The concentration of free protons at pH 5.5 in this small volume however corresponds to only * 0.06 free protons in number. The variables in ODE models can take on fractional values such as these, but noninteger values lack physical meaning.
Recent experiments in the Stamou Lab measured for the first time the acidification of single vesicles by single proton pumps [12]. AHA2, a P-type H + pump from the plant Arabidopsis thaliana, was purified and reconstituted into proteoliposomes, and acidification of isolated vesicles was monitored through the pH sensitive fluorophore pHrodo. Individual traces showed a great deal of variability both in time and between proteoliposomes, in part because they revealed for the first time that AHA2 stochastically switches between active and inactive states. An ODE model was constructed to fit the time-dependent pH traces, and although it provided reasonable parameter estimates and new molecular insight into AHA2-mediated proton leak, the model only describes the mean vesicular pH providing no insight into the experimental fluctuations. Moreover, some proteoliposomes were so small that they likely contained less than 1 free proton even at acidic pH values. To address these shortcomings and build a framework for answering the questions above, we developed a discrete, stochastic (DS) model of vesicular acidification that enforces the discrete nature of particles while also capturing information about fluctuations in the system. We revisited the AHA2 data, which is all taken from our previous publication [12], using a DS model that is equivalent in complexity to our previous ODE model of acidification and report here that both models provide very similar parameter estimates, even when free proton counts are much less than one. The DS model coupled with analytic results also shows that the noise in these experimental traces is very large and primarily reports on the fluorescent dye rather than the true proton fluctuations. We suggest conditions under which the fluctuations in the reporter dye would more closely match those of the free proton count. Finally, the DS model was used to predict pH fluctuations in lysosomal-sized compartments containing different numbers of proton pumps. While the DS model shows that such compartments can achieve different mean pH values by modulating the number of pumps recruited to their membrane, it also shows that the distribution of pH values around these means are relatively narrow and robust to stochastic activity and Poissonlike changes in pump number.

A discrete, stochastic model of acidification
We model the acidification of single vesicles using a system of chemical reactions following our previous ODE work [12]. These reactions account for proton and potassium leaks, a proton pump, lumenal buffer species, and impermeant, negative Donnan particles trapped in the lumen (Fig 1). These reactions determine the change in lumenal concentrations and the membrane potential, and the corresponding ODE for the number of lumenal H + , for example, is given by: H ½H þ �Þ |ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl } leak þ N A Vðk À B ½HB þ � À k þ B ½B�½H þ �Þ |ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl } buffer reactions ; ð1Þ where k P is the time-dependent proton pumping rate into the vesicle, V is the volume of the vesicle (assumed spherical throughout), k þ H and k À H are rate constants describing the passive movement of protons into and out of the vesicle, respectively, k À B and k þ B are rate constants for the dissociation and association of the buffer and proton, respectively, [HB + ] and [B] are the concentrations of the protonated and free buffer, respectively, [H + ] o and [H + ] are the concentrations of protons outside and inside the vesicle, respectively, and N A is Avagadro's number. Throughout, the subscript o refers to an extracellular value and lumenal values do not have a subscript. Many of these quantities are subject to further biophysical and chemical constraints.
Buffer reactions. While Eq 1 only shows a single species, our model has two explicit buffer molecules: solution buffer molecules distributed throughout the vesicle, B 1 , and lipidbound dye molecules called pHrodo, B 2 . In general, detailed balance relates the forward and reverse rates through the experimentally determined pK a value for each buffer (see Table 1) as follows: For all simulations, the value of k À B was 0:1 1 s for both buffer species, and the forward rates then followed from Eq 2 and are given in Table 1. The value of the reverse rate was selected such Protons in small spaces that the forward buffer rates (* 10 4 -10 5 L mol�s ) would be the fastest in the simulation. The change in protonated, N HB , and free buffer molecules, N B , for either species are then given by: where the constraint on the total number of buffer molecules, N HB + N B = N T , eliminates one of the equations. [B 1 ] is constant in all experiments, and N T = V[B 1 ]. pHrodo lipid-dye molecules are incorporated into vesicles at a mole fraction F of 1.5 to 1000, and then the vesicle geometry defines the average number of molecules and concentration: |ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl } where A L is the area per lipid headgroup, and A is the vesicle surface area. Membrane potential. We employ a physical model for the membrane potential [19] determined by summing the total charge in the vesicle and dividing by the membrane capacitance C: where N indicates numbers of specific charged molecules (free protons, protonated buffer, protonated dye molecules, and potassium, respectively), B is the concentration of negatively charged Donnan particles, and F is Faraday's constant (see Table 1 for values). Passive H + & K + permeability. Passive H + and K + membrane leaks were modeled as: where this equation for K + is equivalent to the leak term in Eq 1 with similar definitions, and the rate constants k þ X /k À X for X = K or H depend on Δψ. In the absence of a membrane potential, the flux vanishes when the internal and external concentrations are equal leading to In the presence of a potential, the base rates must be modified to obey the Nernst equation and provide the correct equilibrium condition. The simplest model for a monotonic cation that enforces reversibility is the following: where R is the gas constant and T is temperature. We chose to split the membrane potential dependence equally among the forward and reverse reaction rates: The membrane permeability (cm/s) is related to A (nm 2 ) and the reaction constant, k 0 X (L/s), via the equation: k 0 K for K + is zero prior to valinomycin addition and afterwards is set by Eq 10 and its membrane permeability (Table 1). k 0 H for H + is composed of two terms: The first term is direct membrane mediated leak that is always present, and the second term is proton pump dependent leaking that occurs once the pump stochastically turns off as described previously [12]. The model parameters in Table 1, P H and P AHA2 , refer to the membrane permeability calculated from k 0 lipid and k 0 AHA2 individually using Eq 10. Implementation of the stochastic model. The ODE-based model of vesicular acidification summarized by Eqs 1, 4, and 7 was implemented in the biochemical reaction network software, COPASI [20], using the time course feature. Only the internal concentrations were allowed to vary, and the model held all extracellular concentrations fixed. To simulate a discrete, stochastic (DS) process, COPASI employs a version of the classic Gillespie algorithm called the Next Reaction Method [21]. The reaction equations, rate constants, and starting concentrations from Fig 1 were entered into COPASI, and a trajectory of reactions and corresponding reaction times were stochastically generated. After each reaction was executed, the the membrane potential was updated according to Eq 6 as well as all affected rate constants. Each simulation was run for 750 s, and species counts were recorded every 0.02 seconds. For comparisons to deterministic ODE results, we simultaneously solved Eqs 1, 4, and 7 using the ode45 stiff solver algorithm in MATLAB (The MathWorks, Natick, MA).
Calculating pH. We calculated lumenal pH from the stochastic vesicle model results using three methods. First, post-simulation, the instantaneous free proton count at each time step was converted to a pH value based on the volume of the vesicle using the relation pH = −log 10 ([H + ]), which we simply refer to as pH. Second, we time-averaged the proton count over a time Δt, usually 0.5 s, and then took the À log 10 ð½H þ �Þ to determine the time averaged pH from the DS model. Third, we computed pH from the fluorescence of the pHrodo dye molecules to estimate pH in a manner identical to previously published experiments [12], pH hν . When pHrodo is protonated, it fluoresces, and a CCD camera detects the number of photons emitted in time windows from isolated vesicles. The pH was then determined from the photon count using a predetermined vesicle-size specific calibration curve.
We first solved the stochastic model to determine the number of protonated dye molecules in time. Next, we estimated the average photons emitted from N HB dye molecules using the equation: where c and b are trace-dependent constants described later. Realizing that emission is stochastic, we simulated the actual number emitted photons assuming they obey Poisson statistics: where P is the Poisson distribution with mean � N hn . We determined c and b for each experimental trace by using the same linear model to convert protonated dye molecules to photon count. To solve for the two unknown parameters, we used the average experimental pH (pH exp ) and N hν values prior to ATP addition (T1) and after maximum acidification was reached (T2) to obtain two constraint equations: where N HB 2 was obtained from pH exp using the Henderson-Hasselbach equation: where N T is the total number of dye molecules. After solving for c and b, the stochastic model reproduced the experimental photon count time series very well using Eq 13. Finally, the model pH in this third method is determined as follows: where importantly N HB 2 is not determined here from the model, but rather from Eq 12 as N HB 2 ¼ N hn =c À b, which includes photon shot noise. Our modeling efforts not only match the photon time series, but it also closely matches the experimental pH calculation curves in Ref. [12]. Data fitting. We used our models to fit experimental pH measurements using the Nelder-Mead search algorithm [22] to find the optimum parameter values that minimized the root mean squared deviation (RMSD) between the simulated and the experimental data across the entire time series: RMSD ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi where N is the number of experimental data points and X mod/exp is the model/experimental free proton count or pH, as described next. When computing pH directly from the free proton count in the DS model, we compared calculated free proton values to the experimental free proton count, which was back calculated from the experimental pH values. However, when we computed pH from the photon count in the DS model, we compared pH hν directly to the experimental pH using Eq 17. Finally for the ODE model, computed pH was again compared directly to the experimental pH. We chose to use free proton count in the first case, because pH is undefined when there are no free protons. Because the RMSD value for any given parameter choice varied from run-to-run, due to stochastic variation, we carried out the optimization on the average RMSD from 10 runs: err ¼ 1=10 P 10 j¼1 RMSD j . Optimum fits were determined by varying k P , k 0 lipid , k 0 AHA2 , relating to the pumping rate, passive membrane permeability, and AHA2-dependent proton leak. All other parameters listed in Table 1 were fixed. The initial interior proton concentration was set to the experimental value and was also assumed to be in equilibrium with the exterior pH, because ATP had not yet been added. The ratios of protonated to deprotonated buffers were set based on the total buffer concentrations and known buffer pK a values. Additionally, the radius, r, of each vesicle was measured in previous experiments [12] wherein individual radii were determined by converting diffraction-limited intensity spots to physical proteoliposome size [23]. These vesicle sizes were then calibrated via cryo-EM [12,23].
Experimentally, AHA2 H + pumping (I p ) begins with the addition of ATP to solution and then stochastically halts at some later time. We accordingly divided each trace into three sections: before, during, and after pumping. I p was non-zero only in the middle section, and P AHA2 was only non-zero after pumping due to AHA2.

The model reveals discrete pH jumps not observed in experiment
We first used the stochastic model to explore the acidification of a small vesicle * 50 nm in radius. While still larger than synaptic vesicles, such a small space (5 × 10 −7 pL) poses an interesting question when interpreting proton concentrations since a pH of 5 corresponds to * 3 free protons. To our knowledge, our recent work, led by the Stamou lab, represents the only experimental study of the acidification of isolated vesicles by single transporters, which we call Single Transporter Activity Recordings (STARs). These artificial proteoliposomes reconstituted with 1 or a few AHA2 proton pumps, verified with single molecule quenching of tagged AHA2 [12], ranged in size from tens of nanometers to hundreds of nanometers, and they serve as an excellent standard to compare against our discrete, stochastic model. Fig 2 shows the time resolved trace from a 58 nm vesicle (black curve in all panels). Near 180 s, ATP is added to solution to initiate pumping of the AHA2 proton pump, at which point the vesicle acidifies from the bathing pH of 6.5 to 5.5. Then near 500 s this transporter spontaneously turns off and protons leak out of the vesicle returning the intravesicular pH to the exterior value. We first carried out curve fitting to obtain a set of model parameters that most closely matches the experiment following the optimization scheme outlined in the Methods. The orange curve in Fig 2A is a representative stochastic run of the DS model using these optimized parameters. In this simulation, we are using the instantaneous number of free protons in the vesicle at each point in time and then converting that value to pH using the standard definition of pH (orange traces). It is immediately clear that the model produces a result that is dramatically different from the experimental trace. As expected, the model shows discrete changes in the state of the system when 0, 1, 2, 3, or 4 free protons are present in the lumen-pH 5.7 corresponds to 1 proton, 5.4 to 2, and so on. For a large portion of the beginning and end of the simulation trace, the instantaneous proton count is zero, resulting in an infinite pH value not visualized in Fig 2A. Despite the clear differences, there are similarities between the simulation and experiment. When AHA2 starts pumping, there is a general acidification that occurs in both traces, and when pumping stops both alkalinize. While the vesicle returns to the extravesicular pH of 6.5, the DS model returns to having no protons on average (pH ! 1). While we suspect that the DS model is an excellent representation of the true biological system, and that it likely produces the correct instantaneous pH values, why does it not match experiment?
In a first attempt to address this question, we accounted for the fact that the experiment does not report on the instantaneous pH. The fluorescence emitted from acidifying vesicles was collected over a 0.5 second window and reported every 2 seconds [12]. Thus, we time averaged the free proton count from the DS model over a 0.5 second window prior to converting to pH (pH). As can be seen in Fig 2B, this approach produces a continuous pH trace that no longer shows discrete steps in proton count and more closely matches experiment. The alkaline pH values correspond to fractional numbers of free protons, and they result from many snapshots of the vesicle with no protons and a few with one. Although this method of time averaging produced a simulated pH curve with the same mean pH values as the experimental trace, there are pronounced differences in the magnitude of pH fluctuations between the two curves. In the model, the fluctuations are quite large at alkaline pH levels and are suppressed as the vesicle's mean lumenal pH decreases. This trend is simply a consequence of the log scale of pH measurement (i.e. while variance in the number of free protons grows at acidic pH values, the log 10 suppresses the fluctuations); however, it was puzzling that the experimental traces showed no obvious pH dependent change in the magnitude of the fluctuations, as can be seen in the example curve in Fig 2. We also modeled the pH in the vesicle in a manner that exactly mimics the experimental setup to gain deeper insight into what the experiments are telling us. As is the case with many Protons in small spaces environmental reports, the fluorescence is not a direct readout of the free molecule/ion of interest, but rather it is the result of binding of that molecule/ion to a reporter, which then changes its fluorescent properties. Here, the reporter is the pH sensitive dye, pHrodo, which resides in the membrane leaflet due to conjugation to a DOPE lipid. For a 58 nm radius vesicle, there are approximately 88 pHrodo lipids in the inner leaflet. Given the pK a of 5.72 of this lipidated form of pHrodo (see Ref. [12]), on average, initially 12 lipids are protonated, and this number increases to 51 as the system acidifies. As described in the Methods, we simulated the stochastic emission of photons from each bound pHrodo, which also increases in time (Fig 2C  inset), and we plotted the predicted lumenal pH from this signal in the main panel, which we call pH hν . pH hν not only tracks the mean experimental pH, but also qualitatively the experimental fluctuations. The two approaches for interpreting the DS model exhibit very different fluctuations because one is based on proton count (panel B) and one is essentially based on bound dye count (panel C), and both quantities and their variances scale very differently with pH. Proton count changes exponentially with pH, while the bound dye count, and hence the photon count, increases almost linearly over the experimental range of interest (Fig 2C inset). Although the pH hν best fits the experimental results, it highlights the fact that the experimental measurements are not directly capturing information about the free protons, but rather the fraction of protonated pHrodo molecules. Hence, we assert that the instantaneous pH calculations determined directly from the DS model provide the best indication of the free proton dynamics in these compartments.

Parameter estimates are model independent
Next, we wanted to apply our two methods of fitting the stochastic model to a range of STAR data recorded from vesicles of different sizes, quantitatively compare the parameter estimates from each approach, and importantly, compare the results to solutions derived from ODE calculations. The vesicles in Fig 3 range from 75 to 197 nm in radius (a volume difference of approximately 18×), and for each trace, we fit the DS and ODE models to the data and identified the optimal AHA2 H + pump rate (I P ), membrane proton permeability (P H ), and the AHA2-dependent proton permeability (P AHA2 ) that occurs after the pump stops (Table 2). Since the DS model is stochastic, any individual simulation may provide a poor match to a given STAR even if the parameter values are optimal. Therefore, for each point in parameter space, we ran the stochastic models 10 times and averaged the goodness of fit from each run to provide a representative score, which was then optimized using a Nelder-Mead search algorithm [24].
Model fits (orange curves) are reported in separate columns superposed on the corresponding experimental trace (black). The ODE model (Fig 3 left column) produces piece-wise smooth curves that overlay well with the experimental data, throughout the trace; however, as discussed earlier, the ODE model does not reproduce the fluctuations in the trace. The greatest deviations between the fit and the data occur in traces 2 and 3 around 150-175 seconds. At these points in time, the experimental trace changes rapidly, either due to a stochastic fluctuation or to the pump rate changing in a manner not encoded in the model. For the DS model (middle column), pH produces discrete pH jumps, as discussed earlier, not observed in the data recorded from the small 75 nm radius vesicle (trace 1), but the discrete nature is obscured in the larger vesicles (traces 2-4) because they contain many more protons. Incidentally, the ODE model predicts a free proton concentration less than 1 for most of trace 1. As the vesicle acidifies in the DS method for traces 2-4, the noise decreases for reasons discussed previously, and this is not observed in the data. Lastly, interpreting pH changes in the DS model by explicitly modeling the emitted photons (right column) produces the closest fit to the experimental  Table 2. Each row is a different experimental trace plotted in black, and all simulations are orange. The first, second, and third columns contain the ODE, DS model, and DS model using photon count results, respectively. In the last column, the calculated pH is pH hν . All parameter values can be found in Table 2.
https://doi.org/10.1371/journal.pcbi.1007539.g003 data for the small vesicle (trace 1), and as discussed before, it exhibits similar noise at both high and low pH. However, the magnitude of the fluctuations for the larger vesicles is smaller than what is observed experimentally, which may arise from other environmental sources not included in our model. As we show in Fig A in S1 File, the membrane potential is predicted to be negligible for the DS model due to sufficient clamping by potassium via valinomycin. Importantly, despite qualitative differences in the ODE and DS model, they both predict very similar parameter estimates ( Table 2). Both DS methods predicted pump values (I P ) within 10% of the ODE value, while the passive membrane permeability values (P H ) are all between 0 and 30%. The ODE results for P H match the DS model interpreted from photon counts better, but we are unsure if this is true or the result of analyzing a small number of traces. The greatest variability for the DS method predictions is in the estimated P AHA2 values, which differ for the two DS methods by a factor of 2-4 for traces 2 and 3. We note that estimates of this quantity, which represents the leak of protons through the transporter once it stops pumping, are prone to error. Since the protons leave the vesicle so quickly over a short period of time, as can be seen near 500 s in Fig 2, the value is poorly constrained from above. That is, increasing P AHA2 above a critical value produces nearly identical looking fits with similar RMSD values. Thus, we do not think that these differences are crucial. Rather, the closeness in all predicted values across all models and methods suggests that it is not essential to model the fluorescent response of the dye molecules to capture the mean behavior of the system nor does one need to model the discrete nature of the proton movement, as we elaborate upon next.
One reason why the mean parameter estimates for the ODE and DS models are close is because the ODE model can be interpreted as an average of multiple DS model simulations, as we show in Fig 4. Here, we used the parameter estimates for trace 2 and ran the DS  Table 2, r = 145 nm, is displayed. An average of 10 and 100 DS simulations with the DS method predicted parameter values is displayed as well. As the number of averaged DS simulation increases, the ODE model fit is approached.
https://doi.org/10.1371/journal.pcbi.1007539.g004 model one time (blue), averaged over 10 runs (red), or averaged over 100 runs (green) and plot the resulting pH. As more traces are averaged, the fluctuations are suppressed, and the final curve is indistinguishable from the smooth trace produced from the ODE model (black dash). Evidently, the ODE model accurately describes the mean lumenal pH of the liposome over time. The similarity between ODE and DS model predictions suggests that the deterministic, continuum model is detailed enough to provide parameter estimates, even for very small vesicles with large pH fluctuations, and proton counts that are at times less than one.
Finally, we carried out sensitivity analysis on the ODE model for trace 2 to determine which parameters most greatly impact the dynamics and steady state behavior. Given the similarity between mean parameter estimates for both DS and ODE models, we expect that these results apply equally to both. As with several of our previous pH regulation studies, the pumping rate, proton permeability, and surface area all impact the steady state pH the greatest, while the properties of the primary buffer, AHA2 pumping rate, and volume impact the rate of acidification (Table A in S1 File).

Vesicle radius and mean pH are the greatest determinants of pH fluctuations
A central question in intracellular pH regulation is how the pH setpoint is established in an organelle. A related question is, once set, are vesicles prone to large pH fluctuations that would change the proton motive force (PMF) across the membrane influencing H + -dependent transporters, ion channels, and lumenal pH-dependent enzymes? Since the DS model captures fluctuations, we wanted to address these points computationally by identifying the parameters that most greatly influence the fluctuations in the free proton count, and correspondingly the pH. Our analysis revealed that vesicle radius and mean lumenal pH exert the greatest impact on the number of total free protons in the system, and hence the H + fluctuations. In Fig 5A, we carried out a series of calculations in which we varied the radius of an idealized liposome from 70 to 270 nm with the external pH constant at 5.5. There were no proton pumps, only a passive H + permeability, and the mean lumenal pH mirrors the external value for all sizes. We see as the radius of the vesicle increases the standard deviation of free proton count, σ(N H ), also increases ( Fig 5A). Next, we kept the vesicle radius constant at 150 nm, but varied the bath pH from 4.5 to 6.3, and as the lumen acidifies σ(N H ) increases (Fig 5B). If the free proton count obeys Poisson statistics, we would expect the standard deviation to scale like the square root of the number of protons. For panel A, as the vesicle radius increases at fixed pH, the number of free protons is proportional to the vesicle volume, and we would expect s N H ð Þ � ffi ffi ffi ffi r 3 p ¼ r 3 2 , as it does (dashed curve). Meanwhile, the vesicle size is fixed as the pH changes so the number of free protons scales like 10 −pH , and the noise scales like ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 10 À pH p as expected (dashed curve in Fig 5B). In general, the theory is an excellent match to the numerical results (points), supporting the claim that the H + fluctuations obey Poisson statistics at steady state in our model. In panel Fig 5C, this correspondence is explicitly shown for a single dynamic trace of H + from the DS model at pH 5.5 and r = 150 nm (diamond in Fig 5A and 5B). Binning the counts on the right reveals a histogram with values from 0 to 8 that is again well fit by a Poisson distribution (solid line).
Experimentally, it is more useful to examine the fluctuations in pH units than proton count, which we see scales very differently with radius and mean pH in panels Fig 5D and 5E, respectively. Standard deviation in pH decreases as vesicle radius increases, and it increases as mean pH increases. These trends are captured by the dashed curves in panels D and E, which follow from the transformation to pH: where we note that N H scales like r 3 or 10 −pH , and we recall how σ(N H ) varies with r and pH from the last paragraph. Thus, fluctuations in pH are suppressed in larger, more acidic vesicles-the opposite of what is observed for proton count. Finally, we wanted to determine if active pumping changed the Poisson-like characteristic of the noise, so we carried out fluctuation analysis at steady state for an acidic vesicle with active pumping. For different sized vesicles the proton pump rate was varied to achieve different steady state pH values, and the resulting standard deviation of pH was recorded (panel F). σ(pH) matches the values in panels D and E (diamonds), suggesting that the process of proton pumping-in our simple modeldoes not influence pH fluctuations. Furthermore, we explicitly show for one of the DS simulations with pumping that at the acidified steady state pH, the distribution of free proton counts is still well fit by a Poisson distribution (panel F inset).

Fluorescent reporters cannot reveal the true pH fluctuations in a compartment
Important information regarding the nature of complex physical systems can often be extracted from the natural fluctuations. However, the proton noise analysis presented in  Table 1. is based on a very simple model, and we directly analyzed the noise in the true proton count, which cannot be done experimentally. If real pH regulatory systems harbor more interesting feed back systems that deviate from Poissonian statistics, we would have to extract this from a fluorescent reporter, as is often done in single molecule experiments. Thus, we wanted to determine if the true fluctuations in the free proton count could be related to the fluctuations in the fluorescent reporter. The pH-dependent fluctuations in pH reported by the dye are complex, and we carried out analytic analysis to better understand it. The noise in photon counts emitted by the dye arises from fluctuations in the protonated number of dye molecules, which follows from binomial statistics as well as photon shot noise, which we assumed to be Poisson distributed in our DS model. As derived in the supporting information in S1 File, applying the transformation in Eq 18, we arrive at the expression for the standard deviation in pH hν expected from the model: where f ¼ ½B 2 �=½B T � ¼ 1=ð1 þ 10 pHÀ pK 2 a Þ is the fraction of protonated dye molecules, λ is the mean photon production count per protonated dye molecule per collection time, and N dye is the total number of buffer molecules. While this mathematical expression fits our simulated data very well, it also makes it clear that there is no relationship between the fluctuations in the free proton count and the experimental fluctuations in the measured pH-the former obeys Poisson statistics and the later is dominated by binomial statistics of the bound dye (Fig B in  S1 File). Additionally, the covariance between the number of free protons and the number of bound dye molecules is very weak.
Next, we calculated the noise in the STARs to attempt to validate Eq 19. To do this, we carried out a denoising procedure to extract the fluctuations around the pre-valinomycin bath pH, and around steady-state acidified regions from over 100 single vesicle, single transporter traces. First, we identified by visual inspection a set of vesicle traces which remain acidified for an extended period. Next, we used tvdip.m [25], a total variation denoising [26] routine, to aid in selecting a steady-state time window from the most acidified region of each trace. We rejected traces, or narrowed time windows, which contain pH spikes consistent with transient pump off states. Finally we selected the subset of traces which acidify to 5.71 < pH < 5.81, as this 0.1 pH window contained more traces (15) than any other 0.1 pH window. We then calculated the Fano factor [27] (the variance divided by the mean) of the raw photon count for these 15 traces both in the steady-state acidified region, and also at pH 6.5 using the first 84 data points from before valinomycin addition.
For a Poisson processes, the Fano factor is 1 as the variance equals the mean; however, the fluctuations in raw photon count are super-Poissonian with mean Fano factor values 15-35 (Fig C in S1 File). Additionally, this noise is nearly constant between the two pH values. We calculated an expected Fano factor of 20 (see S1 File), in good agreement with the experimental distributions and confirming our prediction that the experimental noise is dominated by the dye and not related to the true pH.

Time-dependent pH fluctuations in an acidic compartment
For cellular compartments, there are additional considerations that we have not yet included that could cause deviations from the desired pH setpoint, namely the total protein copy number and the realization that pumps-like channels-cycle between on and off states [12]. If during organelle biogenesis a compartment acquires too many or too few proton pumps (or some other critical protein), the pump-to-leak ratio would be skewed pushing the lumen more or less acidic. Likewise, if pumps turn on and off stochastically, like AHA2, then the compartment may experience large fluctuations in time. For eukaryotic cells, the V-ATPase is responsible for acidification [28], and while it is not known if these pumps normally cycle on and off, it is known that they respond to cytoplasmic queues such as low glucose [29], and it is reasonable to assume that they do stochastically inactivate.
Here, we simulated a 340 nm radius vesicle, which is the typical size of a lysosome [17]. Starting at a pH 6, we carried out 6 DS model simulations each with a different number of total V-ATPases ranging from 20 to 640 total pumps (Fig 6). All pumps have the same proton pump rate of 100 H + /s [30], but they can each independently activate and inactivate stochastically. For AHA2, the mean dwell time in the on state is 273 seconds, and the off time is comparable under high ATP concentrations [12]. Based on this information, we modeled V-ATPases cycling between active and inactive states by pulling from an exponential distribution with a mean lifetime of 273 seconds for each pump. We were thus able to simulate a population of proton pumps randomly turning on and off, and explore how the steady state pH was influenced.
As expected, vesicles with fewer active pumps did not acidify as much as vesicles with far more active pumps (Fig 6A). Initially we thought that as the total pump number (N) increased the vesicle would exhibit increasingly smaller noise in the steady state pH because the fluctuations in the active pump number scale like ffi ffi ffi ffi N p , while the mean number of active pumps scales like N. To quantify these fluctuations, we binned the steady state pH values produced by the DS model in panel B. While it is true that the vesicle with 640 pumps exhibits the smallest pH fluctuations, the width is relatively uniform across the entire range of pH setpoints with a total spread of about 0.5 units at neutral values and 0.25 units at acidic values (Fig 6B), which is consistent, but larger, than the more idealized results from Fig 5F. Thus, the other sources of noise in the system, particularly fluctuations in the number of active pumps, contribute to a greater spread in pH. Lastly, we wanted to explore one mechanism of pH regulation in which there is no active sorting of proton pumps during lysosomal maturation, and N is simply equal to the compartment surface area times the pump density. In this case, the expected number of pumps will be Poisson distributed with a mean N. In order to estimate how an entire ensemble of vesicles would behave, we carried out 5 simulations of a 340 nm radius vesicle with 599, 623, 640, 657, and 682 proton pumps corresponding to 5 th , 25 th , 50 th , 75 th , and 95 th percentile of a Poisson distribution with mean number of N = 640 pumps, respectively. We then plotted the frequency of observed pH values, as in panel B, from each simulation scaled by the probability of observing that number (Fig 6C). As expected, the spread in values is larger than the spread observed from N = 640 alone, but the width is still very small on the order of 0.1 pH units.

Discussion
The process of vesicular acidification is one that is biologically essential, but yet still not completely understood. Mathematical models have been useful in interpreting experiments to help understand how different organelles achieve different pH values. The success of past ODE models in particular, in addition to the limitations of these models, motivated us to develop the discrete, stochastic (DS) model of vesicular acidification presented here. The DS model enforces the discrete nature of particles within cellular compartments and describes random fluctuations of pH in cellular compartments-two features that are absent in ODE models. We utilized the DS model to estimate the AHA2 proton pump rate and liposome H + membrane permeability from experimental pH recordings of isolated vesicles (STARs) from the Stamou Lab. One of the most surprising and important results of our present study is that DS and ODE models provide quantitatively similar parameter estimates ( Table 2), suggesting that ODE models are sufficient for describing intracellular ion regulation as long as only mean values are needed and not information regarding the fluctuations.
Measurements of protein and RNA variability have been used to identify unregulated housekeeping genes from Poisson distributed copy numbers [31], to understand how translational efficiency controls noise [32], to reveal how noise can give rise to population variability [33], and to identify how HIV employs a positive feedback loop to control latency [34]. Here, we made an analogy between protein/RNA copy number and proton counts in an attempt to understand the noise in acidic compartments. Unlike many single molecule experiments, the reporter is not the actual tagged molecule of interest, but rather it is related to the number of proton-bound dye molecules. Another difference is that the molecules of interest in gene regulatory network studies are unbounded since the cell can essentially make an unlimited number of molecules, while here the number of reporter molecules is fixed, which means that the expected noise from the reporter is binomial, not Poissonian, and ion fluctuations are not necessarily reflected in fluctuations in the number of bound dye molecules. Therefore, it follows that the fluctuations in pH determined from an experimental trace cannot be related back to the true pH fluctuations in the vesicle, because the experimental noise is dominated by fluctuations in the bound dye molecules as predicted by Eq 19, which is consistent with experiment ( Fig B in S1 File). We also expect that instrumental noise, i.e. from the camera and other sources, will have a small but non-negligible contribution to the fluctuations in pH measurements as well, and if this noise is Poissonian with respect to the incident photon count, it will only be about 3-10% of the total noise. One direct consequence of this realization concerning the fluctuations in the dye is that we cannot tell the true range of pH values in individual vesicles if these excursions are rather short lived, which may impact our understanding of the range of possible cellular reactions that can take place in a given compartment if its mean pH is shifted away from a value needed for catalysis.
Thus, simulation may be one of the best ways to assess the true range of pH fluctuations in isolated compartments. Nonetheless, fluorescent reporters do predict the mean pH in the compartment, and when averaged over times that are sufficiently long compared to the fastest proton dynamics (e.g. hopping on and off of buffer molecules), experimental changes in mean pH are valid and interpretable in terms of our ODE and DS models. In our simple system, the variance in the pH decreases with vesicle size and mean pH (Fig 5D and 5E), but this trend is less pronounced when multiple proton pumps are present and randomly activated and inactivated. The stochastic switching between on and off states produces single vesicle pH distributions that have similar variance across a pH range from 4.5-5.6 ( Fig 6B). We even considered a population of lysosome-sized organelles of fixed radius whose pump numbers were random, but proportional to the surface area, and the pH spread in this ensemble was only about twice the spread observed for a single vesicle (640 pump data in Fig 6B compared to Fig 6C). Thus, based on our initial work here, the relatively large size of the lysosome aids in suppressing noise from stochastic elements in the environment, so that it can achieve an acidic pH with a tight distribution. Now that we have a better understanding of pH fluctuations in a simplified system, we intend to build more realistic DS models based on our previous ODE studies of endosomes [15] and lysosomes [17] with a particular interest in exploring how fluctuations in pH are influenced by specific channels, transporters, and other molecules localized to each compartment. We will also construct models of synaptic vesicles, which we expect to exhibit large variations in pH given their small size. As stated earlier, there is little-to-no experimental evidence for fluctuations or steady state distributions of pH in different organelles, and the biological significance of these variations, if they do exist, is scarce. That said, subtle changes in mean pH can have a drastic impact on physiology, as exemplified by the recent finding that ClC-7 gainof-function mutations that results in lysosomal hyperacidifciation by only 0.2 units are associated with delayed myelination, hypopigmentation, and several other adverse effects [35].
Based on examples like this, we expect that compartment-to-compartment variability coupled with ionic fluctuations in time will emerge as important determinants in the homeostatic regulation of key cellular processes, and it will only be a matter of time for the biochemical tools and experiments to mature to the point where these phenomena can be measured and brought to light. (CSV) S1 Model. COPASI model. The discrete, stochastic model created with COPASI. The file format is XML, and it can be opened and run with COPASI [20]. Parameters such as liposome radius, AHA2 pumping rate, pump on and off time, and membrane H + permeability must be set before running the model (Model > Biochemical > Global Quantities). The model is run by selecting the Time Course feature under Tasks. (CPS) S2 Model. MATLAB ODE model. ODE model of vesicle adicification due to single molecule proton pumping. Description of the inputs can be found in the code. As an example, to reproduce ODE trace #2 from Fig 3, set PH = 3.9; IP = 488; PAHA2 = 65; ph0 = 6.5; radius = 145; pka2 = 6.15; on = 180; off = 506; tfinal = 750; in Matlab (from Table 2), and then run: [t, ph] = S2_Model([PH, IP, PAHA2], on, off, tfinal, ph0, radius, pka2); plot(t, ph). (M)