Research ArticleVIROLOGY

Atomistic dynamics of a viral infection process: Release of membrane lytic peptides from a non-enveloped virus

See allHide authors and affiliations

Science Advances  14 Apr 2021:
Vol. 7, no. 16, eabe1761
DOI: 10.1126/sciadv.abe1761


Molecular simulations have played an instrumental role in uncovering the structural dynamics and physical properties of virus capsids. In this work, we move beyond equilibrium physicochemical characterization of a virus system to study a stage of the infection process that is required for viral proliferation. Despite many biochemical and functional studies, the molecular mechanism of host cell entry by non-enveloped viruses remains largely unresolved. Flock House virus (FHV) is a model system for non-enveloped viruses and is the subject of the current study. FHV infects through the acid-dependent endocytic pathway, where low pH triggers externalization of membrane-disrupting (γ) peptides from the capsid interior. Using all-atom equilibrium and enhanced sampling simulations, the mechanism and energetics of γ peptide liberation and the effect of pH on this process are investigated. Our computations agree with experimental findings and reveal nanoscopic details regarding the pH control mechanism, which are not readily accessible in experiments.


Understanding the molecular mechanism by which viruses invade host cells can provide insights for developing antiviral therapies and drug treatments for diseases caused by these viruses. Enveloped viruses such as HIV, Ebola, and influenza breach host cell membranes through a membrane fusion process catalyzed by fusion proteins in the viral envelope, and this process has become increasingly well characterized (1, 2). However, the mechanisms of cellular membrane penetration by non-enveloped viruses remain enigmatic. Emerging studies suggest that several of these non-enveloped viral families, including Picornaviridae, Polyomaviridae, Reoviridae, Parvoviridae, Adenoviridae, and Nodaviridae, contain a membrane-active component of their capsid, which becomes active under endosomal conditions (36).

Flock House virus (FHV), an insect-infecting member of the family Nodaviridae, serves as a model system for investigating the molecular details of cell entry in non-enveloped viruses (7). FHV can replicate in insects (8), plants (9), and yeast cells (10), and the viral genome consists of two single-stranded positive-sense RNA genomes. The capsid proteins initially assemble into an immature T = 3 icosahedral capsid, composed of 60 identical icosahedral asymmetric units (iASUs), and each iASU consists of three chemically identical α subunits, which have 407 residues. The subunits are arranged in three slightly different quasi-equivalent positions (A, B, and C). The capsid undergoes maturation to an infectious particle involving autoproteolytic cleavage of the capsid proteins between residues Asn363 and Ala364, catalyzed by Asp75. This results in two particle-associated cleavage products: a mature capsid protein β (residues 1 to 363) and a 44-residue γ peptide (residues 364 to 407). The cleavage site (N363/A364) is on the interior of the capsid shell, and mutations to either the catalytic D75 residue or the cleaved residue N363 block the maturation cleavage and render the particles non-infectious (11). The implication is that a dissociated γ peptide is essential for viral infection.

The x-ray crystal structure of the mature capsid was resolved at 2.7-Å resolution and shows that γ peptides are noncovalently associated with the capsid interior (Fig. 1). While the N-terminal residues (364 to 385) of γ form an amphipathic α-helix, the predominantly hydrophobic C-terminal region of γ (residues 385 to 407) is not resolved in the crystal structure, presumably due to a disordered structure and/or conformational heterogeneity across the capsid (12). Biophysical studies have shown that the N-terminal helical region (termed γ1) is capable of membrane disruption (1315), while the C-terminal region of γ plays an important role for specific recognition of the FHV genome and packaging of genomic RNA during assembly (16). Subsequent studies showed that the full-length peptide has significantly greater membrane lytic activity than γ1, and mutations and deletions in the C-terminal region could impair the lytic activity and infectivity (17, 18). The positions of the γ helices are located in significantly different environments in the three subunits of the iASU due to quasi-equivalence. The A subunits form pentamers, and the γ peptides from the A subunits (γA) interact with each other and form a pentameric helical bundle, with the helical axes oriented normal to the capsid surface (Fig. 1C). This organization is stabilized by hydrophilic interactions between peptides and are positioned far from the viral RNA (19). In contrast, γ peptides from B and C subunits are adjacent to the resolved viral RNA segments and form numerous contacts with the genome (12, 19). On the basis of the orientation of the γA peptides and lack of RNA contacts, it has been speculated that the pentameric helical bundle of γA peptides is responsible for membrane lysis (19, 20). This model is supported by studies on Nudaurelia capensis ω virus (NωV; T = 4), where γ peptides are positioned in similar environments (quasi-equivalent) (21). Studies with NωV using liposome-based assays show that maximal membrane disruption is acquired with less than half of the subunits of the capsid cleaved (21), whereas structural studies using time-resolved cryo–electron microscopy (cryo-EM) suggest that subunits surrounding the icosahedral fivefold axis (A subunits) are cleaved much earlier than the subunits positioned in different quasi-equivalent environment of the capsid (22).

