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 ΔP_e(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 [20]:

where k_i (resp. ki(-1)) are binding (resp. unbinding) rate constants of protons to weak bases B_i, 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 βe0 of the endosome. The kinetics equation for the number of free protons P_e(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_∞ [12] and the permeability L decreases with the endosomal maturation [19] as they depend on time, as we analyse below.

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 [21]. We approximate here the kinetics of the ratio Rab5/Rab7 (Figure 4C in Rink et al. [21]) by a sigmoidal function

with the two free parameters: the half-maturation t_1/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

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. [21]) 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 [22]. The steady-state pH is pH∞early=6.0 and pH∞late=5.5 for early and late endosomes respectively [23]. 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 Learly=3.510-3NAcm s-1 and Llate=0.1NAcm s-1, respectively, and the half-maturation time is t_1/2 = 10 min.

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 [24], summarized in the chemical equation

When a proton P_e binds a free M2 protein channel with rates k_e(binding) and ke-1 (unbinding), it is transported inside the virus core with a rate kv-1, while exit occurs with a rate k_v. At steady state, the inward flux in a single virus is computed from Equation (10) (see [24])

where n_M2 is the number of M2 channels per viral particle, P_v 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 [24]. The most abundant internal proteins are M1 (3, 000 copies per virus) and the nucleoproteins (NP, 330 copies per virus) [25] (Figure 1A). Proton binding sites of viral proteins are the ionogenic groups in their amino acid side chains [26], 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) [26]. 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 βv0 is the buffering capacity of the lumen inside the virus, and βvM1, βvNP and βvRNA 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 P_v(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 Pe(t=0)=10-7.2NAVe and Pv(t=0)=10-7.2NAVv, 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 r_e = 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 r_e = 100 nm. In addition viral particles may accumulate during the endosomal journey [27]. 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.

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 n_T, by approximating the number of bound protons X~0(c) with the proton concentration c variable, by a linear function (Figure 3 in Huang et al. [15])

where X~0(10-7mol.L-1) 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⁻⁷mol. L⁻¹. Recently, we also used a similar jump model [11] 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 n_s − 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 c(X)=103Xns-7 for which X₀(c(X)) = X, the mass-action law leads to l(X0(c),c)r(X0(c),c)=1 or equivalently l(X)Kc(X)(ns-X)=1, and we get

In summary, the binding and unbinding rates r and l are given by

To compute the mean time that exactly n_T protons are bound to a single HA, we use a Markov jump process for the number of protonated sites X(t, c) among the n_s = 9 HA1 proton binding sites available. Following a similar approach as used in Lagache et al. [11] for non-enveloped viruses, we scale the variable

where ϵ = 1/n_s 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 < x_T = ϵn_T < ϵn_s [11, 14, 29–31 ]

where x₀(c) is the mean number of HA1 sites that are additionally protonated for a concentration c > 10⁻⁷mol.L⁻¹ (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)

where F(x)=32log(10)x2-log(107c)x.

To conclude, the expression for the MFPT (23) differs from the one computed in Lagache et al. [11] (Equation 5) derived for a linear unbinding rate l(x, c) = k₋₁xn_s. 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 n_T of sites that have to be protonated (among the n_s = 9 total sites) to trigger the HA change of conformation.

The reciprocal of the mean time 1τ(c) has been measured for various pH values [28]: (τ(pH = 4.9))⁻¹ = 5.78s⁻¹, (τ(pH = 5.1))⁻¹ = 0.12s⁻¹,…, (τ(pH = 5.6))⁻¹ = 0.017s⁻¹. To estimate the unknown K and n_T, 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 [28] (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 [28]. We highlight that the fitting of only two parameters K and n_T 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.

We have seen in Section 2.1 that during endosomal acidification, a huge number of protons enter the endosome (more than 20^*10⁶ 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 τbind=1Kc, 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) = N_AcV, while η is the interacting radius between a proton and a binding site and D_p the diffusion constant of a free proton (D_p = 100μm² s⁻¹ measured in the cytoplasm [36]). For η = 1 nm, we find a small ratio

This result indicates that, on average, only 1 out of 10⁴ 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. [28] and confirmed in Figure S1.

To conclude, we found that the threshold for HA1 conformational change occurs when there are n_T = 6 bound proton in a total of n_s = 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.

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 [17] (Figure 4A). Thus, each trimer in the contact zone possesses 6 adjacent neighbors.

We solved numerically Equation (31) and we chose the position of each new fusogenic HA (with n_T = 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.

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.