ORIGINAL RESEARCH article
Stochastic Model of Acidification, Activation of Hemagglutinin and Escape of Influenza Viruses from an Endosome
- 1Applied Mathematics and Computational Biology, Ecole Normale Supérieure, Paris, France
- 2Department of Biology, Molecular Biophysics, IRI Life Sciences, Humboldt-Universität zu Berlin, Berlin, Germany
- 3Institute of Chemistry and Biochemistry, Free University Berlin, Berlin, Germany
- 4Newton Institute and Department of Applied Mathematics and Thearetical Physics (DAMTP), Cambridge, United Kingdom
Influenza viruses enter the cell inside an endosome. During the endosomal journey, acidification triggers a conformational change of the virus spike protein hemagglutinin (HA) that results in escape of the viral genome from the endosome into the cytoplasm. It is still unclear how the interplay between acidification and HA conformation changes affects the kinetics of the viral endosomal escape. We develop here a stochastic model to estimate the change of conformation of HAs inside the endosome nanodomain. Using a Markov process, we model the arrival of protons to HA binding sites and compute the kinetics of their accumulation. We compute the Mean First Passage Time (MFPT) of the number of HA bound sites to a threshold, which is used to estimate the HA activation rate for a given pH (i.e. proton concentration). The present analysis reveals that HA proton binding sites possess a high chemical barrier, ensuring a stability of the spike protein at sub-acidic pH. We predict that activating more than 3 adjacent HAs is necessary to trigger endosomal fusion and this configuration prevents premature release of viruses from early endosomes.
For most viruses, the initial step of infection starts when the viral particles bind to specific receptors and enter the cell through the membrane, inside an endosomal compartment (Figure 1). Viral particles are then transported inside the endosome, from the cell periphery towards the nucleus. Several modeling approaches, including kinetics rate equations , stochastic modeling [2, 3] and mechanics of molecular binding  have been developed to describe how membrane receptors are activated and engaged into endosomal pathways. However, little attention has been devoted to study viral trafficking inside an endosome, which is a critical and limiting step in replication and more generally to gene delivery [5–8].
Figure 1. Structure and endosomal trafficking of the Influenza virus. (A) Influenza is an enveloped virus. Main spike proteins anchored in the envelope are the neuraminidase (NA) and the Hemagglutinin (HA). Protons can access the core of the virus through M2 channels. Main matrix protein is M1 protein. Viral genome of the virus is composed by eight viral ribonucleoproteins (vRNPs). (B) Influenza virus enters the cell via receptor-mediated endocytosis and progress rapidly toward an early endosome. Then, maturation is associated with an acidification of the endosome lumen and a retrograde transport of the endosome along the microtubules toward the nucleus, the destination of vRNPS for virus replication. The final destination of many endosomes is the degradative lysosomes. Thus, the timing of Influenza virus escape has to be tightly regulated to avoid degradation in lysosomes while delivering the genetic material close to the nucleus.
Cytoskeleton retrograde flow plays a key role for the Influenza virus transport inside the endosome toward the cell nucleus and to ensure a safe delivery of its genome near the nucleus, before replication [9, 10]. During this transport, the endosome can fuse with lysosomes, leading in that case to viral degradation. Thus, escaping the endosome at the right time must be tightly regulated to ensure that genes are released as close as possible from the nucleus, while avoiding degradation. Our goal here is to develop a first-principles model and the associated stochastic analysis to study this optimal escape time, and how it is controlled by acidification and conformational changes of viral proteins. The stochastic model developed in Lagache et al.  is inappropriate to compute the viral escape time based on the activation of a single molecule, and we develop here a different approach based on the mass-action equation for acidification.
The genome of Influenza virus is encoded by viral ribonucleoproteins (vRNPs) enveloped in a membrane. These vRNPs must translocate into the nucleus  (Figure 1) for reproduction. Endosomal escape is ensured by fusion between the endosomal and Influenza virus membrane. This fusion is mediated by a low pH conformational change of the viral glycoprotein hemagglutinin (HA) (Figure 1A). We account here for the detailed properties of the glycoprotein HA, composed of two linked subunits HA1 and HA2, the latter anchoring HA to the viral envelope. At neutral pH, HA is not active (in a non-fusogenic state), but as the pH decreases due to acidification (proton entry into the endosome), a partial dissociation of the HA1 subunit results in a spring-loaded conformational change of HA2 into an active (fusogenic) state . Consequently, the residence time of the Influenza virus genome within an endosome before fusion depends on the kinetics of endosome acidification. Yet, the absence of direct in vivo measurements of these parameters makes the endosomal step of virus infection difficult to analyze both theoretically and experimentally. To estimate the timing of the pH-driven fusion of Influenza viruses, the model we develop here accounts for the main kinetic parameters of the fusion process: endosomal acidification, binding of protons to HAs and independent activation of multiple HA neighbors, leading to membrane fusion and release of the genome into the cytoplasm.
The manuscript is organized as follow: Section 2.1 presents the kinetic model for endosomal acidification, calibrated to experimental data (Figures in the SI). The model depends on the buffering capacity of the endosome, membrane leakage and proton pumping rate that controls proton fluxes inside the endosome. In Section 2.2, we model the discrete and cumulative binding of protons to HAs using a Markov jump process . We find an analytical expression for the kinetics of HA conformational change at a fixed proton concentration, by analyzing the mean first passage time (MFPT) equation for the number of bound protons to a given threshold. In our previous work , we developed a jump model for the conformational change of active proteins for the escape of non-enveloped Adeno-Associated Viruses (AAV) from a vesicle. Contrary to the assumptions of , the binding rates of protons to HAs are non-linear  and thus we obtain here a different analytical expressions for the MFPT of bound protons to the critical threshold. Finally, in Section 2.3, we combine the kinetic models of acidification and HA conformational change and estimate the rate of HAs' activation inside the endosome. While the conformational change of a single protein for AAV is sufficient to lyse the endosome and release genes in cytoplasm, the Influenza virus is covered by 400 HAs and several adjacent HAs seems required for fusion.
We predict here the mean number of fusogenic HAs in the endosome and use Monte-Carlo simulations to compute the time needed for neighboring HAs to change conformation in the contact zone between the viral and endosomal membrane: we find that at least three adjacent activated HAs are necessary to trigger membrane fusion [16, 17], a cooperativity process that should prevent premature fusion. Some predictions are tested experimentally using co-labeling viruses and endosomal markers, confirming that intracellular fusion of viruses mainly occur in maturing endosomes (ME).
2.1. Kinetic Model of Endosomal Acidification
The present model of endosomal acidification is based on computing the free number of protons Pe(t) at time t in the endosomal compartment. The protons enter with an entry rate λ(t)S through the V-ATPase proton pumps (S is the endosomal surface area and the rate λ(t) is associated with the proton pumps activity) and can escape with a leakage rate Lext(t), but can also bind to endosomal buffers.
The proton pump rate λ(t) is mainly determined by the membrane potential Ψ(t) (Figure 11 in Grabe et al. ), which depends on the endosomal concentrations of several cations (H+, K+, Na+ …) and (Cl− …). The ionic concentrations inside endosome are tightly regulated by channels, exchangers and leak and in particular, by raising the interior-positive membrane potential, Na-K ATPase exchangers have been proposed to limit the acidification of early compared to late endosomes .
2.1.1. Mass Action Law for Free Protons
To derive the time-dependent equations for the free protons, we use the balance of fluxes: the fast equilibrium between fluxes determines the number of protons ΔPe(t) entering the endosome during the time step Δt
Entering protons are rapidly bound to endosomal buffers that we model using an ensemble of acid-base reactions :
where ki (resp. ) are binding (resp. unbinding) rate constants of protons to weak bases Bi, for 1 ≤ i ≤ n. The binding rates of entering protons to the endosomal basis can be reduced to a single constant that we call the effective buffering capacity of the endosome. The kinetics equation for the number of free protons Pe(t) inside an endosome is (Section 4)
When the proton leakage is counterbalanced by the pump activity, after a time long enough, the pH reaches an asymptotic value pH∞, where the endosome cannot be further acidified. This value is given by
Consequently, the rate λ depends on pH∞ with
and Equation (3) can be rewritten as
To conclude, we derived here a first order kinetic model (Equation 32) for the endosome acidification, based on the rapid equilibration of protons with buffer (see Equation 36). However, Equation (6) alone is not sufficient to account for endosomal maturation, because the final pH∞  and the permeability L decreases with the endosomal maturation  as they depend on time, as we analyse below.
2.1.2. Modeling pH Change and Acidification of an Endosome
Acidification in live cell imaging and the transition from an early endosome (EE) to a late endosome (LE) is monitored by a gradual exchange of Rab5/Rab7 proteins . We approximate here the kinetics of the ratio Rab5/Rab7 (Figure 4C in Rink et al. ) by a sigmoidal function
with the two free parameters: the half-maturation t1/2 and Rab conversion τc times. The steady-state pH∞(t) relative to the amount of Rab7 is given by
Thus we propose that the permeability rate follows the Equation
2.1.3. Acidification Model Calibrated from Live Cell Imaging Kinetics
We now explain the calibration of the acidification model to experimental data: First, we fitted Equation (7) to the experimental data (Figure 4C of Rink et al. ) where the lag time between initiation and termination of the Rab5/Rab7 permutation is estimated to 10 min., leading to a time constant for τc = 100s.
We use data from endosomal acidification in MDCK cells where the pH inside endosomes decreases very quickly within the first 10–15 min (Figure 2) to reach a steady-state pH around 5.5 after 20 min, in agreement with . The steady-state pH is and for early and late endosomes respectively . Thus, we calibrated the permeability constant L and Rab conversion kinetics by solving numerically Equation 6 and fitting the experimental acidification curve (Figure 2). We found that the permeabilities of early and late endosomes are and , respectively, and the half-maturation time is t1/2 = 10 min.
Figure 2. Endosomal acidification. The kinetics of acidification obtained with intracellular fluorescence microscopy (red line. Mean ± SEM) is compared with coarse-grained modeling (Equation 6, black line. Model parameters are summarized in Table 1).
2.1.4. Proton Influx Inside the Viral Core and Buffering
The last step of the kinetic model includes the buffering of protons to viral core components, defining the buffer capacity. Indeed, protein buffering capacity, influx of protons through M2-channels inside the viral core (Figure 1A) and the presence of viruses inside endosomes influences the overall buffering capacity of the endosome and acidification. To compute the influx through each viral M2 channel, we use a first order kinetics , summarized in the chemical equation
When a proton Pe binds a free M2 protein channel with rates ke(binding) and (unbinding), it is transported inside the virus core with a rate , while exit occurs with a rate kv. At steady state, the inward flux in a single virus is computed from Equation (10) (see )
where nM2 is the number of M2 channels per viral particle, Pv is the number of free protons inside the viral core and
To extract the buffer capacity of a virus, we accounted for the viral genome, the internal viral proteins and unspecific buffers that can be reached through the M2 channels . The most abundant internal proteins are M1 (3, 000 copies per virus) and the nucleoproteins (NP, 330 copies per virus)  (Figure 1A). Proton binding sites of viral proteins are the ionogenic groups in their amino acid side chains , and the main ionogenic buffers in the endosome pH range are the aspartic acid (Asp, pKa = 3.9), the glutamic acid (Glu, pKa = 4.32) and the histidine (His, pKa = 6.04) . Closely related binding sites can have strong influences on each other due to electrostatic interactions. In addition, the three-dimensional protein folding can hinder the accessibility of some residues to the solvent and protons.
Consequently, calculations based on the three-dimensional structure of the protein are necessary to determine the buffering capacity to pH. Using the spatial organization (crystal structure) of viral proteins, the overall buffering capacity βi of the viral core is given by
where is the buffering capacity of the lumen inside the virus, and , and are the buffering capacity of M1 and NP proteins, and viral RNA (see Section 4).
Similar to the flux Equation (6), the number of free protons Pv(t) contained in viral core at time t determines the influx of protons through M2 channels (Equation 11) and satisfies equation
By solving numerically Equation (14) with the initial conditions and , we estimate that = 60, 000 protons enter the viral core during endosomal maturation. Using Equation (6) for the endosomal acidification kinetic, we find that more than 20, 000, 000 protons bind to the endosomal buffers during acidification of an endosome with a radius re = 500 nm (Table 1). Thus, the buffering capacity of a single virus should not influence the endosomal acidification. However, the number of protons that bind to endosomal buffers drastically decreases to 175, 000 buffered protons when the endosomal radius is reduced re = 100 nm. In addition viral particles may accumulate during the endosomal journey . Thus, for multiplicity of infection (MOI) and viral accumulation in endosomes, the viral buffering capacity may significantly affect the acidification kinetics of small and intermediate size endosomes.
2.2. Markov Jump Model of HA Conformational Change
Although the number of protons entering in the endosome is quite huge, as discussed in the previous section, the actual number of free protons defining endosomal pH is surprisingly low (~300 at pH 6 in an endosome with a radius of re = 500 nm). In addition, there are few proton binding sites on a single HA that trigger a conformational change , which is the event of interest. This change of scale between many entering protons and few free protons and HA binding sites requires a different description than the previous continuous model.
To compute the mean time for HA conformation to change as the pH drops, we first extracted the forward and backward proton binding rates by converting the HA conformational change kinetics, obtained from experimental data at various pH  into rate constants.
At temperature T = 300K, when the pH decreases from 7 to 4, the number of protons bound to HA1 increases approximatively from 123 to 132 (Figure 3 in Huang et al. ), suggesting that the number of available number of binding site is ns = 9 at acidic pH. The Influenza virus carries nHA = 400 HA trimers  (Figure 3A) and thus there are exactly nHAns sites that can competitively bind protons. In this section, we compute the mean time that a threshold nT of bound protons to HA1 is reached, which is a model of fusogenic state, where proteins engage into the generation of a fusion pore with the endosomal membrane.
Figure 3. Free protons in the endosome triggers HA conformational change. (A) The right-hand side shows a scheme of an isolated HA trimer. Free protons in the endosome can bind to HA trimers. The protons binding rates r(X, c) and l(X) depend on the number of occupied sites X and on the concentration c of free protons in the endosome. When the number of bound protons reaches a given threshold, the HA trimer changes conformation into a fusogenic state. (B) MFPT of a discrete Markov jump process with transition rates r(X, c) and l(X) to a critical threshold nT = 6 is estimated with Monte-Carlo simulations (N = 1000) (solid red line), and compared with WKB approximation (solid black line) and experimental data points .
2.2.1. Modeling HA Conformational Change
To analyse the conformational change of a single HA trimer, we follow the occupied proton sites X(t, c) at time t, for a fix proton concentration c. During time t and t + Δt, the number of specific bound sites can either increase with a probability r(X, c)Δt, when a proton arrives to a free site or decreases with probability l(X, c)Δt when a proton unbinds or remains unchanged with probability 1 − l(X, c)Δt − r(X, c)Δt (Figure 3A).
We estimate hereafter the rates l(X, c) and r(X, c) and the critical threshold nT, by approximating the number of bound protons with the proton concentration c variable, by a linear function (Figure 3 in Huang et al. )
where is the mean number of bound protons at pH = 7 and
is the mean number of HA1 sites that are additionally protonated for a proton concentration c > 10−7mol. L−1. Recently, we also used a similar jump model  to study the conformational change of active proteins for the escape of non-enveloped viruses, based on the assumption that proton binding and unbinding rates were linear functions of the proton concentration. Here however, the mean number of bound sites depends linearly on the endosomal pH (log of the proton concentration) (Equation 16), confirming the non-linear binding and unbinding rates of protons to HA.
To account for the non-linearity of the mean number of bound protons (Equation 16), we derived the expressions of the binding r and unbinding l rates of protons to HA binding sites. First, we assume that the binding rate r(X, c) depends on both the proton concentration c and the number of free binding sites X, whereas the unbinding rate l(X) depends only on X. Indeed, an increased concentration of protons inside the endosome favors the encounter and binding between protons and HA sites, but do not influence the unbinding rate of bound protons. Moreover, we assume that the binding rate r(X, c) depends linearly on the proton concentration c and the number of free sites ns − X of the HA trimer leading to
where K is the forward binding rate of a proton to a binding site.
To determine the non-linear proton unbinding rate l(X, c), we use the mean number of protons bound to HA at different pHs (Equation 16). Using at equilibrium the concentration for which X0(c(X)) = X, the mass-action law leads to or equivalently , and we get
In summary, the binding and unbinding rates r and l are given by
2.2.2. Rate of HA Conformational Change
To compute the mean time that exactly nT protons are bound to a single HA, we use a Markov jump process for the number of protonated sites X(t, c) among the ns = 9 HA1 proton binding sites available. Following a similar approach as used in Lagache et al.  for non-enveloped viruses, we scale the variable
where ϵ = 1/ns and we use the Wentzel-Kramers-Brillouin (WKB) expansion of the mean first passage time (MFPT) τ(c) of the scaled number of protonated sites x(t, c) to the (unknown) critical threshold 0 < xT = ϵnT < ϵns [11, 14, 29–31 ]
where x0(c) is the mean number of HA1 sites that are additionally protonated for a concentration c > 10−7mol.L−1 (Equation 16). The function ϕ(x, c) is given by
Replacing the transition rates r(x, c) and l(x) by their expressions (19) in Equation (22), we obtain that (see Section 4)
To conclude, the expression for the MFPT (23) differs from the one computed in Lagache et al.  (Equation 5) derived for a linear unbinding rate l(x, c) = k−1xns. Equation (23) links the affinities between the ligand (concentration c) and the binding sites of a trimer to the conformational change mean time τ(c) of the trimer. The two unknown parameters of the conformational MFPT (Equation 23) are the binding rate K of protons to HA sites and the number nT of sites that have to be protonated (among the ns = 9 total sites) to trigger the HA change of conformation.
The reciprocal of the mean time has been measured for various pH values : (τ(pH = 4.9))−1 = 5.78s−1, (τ(pH = 5.1))−1 = 0.12s−1,…, (τ(pH = 5.6))−1 = 0.017s−1. To estimate the unknown K and nT, we thus use formula (23) to fit these data by a least square optimization procedure, and obtain that
and the forward rate
These two estimations are the predictions of the present model.
We reported here a good agreement between the WKB approximation (Equation 23) with the Monte-Carlo simulations of the Markov jump process, with the transition rates r(X, c) and l(X) given by Equation (19), and the experimental values of  (Figure 3B, model parameters are summarized in Table 2). The WKB solution is very close to the Markov jump simulations, especially for pH values ≥ 5.8, where the fusion takes place (Section 2.3 below). For lower values of the pH, the MFPT to threshold that triggers the conformational change of HA decreases drastically and a small discrepancy between discrete Monte-Carlo simulations and continuum WKB approximation can be seen. For these lower pH values, we observe that the WKB and the Monte-Carlo simulations agrees with the conformational change rates of HA, measured experimentally . We highlight that the fitting of only two parameters K and nT of the Markov model lead to a very good fit to all experimental data points, which indicates that the Markov jump approach with the WKB approximation is suitable to model the HA conformational changes.
2.2.3. A High Potential Barrier of HA Binding Sites Ensures HA Stability at Neutral pH
We have seen in Section 2.1 that during endosomal acidification, a huge number of protons enter the endosome (more than 20*106 that bind mostly to endosome buffers, leaving very few free protons (around 300 at pH 6)). To test whether HAs buffer entering protons or interact with the remaining few free protons, we estimate the potential barrier generated at each HA binding site. For this purpose, we compare the reciprocal of the forward rate constant K (Equation 25), which is the mean time for a proton to bind a HA protein, with the free Brownian diffusion time scale. For a fixed proton concentration at a value c, the proton binding time is , while the mean time for a proton to diffuse to the same binding site is [32–35]
The number of endosomal protons at concentration c is n(c) = NAcV, while η is the interacting radius between a proton and a binding site and Dp the diffusion constant of a free proton (Dp = 100μm2 s−1 measured in the cytoplasm ). For η = 1 nm, we find a small ratio
This result indicates that, on average, only 1 out of 104 encounters between a proton and a HA binding site lead to a binding event. Thus, the binding of protons to HA is strongly reaction-limited, dominated by a very high activation energy barrier at the HA binding sites. This high barrier prevents rapid proton binding and consequently, the buffering capacity of HAs can be neglected compared to the high capacity of other endosomal buffers. In addition, the activation energy barrier of HA binding sites ensures a high stability of the protein at pH above 6, as previously characterized in Table 2 of Krumbiegel et al.  and confirmed in Figure S1.
To conclude, we found that the threshold for HA1 conformational change occurs when there are nT = 6 bound proton in a total of ns = 9 binding sites. The binding is characterized by a very high potential barrier. Thus, when protons enter an endosome, they will first be captured by endosomal buffers. The remaining free protons can bind to HA1 sites after passing across the high potential barrier to trigger HA conformational change.
2.3. A Closed Model of Virus-Endosome Fusion
Combining the kinetic model of endosome acidification with the Markov jump model of HA conformational change, we now derive a kinetic model of HAs conformational change inside an endosome. We account for the nT = 6 protons activating a HA1 trigger leading to HA conformational change. We now estimate the numbers HA0(t), HA1(t)…HA6(t) of viral HAs that have 0, 1…6 bound protons at time t, and compute the number of fusogenic (active) HA6(t), responsible for membrane fusion. From relation (17), the forward rate of a proton to a free HA1 binding site is
and the backward rate l(X) is given by relation 19, thus the chemical Equations for protons Pe and HA proteins are summarized by
where the rate constant depends on each stage as given by relation 28. The stage HA6 is irreversible and the kinetic rate equations are
Once the proton entry rate (Equation 6) is known, these equations can be solved numerically.
2.3.1. Modeling the Onset of Fusion between Virus and Endosome Membranes
Membrane fusion is triggered by the conformational change of multiple adjacent trimers located in the contact zone between the viral and endosomal membranes [16, 17]. However, the exact number of fusogenic HAs involved in formation and fusion pore enlargement is still unclear. To estimate this number, we model the contact zone between the virus and endosome membranes by 120 HAs among 400 covering the viral particle  (Figure 4A). Thus, each trimer in the contact zone possesses 6 adjacent neighbors.
Figure 4. Model and fluorescence experiments of the intracellular onset of virus-endosome fusion. (A) The fusion between virus and endosome membranes is triggered by the conformational change of Na adjacent HAs in the contact zone between virus and endosome (= 120 among the 400 HAs covering the virus envelope ). (B) Solving Equation (31) we estimated the time window (95% confidence interval) of intracellular fusion for 1 ≤ Na ≤ 6. (C) Using time windows of fusion onset and endosome maturation kinetics (Equation 7), we estimated the localization (EE, ME, or LE) of fusion onset as function of the number Na of adjacent fusogenic HAs needed for the fusion onset.(D) MDCK cells expressing Rab5-CFP and Rab7-GFP were incubated with R18-labeled Influenza A viruses. Fusion was observed as a strong increase of R18 signal due to de-quenching after dilution. Scale bar = 1 μm (E) In vivo localization of fusion events. Fusion events were counted and categorized regarding their localization in early endosome (EE) (Rab5), ME (Rab5 + Rab7) or LE (Rab7).
We solved numerically Equation (31) and we chose the position of each new fusogenic HA (with nT = 6 bound sites) randomly among the 400 HAs covering the virus envelope. For each simulation, we defined the onset of virus endosome fusion by the activation of Na adjacent HAs in the contact zone (Figure 4A). Using 1, 000 Monte-Carlo simulations, we estimated the mean and confidence interval at 95% of the fusion onset time for different Na. We found that for Na = 1 − 2, most viruses fuse in EE, whereas for Na = 3−4, they fuse in ME. Finally, for Na = 5 or 6, viruses mostly fuse in LE (Figures 4B,C). To conclude, viruses shall fuse in ME [16, 17] and thus Na = 3 − 4.
2.3.2. Intracellular Localization of Fused Viral Particle with Live Cell Imaging
To experimentally determine the localization of virus fusion, we used the fluorescent endosomal markers Rab5 (EE) and Rab7 (LE) in combination with an intracellular fusion assay to detect virus-endosome fusion so that the localization to a specific compartment can be assigned. Single virus spots were analyzed, where fusion was indicated by a pronounced increase of spot signal (Figure S2). To determine the cellular localization of virus fusion, we analyse infected Rab5- and Rab7-expressing cells with R18-labeled viruses (Figure 4D). We classified single endosomes based on the presence of the two Rab proteins into three classes (Figure S3). Early endosomes (EE) do not show Rab7 association, such as late endosomes (LE) do not posses Rab5 signal. When endosomes possess both signals, they were counted as maturing endosomes (ME).
We observe a gradual increase of Rab7 along with a decrease of Rab5 (Figure 4D). After 5 min, we rarely observe fusion events in Rab5-only endosomes. The majority of fusion events (61%) are detected in maturing endosomes between 10 and 20 min post infection (Figure 4E). At later time points, the localization of fusion events shifted toward late endosomes. However, de-quenching kinetics show that fusion mostly occurs between 10 and 20 min (Figure S2).
We conclude that virual fusion was essentially associated with maturing endosomes confirming that the number of adjacent fusogenic HA required to mediate fusion are Na = 3 or 4.
3. Discussion and Conclusion
Influenza viruses are internalized into endosomes via receptor-mediated endocytosis. During their transport along microtubules, endosomes accumulate protons, which eventually enable virus-endosome fusion mediated by the influenza HA. This fusion mediates the release of the viral genome in the cell cytoplasm. The duration of endosomal transport as well as the localization of fusion critically depend on endosomal acidification and HA conformational change at low pH. Here we presented a new model to investigate the role of key parameters that shape the endosomal residence time of influenza viruses.
The Markov-jump process model of HA conformational change that we have developed here shows that 6 bound protons are enough to trigger conformational change for a total of 9 binding sites. The model also reveals that the HA activation is characterized by a high potential barrier with a forward rate constant K = 7, 500L.mol−1s−1. Interestingly the unbinding rate depends on the number of binding sites, suggesting a modification that depends on history. This is in contrast with the large number of protons of the order of millions that enters and contributes to acidification. This confirms that buffers play a critical role to reduce this enormous quantity of protons into a countable number that will trigger HA activation. The multiscale process involved here is very different from the previous analysis and results  we found for the Adeno-associated virus.
Finally, associating the kinetic model of endosomal acidification with a Markov-jump process model of HA conformational change, we estimated how the number of fusogenic HAs evolved in time inside endosomes, and we modeled the onset of fusion with the stochastic activation of Na adjacent HAs. Using the model, we predict a high HA thermal stability at neutral pH due to a high activation barrier of proton binding sites. For Na ≥ 3, we show that fusion should occur in ME, preventing a premature fusion in EE. Because endosomal maturation is associated with retrograde transport of endosomes along MTs, we conclude that mature virus should accumulate near the nuclear surface. Moreover, when an endosome contains multiple copies of Influenza viruses, we found that the cumulative buffering capacity of viruses might delay the acidification kinetics and viral escape into ME and LE.
4. Materials and Methods
We describe here the experimental methods for extracting parameters used to validate the modeling approach. In the last section, we present the computation for the WKB computation for formula Equation (23).
4.1. Materials, Cell and Virus Culture
Madine Darby Canine Kidney (MDCK) cells were cultured in Dulbeccos Modified Eagles Medium (DMEM) without phenol red, supplemented with 1% penicillin/streptomycin and 10% fetal calf serum (FCS). The cells were passaged every 3–4 days. One day prior to the experiment, the cells were detached from the cell culture flask using 0.5 % Trypsin/EDTA for about 10 min. The cells were diluted in DMEM and 10–50 * 105 cells were seeded in 35 mm poly-L-lysine coated glass bottom petri dishes (MatTek Corp.). Influenza A (H3N2) X-31 was propagated in chicken eggs and A/Panama/2007/99 on MDCK cells. Prior to the experiment the virus was diluted to 1 mg/ml protein concentration. Octadecylrhodamine B (R18) was purchased from Molecular Probes (Life Technologies, USA). Phosphate buffered saline (PBS) was used for all dilutions during the experiments. Double labeled FITC/Rhodamine dextran was purchased from Life Technologies (USA). FITC-dextran was purchased from Sigma-Aldrich, Germany. Rab5-GFP (in pCDNA3) and Rab7-GFP (pEGFP-C1) were kindly provided by Volker Haucke (FMP, Berlin). Rab7 was then cloned into pECFP. MDCK cells were transfected using Turbofect (Fermentas, USA) according to the manufacturers manual. To disrupt microtubules, the cells were incubated with medium containing 50 μM nocodazole (Sigma Aldrich, Germany) for 30 min before experiments.
MDCK were washed in PBS buffer and fixed in PBS containing 2% paraformaldehyde and 0.2 % glutaraldehyde for 20 min. The cells were permeabilized with PBS containing 0.2% Triton X-100 and 0.2 % BSA for 20 min, washed in PBS and incubated in anti-Nucleoprotein (Millipore, USA) antibody for 1 h. The cells were washed in PBS and incubated with the secondary anti-mouse Cy2 conjugate antibody for 1 h (Amersham, GE, USA). Finally, the cells were counterstained using PBS containing 0.2 μ g/ml DAPI for 10 min.
4.3. Fluorescence Microscopy
For fluorescence microscopy, we used an Olympus FV1000-MPE confocal microscope (Olympus, Japan) equipped with 405 nm (DAPI), 440 nm (CFP), 488 nm (GFP), 559 nm (R18) and 635 nm (A647) laser lines, an Olympus 60x/1.2 water UPlanSApo objective and 405-458/515/559/635 405/488/559/635 dichroic mirror filter sets.
4.4. Endosomal pH Determination
One day prior to the experiment, MDCK cells were seeded into 35 mm poly-L-lysine coated glass bottom petri dishes (MatTek Corp.). For dextran labeling, the cells were washed with PBS and incubated in serum free medium for 30 min at 37° C, followed by 5 min with 10 mg/ml dextran at 37° C (pulse). After the pulse, the cells were immediately washed and image acquisition was started. For plotting the pH evolution (experiments and model), we considered a time delay of 5 min due to technical limitations and set the starting point to pH 7.2.
4.5. Determination of a pH Standard Curve
The pH standard curve for intracellular pH measurements was done as previously described by others (). In short, MDCK cells were detached from the culture dish, washed and pelleted with 2,000 g for 5 min. The cells were divided into eight fractions, pelleted again and resuspended pH standard buffer (obtained by mixing 50 mM HEPES buffer with 50 mM MES buffer (both containing 50 mM NaCl, 30 mM ammonium acetate, 40 mM sodium azide and 10 μM nigerizin). The samples were left on ice for 5 min and analyzed by flow cytometry.
4.6. Calculation of the Proton Binding Capacity of Viral Proteins
Proton binding and the total charge of the proteins was calculated as follows. First, the pKa values of all titrable residues in the proteins were determined with electrostatic energy calculations using the software Karlsberg+ . The calculations are based on the crystal structures with PDB IDs 1hgg, 1hgd, 2bat and 2q06. All non-protein molecules were removed from the structures and therefore not considered in the calculations. Except for hydrogen atoms, all atomic coordinates were kept as found in the crystal structure in the Karlsberg+ calculations. If a titratable residue was missing in the structure, it was assumed that its pKa value equals the model pKa of this residue. The model pKa is the experimental pKa value of a residue in aqueous solution. The protonation state and therefore the charge of a residue at a certain pH value are determined by its pKa value. Here a residue is considered to be protonated if its pKa value is larger than the corresponding pH value. The number of bound protons at a given pH value was then obtained by counting the protonated titratable residues in a protein, the total charge by summing up all individual residue charges.
4.7. Influenza Virus—Ghost Membrane Fusion Assay
To check that Influenza virus are fusion competent and calibrate intracellular fusion experiment we first performed an in vitro ghost-membrane assay . We labeled Influenza viruses with lipophilic dye R18 and measured virus-membrane fusion by monitoring the fluorescence de-quenching (FDQ) of the lipid-like fluorophore R18 upon fusion of R18-labeled viruses with membranes. To this end, 10 μl of labeled virus suspension (1 mg/ml) were mixed with 40 μl ghost suspension (≈ 2*105 cells) and incubated for 20 min at RT. Unbound virus was removed by centrifugation (5 min, 1,200 g). The virus-ghost suspension was transferred to a glass cuvette containing pre-warmed fusion buffer (pH 7.4), and the fluorescence was detected (λex = 560 nm; λem = 590 nm) by using a Horiba Yobin Yvon FluoroMax spectrofluorometer. Fusion was triggered by the addition of citric acid (0.2 M). The suspension was stirred continuously with a 2 by 8 mm Teflon-coated magnetic stirring bar. After 600 s the fusion was stopped by adding Triton X-100 (50 μl, final concentration 0.5%) to obtain maximum R18 de-quenching. The final pH in the cuvette was measured using a standard pH meter. The percentage of FDQ was calculated as:
where F(0) and F(t) are the fluorescence intensity before starting fusion and at a given time (t), respectively. We observed that fusion starts around pH 6 as predicted by our model and already shown experimentally  (Figure S1).
4.8. Intracellular Fusion Assay
Influenza virus was diluted to a final concentration of 1 mg/ml in PBS and incubated with 20 mM R18 for 30 min at room temperature. Unbound R18 was removed by centrifugation at 25.000 g for 5 min or gel filtration (G25 sephadex in PBS). The virus was resuspended/eluted in PBS. Immediately before the experiment, the virus was diluted to 40 μg/ml and viral aggregates were removed with a 0.2 μm sterile filter. The virus was applied to the cells and allowed to bind for 10 min at 4°C. The temperature was elevated to 37°C and the R18 fluorescence was monitored by confocal microscopy. 40 min after start of acquisition, the cells were fixed with PBS containing 2% paraformaldehyde and 0.2 % glutaraldehyde and the DNA was stained with Hoechst 33342 (Invitrogen, USA). The boundaries between cells were determined from the bright field image. Summed z-stacks were analyzed using an IDL-based particle identification software .
4.9. Derivation of the Acidification Kinetics Equation (3)
Protons that enter inside the endosome interact with endosomal buffers, and the number of free protons Pe(t) inside an endosome is given by the mass-action law
where Pe-Bi(t) and Bi(t) are the number of weak acids and bases inside the endosome at time t, NA is the Avogadro constant and Ve is the volume of the endosome. Assuming that the membrane potential Ψ(t) reaches rapidly its steady state value Ψ(∞) compared to the acidification kinetics , we approximate the pumping rate λ(t)S with its steady state value λ
where the parameter λ is related to the membrane potential Ψ(∞). In addition, the protons leak Lext(t) is proportional to the endosomal concentration and the endosomal surface 
where L is a permeability constant. Consequently, using approximations 33 and 34 in Equation (32), we obtain the general dynamics of free protons:
When the protons enter the endosome they can interact with buffers which is much faster than acidification. Thus, we assume that ΔPe(t) protons entering the endosome during dt bind instantaneously to bases, which leads to step decrease −ΔBi(t) of the number of proton binding sites of each base Bi(t):
To estimate the pH change ΔpH(t) associated with the entry ΔPe(t) of protons and the corresponding decrease −ΔBi(t) of each base, we use Equation (2) at equilibrium
where and Ci = Pe-Bi(0) + Bi(0) are constant. Consequently,
where pKi = −log(Ki)/log(10). By differentiating Equation (39) with respect to each Bi(t), we obtain the infinitesimal variation ΔpH(t) of the endosomal pH at time t
Using Equation (38), we get
is the buffering capacity of the weak acid-base couple (Pe-Bi, Bi). Finally, using Equations (36) and (42), the infinitesimal change ΔpH(t) of the endosomal pH associated with the entry of ΔPe(t) protons is
where is the total buffering capacity of the endosome, which is approximately constant . Finally, using the proton extrusion and pumping rates (Equations 33 and 34), we obtain the kinetic Equation
The endosomal pH is related to the number of free protons Pe(t) by
and replacing the pH derivative in Equation (46), we obtain that the accumulation of free protons Pe(t) inside the ensosome during acidification is given by the kinetics Equation
4.10. Estimation of the Viral Buffering Capacity (Equation 13)
First, we estimated the buffering capacity of NP proteins  by computing the pKa values of all titratable residues in the proteins with electrostatic energy calculations using the software Karlsberg+ . We then determined the mean number of protonated residues of NP proteins and we found that increases almost linearly with pH:
indicating that the buffering capacity of NPs is approximatively constant between pH 7 and 5 (Equation 42)
where is the volume of the viral internal lumen, for a spherical viral particle with radius rv = 60 nm . The structure of the matrix M1 protein is unknown and consequently, we use the cumulative contributions of Asp, Glu and His residues to estimate the number of M1 proton binding sites. We thus estimate the fraction Pi(pH) of occupied residues for a fixed pH using the equilibrium constant pKai for any residue i (Asp, Glu or His) to be
The mean number of protonated site is then given by
where the number of residue for each group is , and . Using Equation (53), we plotted as function of the pH and observed that is almost a linear function
and obtain that
Additionally to internal M1s and NPs proteins, protons entering the viral core through M2 channels can also bind to viral nucleic acids and in particular to basic groups in the guanine, adenine and cytosine nucleotides . In particular, the buffering capacity βRNA of oligonucleotides in solution, for a concentration cmonomers of monomers, has been estimated to be in the pH range 5–7 (Figure 3D in Stoyanov and Righetti, ). Consequently the buffering capacity of the = 12, 000 viral nucleotides  is approximatively equal to
Finally, the viral core lumen should also contain other unspecific buffers such as cytoplasmic buffers enclosed during the viral assembly, leading to an unspecific buffering capacity inside the viral lumen that has to be added to the buffering capacities and of internal proteins. Due to possible ionic exchange between viral and endosomal lumens, we approximate with the endosomal buffering capacity , which is independent of the pH and has been estimated to be 
4.11. Derivation of Equation 23
Using the WKB approach , we obtain that the mean first passage time (MFPT) τ(c) of the scaled number of protonated sites x(t, c) to the (unknown) critical threshold 0 < xT = ϵnT < ϵns is given by [11, 29–31]
where x0(c) is the mean number of HA1 sites that are additionally protonated for a concentration c > 10−7mol.L−1 (Equation 16) and ϕ(x, c) is given by
We first compute
Using expressions for the binding and unbinding rates 19, we get
which reduces to
Finally, re-injecting expression of potential ϕ(xT, c) in the MFPT τ(c), we obtain that
Using ϵ = 1/ns and xT = nT/ns, we get
TL, AH, and DH: Designed research. TL and CS: Performed experimental and simulations works. TM: Contributed analytic tools. TL, CS, and DH Analyzed data. TL, CS, AH, and DH: Wrote the paper.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
This research was supported by a Marie Curie grant (DH), by the Deutsche Forschungsgemeinschaft (HE 3763/15-1) (AH) and the Bundesministerium fur Bildung und Forschung (eBio: ViroSign) (CS and AH). TL is partially funded by a fellowship from the Fondation pour la Recherche Medicale and a grant from the Philippe Foundation.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/article/10.3389/fphy.2017.00025/full#supplementary-material
3. Mistry B, D'Orsogna MR, Webb NE, Lee B, Chou T. Quantifying the Sensitivity of HIV-1 Viral Entry to Receptor and Coreceptor Expression. J Phys Chem B (2016) 120:6189–99. doi: 10.1021/acs.jpcb.6b02102
7. Dauty E, Verkman AS. Actin cytoskeleton as the principal determinant of size-dependent DNA mobility in cytoplasm: a new barrier for non-viral gene delivery. J Biol Chem. (2005) 280:7823–8. doi: 10.1074/jbc.M412374200
8. Xiao PJ, Samulski RJ. Cytoplasmic trafficking, endosomal escape, and perinuclear accumulation of adeno-associated virus type 2 particles are facilitated by microtubule network. J Virol. (2012) 86:10462–73. doi: 10.1128/JVI.00935-12
10. Schelker M, Mair CM, Jolmes F, Welke RW, Klipp E, Herrmann A, et al. Viral RNA Degradation and Diffusion Act as a Bottleneck for the Influenza A Virus Infection Efficiency. PLoS Comput Biol. (2016) 12:e1005075. doi: 10.1371/journal.pcbi.1005075
13. Huang Q, Sivaramakrishna RP, Ludwig K, Korte T, Böttcher C, Herrmann A. Early steps of the conformational change of influenza virus hemagglutinin to a fusion active state: stability and energetics of the hemagglutinin. Biochim Biophys Acta (2003) 1614:3–13. doi: 10.1016/S0005-2736(03)00158-5
16. Danieli T, Pelletier SL, Henis YI, White JM. Membrane fusion mediated by the influenza virus hemagglutinin requires the concerted action of at least three hemagglutinin trimers. J Cell Biol. (1996) 133:559–69. doi: 10.1083/jcb.133.3.559
17. Ivanovic T, Choi JL, Whelan SP, van Oijen AM, Harrison SC. Influenza-virus membrane fusion by cooperative fold-back of stochastically induced hemagglutinin intermediates. Elife (2013) 2:e00333. doi: 10.7554/eLife.00333
22. Zaraket H, Bridges OA, Duan S, Baranovich T, Yoon SW, Reed ML, et al. Increased acid stability of the hemagglutinin protein enhances H5N1 influenza virus growth in the upper respiratory tract but is insufficient for transmission in ferrets. J Virol. (2013) 87:9911–22. doi: 10.1128/JVI.01175-13
23. Bayer N, Schober D, Prchla E, Murphy RF, Blaas D, Fuchs R. Effect of bafilomycin A1 and nocodazole on endocytic transport in HeLa cells: implications for viral uncoating and infection. J Virol. (1998) 72:9645–55.
24. Leiding T, Wang J, Martinsson J, DeGrado WF, Arsköld SP. Proton and cation transport activity of the M2 proton channel from influenza A virus. Proc Natl Acad Sci USA. (2010) 107:15409–14. doi: 10.1073/pnas.1009997107
25. Lamb R, Krug R. Orthomyxoviridae: The viruses and replication. In: Knipe D, Howley P, Griffin D, editors. Fields Virology. 4th Edn. Philadelphia, PA: Lippincott Wiliams and Wilkins (1996) 1487–1531.
26. Stoyanov AV, Righetti PG. Buffer properties of biopolymer solutions, as related to their behaviour in electrokinetic methodologies. J Chromatogr A. (1999) 838:11–8. doi: 10.1016/S0021-9673(99)00090-4
28. Krumbiegel M, Herrmann A, Blumenthal R. Kinetics of the low pH-induced conformational changes and fusogenic activity of influenza hemagglutinin. Biophys J. (1994) 67:2355–60. doi: 10.1016/S0006-3495(94)80721-0
30. Matkowsky B, Schuss Z, Knessl C, Tier C, Mangel M. Asymptotic solution of the Kramers-Moyal equation and first-passage times for Markov jump processes. Phys Rev A. (1984) 29:3359–69. doi: 10.1103/PhysRevA.29.3359
41. Ng AKL, Zhang H, Tan K, Li Z, Liu Jh, Chan PKS, et al. Structure of the influenza virus A H5N1 nucleoprotein: implications for RNA binding, oligomerization, and vaccine design. FASEB J. (2008) 22:3638–47. doi: 10.1096/fj.08-112110
Keywords: modeling, first passage time, asymptotic analysis, conformational change, endosomal acidification, influenza virus, trafficking, Kramers-Moyal approximation
Citation: Lagache T, Sieben C, Meyer T, Herrmann A and Holcman D (2017) Stochastic Model of Acidification, Activation of Hemagglutinin and Escape of Influenza Viruses from an Endosome. Front. Phys. 5:25. doi: 10.3389/fphy.2017.00025
Received: 30 March 2017; Accepted: 06 June 2017;
Published: 23 June 2017.
Edited by:Mario Nicodemi, Universita' di Napoli “Federico II”, Italy
Reviewed by:Luis Diambra, National University of La Plata, Argentina
Haiguang Liu, Beijing Computational Science Research Center, China
Copyright © 2017 Lagache, Sieben, Meyer, Herrmann and Holcman. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: David Holcman, email@example.com
†Present Address: Thibault Lagache, Department of Biological Sciences, Columbia University, New York, NY, United States
Christian Sieben, School of Basic Sciences, École Polytechnique Fédérale de Lausanne, Lausanne, Switzerland