Fig. 1 FHV structure.

(A) T = 3 capsid structure. (B) Asymmetric unit. (C) Interior view of fivefold vertex. The A proteins are blue, and the yellow helices are the γ peptides. All images are generated from PDB ID: 4FTB.

While the FHV capsid crystal structure shows that γ peptides are sequestered inside the capsid (12), it is believed that the role of the γ peptides in the infection process is to disrupt endosomal membranes and allow the capsid to escape these membrane-bound compartments. The proposed model is that low pH induces conformational and/or dynamic changes to the capsid, which allow or promote γ peptides to be externalized from the capsid. This model is supported by several lines of evidence. (i) Mass spectrometry combined with limited proteolysis suggests “occasional breathing” or transient exposure of γ peptides on the viral surface (23). (ii) Liposome dye leakage assays showed a pH-dependent lytic activity of FHV capsids, with maximal leakage at pH 6.0, consistent with the pH of endosomes (20). (iii) FHV infectivity was also correlated to endosomal pH, as increasing the endosomal pH through the addition of bafilomycin A1, an H+ adenosine triphosphatase (ATPase) inhibitor, decreased infectivity (20). (iv) When liposome dye leakage assays were performed at pH 7.0 following preincubation at pH 6.0, for as little as 1 hour, the dye leakage activity was equivalent to when the leakage assay was performed at pH 6.0 (20). (v) Bis-ANS (4,4′-Dianilino-1,1′-binaphthyl-5,5′-disulfonic acid), a fluorophore that binds aromatic and hydrophobic residues, had maximal binding at pH 6.0 (20).

While these previous findings provide strong support that FHV has an acid-dependent membrane disruption mechanism involving capsid conformational changes, the nature of the structural transitions and molecular details are not clearly understood. All-atom molecular dynamics (MD) simulations, sometimes termed a “computational microscope” (24, 25), have emerged as a powerful tool for studying viral systems, as this technique can provide substantial insight into the structure and dynamics at a high spatial and temporal resolution, which is typically inaccessible in experiments. Furthermore, MD simulations can go beyond observation of structural features, and thermodynamic quantities can be obtained, which can allow a structure-function model to be developed or validated. With ever increasing computational power, several simulation studies have been performed on full virus capsids to investigate the capsid dynamics and physicochemical properties (2629). In this present study, we use fully atomistic MD simulation to investigate the details of a stage of a viral infection process: the liberation of γ peptides from the interior of a non-enveloped capsid and the effects of pH. Previous simulation studies have examined viral life cycle processes, namely, assembly (3035) and maturation (3639), using various forms of coarse-grained models and/or implicit solvent models. We believe that the current study provides a major advance in the field by examining a viral infection process and its extensive use of atomistic simulations of a large portion of the capsid in explicit solvent through both equilibrium and enhanced sampling simulation methods.


pH effects on the structure and equilibrium dynamics of FHV capsid

In our study, we begin by analyzing equilibrium simulations of the FHV pentameric capsid system (Fig. 2A) at neutral and low pH. These systems were simulated in triplicate for 500 ns each (see Table 1). The systems were relatively stable [as measured by Cα root mean square deviation (RMSD) and root mean square fluctuation (RMSF)] during this time scale, and similar levels of structural relaxation were observed across the independent trials for each pH (figs. S1 and S2). The RMSD differences between the neutral and low pH systems are modest, but the low pH system displays a statically significant greater deviation than neutral pH. The averaged RMSD values over the final 300 ns and averaged over the three trials are 0.385 (±0.05) nm and 0.411 (±0.15) nm for neutral and low pH, respectively. We next analyzed the exposure of hydrophobic residues on the capsid surface through solvent-accessible surface area (SASA) calculations. It is experimentally known from bis-ANS fluorescence measurements that low pH induces increased hydrophobic exposure (20); therefore, we wanted to observe whether our simulations captured this feature. In Fig. 3A, the distributions of SASAs of hydrophobic residues are presented, and it is observed that, at low pH, there is an increase in hydrophobic exposure.

Fig. 2 FHV simulation system and titratable residues.

(A) The simulation system consists of five asymmetric units surrounding a fivefold symmetry axis. The γ peptides associated with the A subunits are retained and can be faintly seen in yellow beneath the fivefold vertex. (B and C) Location of residues that change protonation states in low pH model. The residues consist of D75 (cyan), H215 (yellow), and H334 (orange), which are shown in van der Waal representation. Both a single asymmetric unit (B) and two adjacent asymmetric units (C) are shown to illustrate the location of titrated residues related to neighboring proteins. The H215 residues lie at the center of the asymmetric unit (quasi-threefold symmetric axis), while H334 and D75 lie at the interface between asymmetric units, which is highlighted by the dashed black line. All images are generated externally to the capsid perspective (looking down on capsid) and are from PDB ID: 4FTB.

Table 1 List of simulation types and lengths.

View this table:
Fig. 3 Structural changes induced by low pH.

(A) Exposure of solvent-accessible hydrophobic residues at neutral and low pH. (B to E) The pore radius at the fivefold axis is measured for the replicates at both neutral (B) and low pH (D). The structure with the maximum pore radius is shown for neutral (C) and low pH (E) simulations.

While the SASA and RMSD calculations indicate increased conformational changes in the capsid at low pH, we wanted to further connect structural changes to potential functional motions of the capsid. Based on the structural organization of the γ peptides, it has been proposed in the literature that γs will externalize through the fivefold symmetry axis (19, 40), and it has been shown in NωV that pentameric γs are primarily responsible for membrane disruption (21). Therefore, we examined the dynamics of the capsid at the fivefold symmetry axis to evaluate if low pH induces increased dynamics in the fivefold vertex region. The pore radius at the fivefold vertex was measured during the simulations, and the profiles (Fig. 3, B and D) show a consistent increase in pore radius at low pH, which is consistent across the multiple trials. The maximum observed pore radius at low pH is 7.3 Å (Fig. 3E), which is more than double the maximum pore radius at neutral pH (3.5 Å; Fig. 3C). Furthermore, averaging over the last 300 ns and over the replicates shows that the average pore radius at low pH is 4.0 (±0.6) Å, while the neutral pH simulations display an average pore radius of only 1.3 (±0.2) Å. Representative equilibrium trajectories are visualized in movies S1 (neutral pH) and S2 (low pH) to demonstrate the fivefold pore dynamics. These pore measurements reveal a potential γ liberation pathway through the fivefold axis, and this pathway is consistent with the previous hypothesis based on the FHV structure-based model.

While we observe structural changes at the fivefold axis, which we believe are functionally important, the residues that have undergone protonation at low pH are distal to the fivefold vertex (Fig. 2B). We evaluated fivefold pore dynamics in two additional systems, one in which only D75 was protonated and one in which H215 and H334 were protonated. In the D75 protonated system, the pore radius did not exceed 3 Å, whereas the system with histidines protonated reaches a pore radius of 5.3 Å (fig. S3), consistent with the low pH systems. This result supports titration of histidine residues as being critical for altering capsid dynamics. Furthermore, to identify correlations between the titratable residues and capsid functional motions, we examined the structural dynamics of the systems by calculating secondary structural (SS) changes and linear mutual information (LMI). In Fig. 4A, we present the change in SS of the A subunits relative to the neutral pH capsid, collected over the last 300 ns of the simulation trajectories. In Fig. 4 (B and C), we show the location of residues with significant (>10%) changes in SS. The residues that display significant changes in α-helix or β-sheet content are 91 to 97, 169 to 171, 178, 203 to 204, 210 to 211, 244, 328 to 332, 334 to 335, 352 to 353, and 357 to 358. We see that these residues are proximal to the two histidine residues (H215 and H334) and participate in the protein-protein interfacial contacts. These subtle conformational changes could have a profound effect on the collective dynamics of the capsid. Furthermore, it appears that there is a connection of residues emanating from H334, which displays differential SS changes, and propagating toward the fivefold vertex (Fig. 4C). We did not find any change in SS in the neighboring regions of the other protonated residue, D75.

Fig. 4 pH effects and dynamic couplings.

(A) Change in SS in the A subunit in going from neutral to low pH. (B and C) The Cα atoms of residues with SS changes greater than 10% are displayed as van der Waal spheres on the A subunit structure. Color scale is from red (largest decrease) to blue (largest increase). Top down (B) and internal (C) views are presented. The black star indicates the location of the fivefold vertex, and the orange dashed line is drawn to illustrate a pathway from H334 to the fivefold vertex. Residues that undergo protonation changes are labeled. (D and E) LMI at neutral (D) and low pH (E). The correlations with less than 0.4 are set to zero to highlight stronger correlations. The blue vertical boxed regions are surrounding the titrating residues (D74, H215, and H334), and the red horizontal box is around residues 178 to 186, which is a loop at the fivefold vertex. (F) The boxed red region in (D) and (E) corresponds to the red loop shown on the pentameric structure.

We performed LMI analysis to characterize the correlations inside the capsid, and the corresponding matrices are shown in Fig. 4 (D and E) (41, 42). In general, more mutual information is observed across many different regions of capsid at neutral pH. There is significant mutual information between the innermost (“gatekeeper”) loop of the fivefold axis with different regions of the capsid, but these correlations are diminished at low pH. The region around H215 also shows strong correlations across the capsid, including with the gatekeeper loop, and again, most of those correlations are diminished at low pH. H215 is located at the center of the asymmetric unit (quasi-threefold axis; Fig. 2B) and may serve as a hub residue capable of modulating protein-protein dynamics throughout the capsid. Our results lead us to hypothesize that, under low pH conditions, the capsid motion of the fivefold axis is not concerted with the other regions of the capsid, which we suspect may lead to pore formation at the fivefold vertex. Our interpretation of the LMI and SS analyses is that the capsid is more connected and possibly more rigid at neutral pH, and at low pH, the protonation of the histidine residues leads to a disruption of the network. In particular, these results point toward a model in which, under low pH conditions, the capsid motion of the fivefold axis is not concerted with the other regions of the capsid, which we suspect leads to weaker/softer interactions at the fivefold vertex that is supported by the pore analysis (Fig. 3, B to E) and further investigated in the following sections.

γ Liberation pathways from nonequilibrium simulations

To gain insight into the molecular details of γ liberation, we performed biased [steered MD (SMD)] nonequilibrium simulations to generate translocation of the γ peptides from the interior (under) of the capsid to the exterior (above). We performed the biasing in two separate approaches. In one approach (referred to as “coupled CV”), the collective center of mass (COM) of the five peptides was biased; in the second approach (“single-peptide CV”), the COM of a single peptide was biased. A representative trial for each biasing method is presented in movies S3 (coupled CV) and S4 (single-peptide CV) for the low pH system. Coupled CV was partially motivated by speculation in the literature that the five peptides would move as a unit/bundle and externalize simultaneously (19, 40). Using coupled CV was an attempt to evaluate the bundle hypothesis, where if interpeptide interactions were sufficiently strong the peptides would move collectively (stay as a bundle) and not separate. We conducted 10 independent pulling simulations for each system (neutral and low pH capsid), and all simulation trajectories show that γ peptides were not liberated simultaneously. This finding contradicts the previous bundle hypothesis. We have calculated the nonequilibrium work done by the pulling force required for liberation of the first γ peptide using coupled CV (fig. S4, A and B). We observe that the low-pH environment reduces the work required for externalization in a statistically significant manner (fig. S4B). We also examined the maximum force observed during the SMD trials but did not find that the pH significantly affected this quantity (fig. S4C). For the trials at each pH that displayed the lowest work values, we continued the externalization until all five peptides were externalized. For this calculation, we removed the γ peptides from the system once they were liberated from the capsid and then restarted the SMD run with coupled CV based on the remaining peptides’ COM. The work required to liberate all five peptides is reduced at low pH compared to neutral pH (fig. S5). The largest amount of work is required to release the first peptide, and then the amount of work subsequently decreases as the following peptides are liberated. This result indicates that the release of the first peptide is critical and possibly the rate-limiting step in activating peptides to perform membrane disruption, which motivates further investigation into the mechanism and energetics of the release of a single peptide.

To further explore the γ liberation pathway, we performed 20 independent SMD simulations at each pH, using single-peptide CV. The temporal evolution of external work done by the pulling force for the individual trajectories is shown in Fig. 5A, and a clear separation (nonoverlapping distribution) is observed. The final work values are compared in Fig. 5B. At neutral pH, the average work is 403 (± 8) kcal/mol, and at low pH, the average work is 291 (± 5) kcal/mol, a reduction of over 100 kcal/mol when the pH is lowered. The SMD bias velocity was 0.1 nm/ns in both biasing approaches (coupled- and single-peptide CVs). To check that the work trends are not highly sensitive to the pulling velocity, a series of SMD trials were performed at different pulling rates (0.2, 0.4, and 0.8 nm/ns; see fig. S6), which demonstrate that the results are not unique to the pulling velocity chosen. The average maximum force is significantly higher (1187 ± 21 pN) at neutral pH compared to low pH (948 ± 35 pN), indicating that stronger (or more) interactions resist the γ liberation process in neutral pH conditions (Fig. 5C). The maximum forces are large; however, other viruses, including T = 3 non-enveloped capsids, have been shown to be able to withstand applied forces in the range of 1000 pN without undergoing irreversible deformation. This is based on atomic force microscopy (AFM) nanoindentation experiments and CG MD simulations (4346). The SMD results appear consistent with the LMI analyses (Fig. 4, D and E), which indicated that the system is more strongly correlated and rigid at neutral pH. Hence, our multiple SMD simulations confirm that the energetic requirements to release a γ peptide are significantly reduced under low pH conditions, which is in qualitative agreement with the previous experimental work (20). While the nonequilibrium work values can be related to the equilibrium free energy changes through the Jarzynski equality (47), this method is difficult to converge and can require greater than 100 trials (48, 49). Therefore, we next attempted to estimate the equilibrium free energy profiles of γ liberation by performing extensive umbrella sampling simulations along the minimum work pathways under the different pH conditions.

Fig. 5 SMD results using single-peptide CV.

(A) Work versus Z for the 20 trials at each pH. (B) The final work values are sorted and presented from smallest to largest. The average over the trials is shown at the far right. (C) The maximum force values are presented and ordered from smallest to largest. The average over the trials is shown at the far right. In (B) and (C), the error bars on the average values are the SEM.

pH effects on energetic barriers for γ liberation

To calculate a free energy profile of γ liberation under the different pH conditions, the lowest work trials at each pH (Fig. 5A) were selected as initial pathways. These pathways were then subjected to extensive umbrella sampling simulations to allow the systems to relax from the nonequilibrium pathways to a locally equilibrated pathway. To detect that each window had properly equilibrated and was well sampled, autocorrelation-based time series analysis was performed (50). Umbrella windows were run for up to 322 ns to retain approximately 20% (18.9% was the minimum) of the dataset and have a sufficient (>35) Neff in each window (Fig. 6, A and B). A wide range of equilibration times are observed, but we consistently observe fast equilibration times in the high numbered windows (>60), where the peptide is externalized and dissociated from the capsid. Because the peptide is effectively free in solution at this stage of the pathway, it is reasonable to expect that these simulations would rapidly equilibrate. As an additional check of sampling quality, we measured the helicity of the externalizing γ peptide during the equilibrated phase of the umbrella sampling windows (fig. S7). The helicity of γ has been studied previously as a function of membrane composition, and orthogonal sampling was evaluated using time-lagged independent components analysis (TICA) (51). Here, we find that there is good sampling in this orthogonal degree of freedom, with all neighboring windows overlapping at low pH and only one window (window 0) not overlapping with its neighbor windows at neutral pH. We do observe that, at low pH, γ samples less helical conformations both inside and outside of the capsid. Several studies have shown that γ is disordered in solution (14, 15, 52), and it is believed that the arrangement of peptides and interactions with the capsid stabilize the helical conformation in the capsid. At low pH, we observe reduced contacts between γ and the capsid (fig. S8), which presumably leads to destabilization of the helical conformation and may facilitate γ release.

Fig. 6 Energetic analysis of γ liberation pathways.

Time series analysis of umbrella sampling data at neutral (A) and low pH (B). The top plot shows the breakdown between equilibration and production data for each window. The middle panel shows the ratio of production to total simulation data, with dashed red line drawn at 20%. The bottom panel shows Neff, with dashed red line drawn at Neff = 20. The y axis is terminated at 200 to highlight lower sampled windows, although many windows have Neff > 200. (C) PMF force at neutral and low pH. Shaded region represents the SD obtained from error analysis using 1000 bootstrapped profiles. (D and E) Structures along the γ liberation pathway at neutral (D) and low pH (E) pathways. Snapshots are the final frame from a given umbrella window.

The free energy profiles were computed from the equilibrated portion of data in each window, which totaled 5.0 μs at neutral pH (window minimum sampling time, 25 ns) and 4.2 μs at low pH (window minimum sampling time, 18 ns). The resultant profiles (Fig. 6C) show a significant reduction in the barrier (ΔGbarr) and overall (ΔGext) free energy change when the pH is reduced. At neutral pH, a sharply increasing profile is observed in the first 1 nm, followed by a more rugged and more gradual increase until reaching a barrier of ~36 kcal/mol at ~3.4 nm. This is followed by a slight decrease in the potential of mean force (PMF) reaching a local minimum of ~34 kcal/mol at ~4.0 nm. The PMF then increases, reaching ~40 kcal/mol when the peptide fully dissociates from the capsid around 6.0 nm. At low pH, a more gradual and jagged profile is observed, reaching a barrier of ~18 kcal/mol at 2.9 nm. This barrier is followed by a sharp decrease in the PMF to ~12 kcal/mol, where the system reaches a metastable state at ~3.4 nm. The metastable state is shown in fig. S9 and is characterized by partial exposure of the γ peptide onto the capsid exterior surface. This is followed by a gradual increase in energy to ~28 kcal/mol at the fully dissociated state at ~6.0 nm. Snapshots corresponding to the initial, barrier, metastable, and final states are shown in Fig. 6 (D and E) for both systems. Overall, we observe that ΔGbarr is reduced by half (36 kcal/mol at neutral pH, 18 kcal/mol at low pH) and ΔGext is reduced from 40 to 28 kcal/mol when the system pH is reduced.


Several non-enveloped viruses including rhinovirus, hepatitis A, reovirus, parvovirus, adenovirus, and FHV contain membrane-disrupting components of their capsids, which are triggered to release/externalize under low pH (endosomal) conditions (4, 5). However, the molecular mechanisms that underlie this process remain unclear and are challenging to address using current experimental techniques. The rapid progress in computational power and simulation algorithms now permits MD to investigate these complex biomolecular processes and can provide important insights at atomic resolution into the mechanism of action, where experimental characterization is rather difficult, if not impossible (53). Our aspiration is that the knowledge obtained in studying the FHV model system may impart insights into other systems and infection mechanisms. One system in particular where this knowledge may translate is to NωV, which also contains γ peptides that have pH-controlled membrane disruption activity (21). While the sequence conservation is low between FHV and NωV (~12%), there is strong structural conservation (fig. S10A). When the structures are aligned, we have observed that H334 of FHV is aligned with H541 of NωV. Furthermore, there is motif QFGH from residue 331:334 in FHV, which is closely matched by the motif QFAH from residue 538:541 in NωV (fig. S10B). This motif lies at the subunit interface, and we speculate that H541 could be a critical residue for controlling NωV γ peptide release. However, in NωV, the capsid matures at pH 5 and then is maximally lytic at pH 8. Therefore, we propose that H541 would protonate to stabilize the mature state of the capsid but would then deprotonate to facilitate γ release in an alkaline environment.

In the present study, we have investigated the mechanism of lytic peptide liberation from the FHV capsid and the role of pH in affecting the structural dynamics and energetics of the liberation process through a multiscale atomistic simulation approach including equilibrium simulations, SMD, and extensive umbrella sampling. On the basis of the FHV structure-based model, it is hypothesized that γ peptides located under the fivefold axis are liberated from the capsid to disrupt cellular membranes and permit infection. It is known from in vitro studies that γ peptides disrupt membranes (1315) and do so with maximal efficiency at pH 6 (20). Furthermore, when the γ peptides are not cleaved from the capsid protein, the particles are effectively noninfectious (11), substantiating the critical role these peptides play in viral infection. From our 0.5-μs-long equilibrium simulations, we are able to observe pore opening at the fivefold axis under low pH conditions (Fig. 3 and movie S2). Our analysis demonstrates that at low pH, protonation of histidine residues affects (decrease) allosteric communications across the capsid, which may lead to the gatekeeper loop opening (pore formation) in the capsid (Fig. 4). The role of the two histidines as “pH sensor” residues, which can affect capsid dynamics, is supported by conservation of the histidine positions in closely related nodaviruses (see multiple sequence alignment in fig. S11). Although the formation of a pore in the low pH capsid is an exciting observation, we do not observe spontaneous translocation of the γ peptides, which would presumably require simulation time scales beyond our current computational capabilities.

It is therefore informative to understand how pore formation in the low pH capsid affect the energetics of the γ liberation process. In this study, the γ liberation process is characterized through a series of SMD runs, followed by extensive umbrella sampling along the most favorable pathways. Similar approaches combining SMD with umbrella sampling have been successfully applied to study processes of other large biological complexes such as the chromatophore (54) and membrane transporters (55). The SMD simulation results provide a qualitative but conclusive demonstration that the work necessary for externalizing γ peptides is distinctly lower at pH 6 than at pH 7, which is well corelated with the optimal pH of FHV membrane lytic activity (20). In our umbrella sampling approach, we allow long relaxation times (>100 ns) and select only well-equilibrated data, based on autocorrelation analysis, to reconstruct the free energy profiles. We observe significant reductions in the barriers and overall free energy change required to externalize a peptide at low pH.

While the overall free energy changes we observed, even at low pH, are positive, there are a few factors that could contribute to this result. One is that the initial pathways based on SMD are far from equilibrium. Although the umbrella sampling simulations are quite long compared to previous all-atom studies on capsids (36, 56), there could be residual strains in the system that are not fully relaxed, leading to overestimation of the free energy profiles along the reaction coordinate. Furthermore, events such as localized capsid protein unfolding, could further facilitate the peptide release, but are too rare (or slow) to be observed on the time scales of our simulations. Another more physical consideration is to contemplate the closed system of the capsid within a membrane-bound vesicle/endosome. The process of liberating the peptide would involve not only the free energy cost of externalizing the peptide but also the free energy of partitioning the peptide into the membrane, which is known to be favorable. Experimental measurements indicate the free energy change γ1 partitioning into lipid bilayers is approximately −8 kcal/mol (14), while CG simulations have estimated ΔG of γ1 binding to different bilayers is in the range of −15 to −22.5 kcal/mol when the peptide is in a helical state (52). This computationally computed free energy of binding is comparable to the computed free energy of externalization (28 kcal/mol) at low pH. Furthermore, the peptide is partially exposed in the metastable state (fig. S9), and if the virus was closely approaching the membrane, the peptide might be able to initiate contact with the bilayer before fully dissociating from the capsid. In this scenario, the overall process could be favorable as the externalization free energy cost to reach the metastable state is only 12 kcal/mol at low pH, but is far too high (34 kcal/mol) to be favorable at neutral pH. It remains unclear if the peptides fully dissociate from the FHV capsid or if they remain in contact with the capsid. There would be a functional benefit for γs to remain in contact with the capsid, as this would allow the capsid to be proximal to the location of membrane disruption. In poliovirus, there is evidence to support that membrane-disrupting peptides are involved in umbilical connections between the capsid and the bilayer (57).

In this work, we have explicitly modeled a curvature restrained 15-subunit subsection of the FHV capsid, which represents 1 of 12 of the complete capsid. It is known that large biomolecular assemblies, such as virus capsid, exhibit long-range collective motions, which can often be low-frequency (low-energy) motions that can couple to functional dynamics of the systems (26, 37, 58). It is possible in the context of the complete capsid that collective motions would facilitate the peptide release, which are damped in our reduced system; however, our LMI analysis indicates that the capsid becomes less collective at low pH (Fig. 4). Another component to the complete virus, which is neglected in the current study, is the viral genome. While we cannot discount the possibility that low pH induces conformational changes in the genome, which then influence the pathway and energetics of γ liberation, there is evidence that γ release and genome release for FHV are separate events. This is based on a recent study that used heat to trigger genome release from the FHV capsid (59). During the heating, they observe multiple transitions measured by differential scanning calorimetry (DSC) for the wild-type virus, but only a single transition for maturation defective mutants. The first transition for the wild-type particle corresponds to γ release as detected by anti-γ immunolabeling and is independent from the genome exposure transition.

There is increasing interest in understanding asymmetric dynamical processes in highly symmetric virus structures, as these are functionally relevant processes and now both MD and cryo-EM are capable of characterizing these phenomena (60, 61). The present study provides the first free energy profile for membrane lytic peptide liberation from a non-enveloped virus at different pH values. The results provide strong evidence that lowering the pH induces increased fluctuations at the capsid fivefold vertex and reduces the energy requirement to externalize the lytic peptide. We identify a metastable state at low pH that positions the peptide in an externalized but capsid-associated state. When the membrane binding step is taken into consideration, the overall process of transitioning a peptide from the capsid interior to the membrane may approach a negative ΔG, especially if membrane binding was from the metastable state. Overall, our simulation study results qualitatively agree with the previous experimental observation and elucidate several molecular-level aspects of how low pH can trigger capsid structural rearrangement, facilitating externalization of membrane lytic peptides that can disrupt the endosomal membrane and allow viral infection to proliferate.


System preparation

In this study, we modeled a subsection of the complete FHV capsid, comprising five asymmetric units surrounding a fivefold symmetry axis (Fig. 2A). The initial coordinates used were taken from the 2.7-Å resolution x-ray crystal structure [Protein Data Bank (PDB) ID: 4FTB]. The N-terminal tail of the A, B, and C subunits and the C-terminal regions of the γ peptides are unstructured and not resolved in the crystal structure. The unstructured regions of the capsid—residues 1 to 55 of chain A, residues 1 to 52 of chain B, residues 1 to 54 of chain C, and residues 385 to 407 of γ peptides—were not modeled into this structure to avoid artifacts from generating initial conformations of large undetermined regions. Four Ca2+ ions are bound to each iASU, and these ions were included in our starting structure. To maintain the capsid curvature and mimic a full-capsid model, weak positional restraints (force constant of 500 kJ mol−1 nm−2) were applied to several Cα atoms, positioned along the edge of the pentameric structure. Restrained atoms are listed and shown in fig. S12A. The restraints are able to maintain the capsid structure (fig. S1), whereas a control simulation performed with no boundary restraints displayed substantial conformation changes (RMSD > 7 Å) in a 350-ns simulation (fig. S12B). The low pH system is attempting to model pH ~6, which is the environment where FHV capsids have maximal lytic activity (20). The low pH model is constructed where all D75, H215, and H334 are protonated at t = 0 in all subunits and all asymmetric units. An additional control simulation was conducted in which the residues were incrementally protonated in nine stages over 40 ns. Further details of the incremental protonation simulation are described in Supplementary Methods, and analysis of that simulation is presented in fig. S13. There are only two histidine residues in each subunit, and the choice to protonate the histidines is based on assuming a pKa (where Ka is the acid dissociation constant) of ~6.5 (62). D75 was protonated based on a previous pKa prediction of ~6.0, and the proposed reaction mechanism for γ cleavage involves a protonated D75 acting as a hydrogen donor (63). pKa calculations were performed using constant pH-MD (64) and H++ server (65), as described in Supplementary Methods. Each system was solvated with the TIP3P water model (66) and was neutralized by the addition of Na+ ions, and additional Na+ and Cl ions were added to generate a physiological ion concentration of 0.15 M NaCl. The total system sizes were approximately 390,000 atoms.

It is worthwhile to discuss the use of fixed protonation states in modeling different pH environments as we have done in this work. Sophisticated methods have been developed to model constant pH environments in atomistic MD simulations, where residue titration states are coupled to the system dynamics (64, 67). This approach has been applied to the study of virus capsid systems (36, 68); however, those studies were performed on a single asymmetric unit using rotational symmetry boundary conditions and in an implicit solvent environment. The constant pH MD method has been recently extended to be compatible with explicit solvent either through a hybrid implicit/explicit solvent approach (69) or through the multisite λ dynamics formalism (70). While the explicit solvent constant pH approaches are very promising, especially for membrane-protein systems, they still come with a high computational demand. To our knowledge, the application of these methods has been limited to system sizes in the range of 50,000 atoms and with individual trajectory lengths no longer than 62 ns (7174). In our study, we are able to achieve simulation times up to 500 ns and total sampling for over 23 μs for an explicitly solvated system of nearly 400,000 atoms, which enables us to uncover principals that underlie acid-induced externalization of γ peptides from the capsid and qualitatively agree with previous experimental observations.

General simulation protocol

All simulations in this study were performed with the GROMACS 2018 simulation package (75) using the CHARMM36m all-atom force field (76). All simulation types, lengths, and trials are listed in Table 1. The systems were first energy-minimized using the steepest descent algorithm followed by 100 ps of constant volume and temperature (NVT) (300 K) followed by 100 ps of constant pressure and temperature (NPT) (300 K, 1 atm) equilibrations. In the NVT and NPT equilibrations, the nonhydrogen (heavy) atoms of the protein were restrained with a force constant of 1000 kJ mol−1 nm−2. Production MD simulations were performed for 500 ns at a temperature of 300 K and a pressure of 1 atm. Constant temperature was maintained using the Nose-Hoover thermostat (77, 78) with a coupling time constant of 1 ps, and constant pressure was maintained using the Parrinello-Rahman barostat (79) with a coupling constant of 5 ps. Simulation timestep was 2 fs, and waters were kept rigid using the SETTLE algorithm (80). Nonwater bonds involving hydrogen atoms were held fixed using the LINCS algorithm (81). Electrostatic interactions were calculated by using the particle mesh Ewald method (82) with a real-space cutoff of 12 Å and a Fourier grid spacing of 1.6 Å. The cutoff for van der Waal interactions was set to 12 Å, with smoothing starting at 10.5 Å. Additional method details related to SMD, umbrella sampling, trajectory analysis, pKa predictions, and incremental protonation simulation are presented in the Supplementary Materials.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We would like to thank K. Boyd for providing guidance and scripts for the time series analysis. We would also like to acknowledge A. Brice, M. Ward, and E. Moharreri for efforts working on this project in earlier stages. We also appreciate discussions with M. Banerjee about the integrity of the capsid following peptide release. Funding: This research was supported by the NIH through grant number R35-GM119762 (to E.R.M.). This work used the Extreme Science and Engineering Discovery Environment (XSEDE) for HPC resources at the Texas Advanced Computing Center (TACC) through allocation TG-MCB140016. Additional computational resources for this work were provided through the University of Connecticut Storrs HPC center. Author contributions: E.R.M. designed the research plan. A.K.J. ran MD simulations. E.R.M. and A.K.J. analyzed simulation data. E.R.M. and A.K.J. wrote and revised the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article