Introduction

Human immunodeficiency virus type 1 (HIV-1) is the most prevalent strain of HIV responsible for the development of acquired immunodeficiency syndrome (AIDS). The HIV-1 envelop (Env) consists of a host cell-derived lipid membrane and viral glycoproteins that play a crucial role in mediating viral entry into host cells. The Env glycoprotein is initially synthesized in the endoplasmic reticulum (ER) as a precursor gp160 and cleaved by furin into two subunits, gp120 and gp41. The non-covalently associated gp120–gp41 complex is transported to the cell surface in the form of a trimer, where it is subsequently incorporated into the envelope of nascent virions during viral assembly. The exposure of Env protein is essential for binding to the primary receptor CD4 and the co-receptors CCR5 or CXCR4, triggering membrane fusion and viral entry. However, this exposure also renders the virus susceptible to immune attack. In response to host immune pressure, Env is densely coated with N-linked glycans added during post-translational modification in the ER and Golgi apparatus, which effectively shield vulnerable epitopes from immune recognition.

Since HIV-1 was identified as the cause of AIDS in the early 1980s, extensive research has been conducted to elucidate the mechanisms of viral infection and immune response. A major advance was the development of soluble gp140 trimers, composing gp120 and the ectodomain portion of gp41, designed to stabilize the prefusion Env trimer for structural and immunological characterization. Native-like Env trimers, such as SOSIP (Sanders et al., 2013), native flexibly linked (NFL) (Sharma et al., 2015), and uncleaved prefusion-optimized (UFO) constructs (Kong et al., 2016), mimic the structure of viral spike and serve as valuable antigen targets for developing small molecule inhibitors and broadly neutralizing antibodies (bNAbs). With advances in structure determination techniques, such as X-ray crystallography and cryo-electron microscopy (cryo-EM), numerous high-resolution structures of gp140, both unliganded and antibody-bound, are now available in the Protein Data Bank (PDB). Even with such progresses, however, the remaining portions of gp41, including the membrane-proximal external region (MPER), transmembrane domain (TMD), and cytoplasmic tail (CT), are relatively understudied. The MPER is a highly conserved region targeted by several bNAbs, including 10E8, 2F5, 4E10, and Z13e1 (Ofek et al., 2004; Cardoso et al., 2005; Pejchal et al., 2009; Huang et al., 2012). However, in most studies, MPER is examined as a monomeric peptide bound to antibodies or as part of the Env trimer in membrane mimetics such as bicelles and nanodiscs (Rantalainen et al., 2020; Yang et al., 2022; Qi et al., 2025), which do not capture the structure and dynamics of unliganded MPER in the context of the gp120–gp41 trimer embedded in a membrane bilayer. Despite the critical roles of the MPER and TMD in mediating the fusion of viral and host cell membranes (Salzwedel et al., 1999; Miyauchi et al., 2005), they are often excluded from structural studies due to difficulties in crystallization of hydrophobic TM segments while maintaining their native conformations. Therefore, nuclear magnetic resonance (NMR) spectroscopy remains one of the few viable methods for obtaining structural information of MPER and TMD. However, the NMR studies from different research groups have yielded conflicting conclusions regarding the oligomeric state of the TMD. Reported TMD structures include trimeric coiled coils (Dev et al., 2016; Kwon et al., 2018; Piai et al., 2020; Piai et al., 2021), monomeric helices (Apellaniz et al., 2015), and trimers that are not tightly bundled (Reardon et al., 2014). Similarly, there are different conclusions about the conformation and orientation of MPER. Some studies suggest that gp41 bends at the MPER-TMD boundary (around residue 673), with the entire MPER adopting a helical conformation that lies nearly parallel to the membrane (Sun et al., 2008; Kwon et al., 2018), while others show that the C-terminal residues of MPER forms a continuous helix with the TMD (Dev et al., 2016; Chiliveri et al., 2018; Piai et al., 2020; Piai et al., 2021). In addition, the exceptionally long CT plays an important role in facilitating the incorporation of the Env glycoprotein into virions (Checkley et al., 2011). While it is known to contain three conserved amphipathic α-helical segments, referred to as lentiviral lytic peptides (LLPs), proposed models differ in the arrangement of LLPs (Murphy et al., 2017; Piai et al., 2021), and the complete structure of CT and its location in the membrane remains inconclusive.

Molecular dynamics (MD) simulations have been employed to investigate the stability and conformational properties of monomeric and trimeric helical TMD in both aqueous and lipid bilayer environments since late 2000s (Kim et al., 2009; Gangupomu et al., 2010; Baker et al., 2014; Baker et al., 2014; Hollingsworth et al., 2018). Early studies were constrained by limited computational resources and therefore the simulation times are relatively short. Subsequent work employed metadynamics to probe rare events (Gangupomu et al., 2010; Baker et al., 2014), and simulations performed on Anton supercomputers extended sampling to multi-microsecond time scale (Baker et al., 2014). Piai and coworkers determined the NMR structure of a construct comprising the MPER, TMD, and CT, and carried out MD simulations to access the structural stability of the trimeric MPER–TMD–CT complex (Piai et al., 2021). Majumder et al. subsequently simulated the same MPER–TMD–CT complex and applied a machine learning-based approach to classify its conformational ensemble (Majumder et al., 2025). Maillie et al. combined conventional MD, steered MD, and coarse-grained simulations to examine interactions between MPER-targeting antibodies and membrane lipids (Maillie et al., 2025). In addition, MD simulations have been extensively applied to the well-studied ectodomain. Despite these advances, it remains challenging to investigate the gp120–gp41 trimer as an intact entity considering its structural complexity. In this work, we built a model of full-length gp120–gp41 trimer embedded in a lipid bilayer mimicking the lipid composition of the mammalian plasma membrane (Pogozheva et al., 2022) (Figure 1). However, structural information for the CT remains limited, leading to uncertainty in its conformational organization. To reduce potential bias arising from this uncertainty, we also generated a CT-truncated model and used it as the primary system for analysis. We prepared simulation systems varying in the presence of the cleavage site and CT, as well as the initial position of protein in the membrane. Multiple microsecond-long all-atom MD simulations were performed for each system to explore the motions of individual protein domains and the membrane, and to examine how their conformational variability is affected by the difference in the initial configurations. Our results show that the ectodomain undergoes substantial tilting relative to the membrane plane while maintaining a rigid internal structure. In contrast, the MPER and TMD display highly diverse conformations, and their structural variations are influenced by the presence of the CT and the initial TMD position in the membrane. Moreover, we selected several bNAbs targeting the epitopes across different regions of the Env protein and demonstrate that the simulation trajectories can be used to assess the epitope accessibility.

Figure 1.

Results

The ectodomain maintains a rigid internal structure and tilts independently of the TMD

As described in the Methods section, multiple simulation systems were constructed, varying in cleavage, TMD positioning, and CT truncation. Simulations trajectories are denoted as CH^CT1, UL^CT2, CL^ΔCT3, UH^ΔCT3, etc., where the first letter (C/U) indicates cleaved or uncleaved, the second letter (H/L) indicates the high or low initial TMD position, CT/ΔCT indicates the presence or absence of the CT, and the numeric suffix specifies the trajectory index among three replicas (Supplementary file 1–Supplementary Table 1). In all simulations, both the ectodomain and TMD adopted variable orientations relative to the bilayer plane. To quantify the tilt of two domains, we defined the tilt angles for the ectodomain (θ_EC) and TMD (θ_TM) (Figure 2A). θ_EC is measured between the bilayer normal and the vector from the center of mass (COM) of G594 (on all three protomers) to the COM of D167, while θ_TM is measured between the bilayer normal and the vector from the COM of I684 to the COM of V705. Across trajectories, θ_EC typically ranges from 0° to 40°, with only 0.7‰ exceeding 40°. Representative structures for different θ_EC values are shown in Figure 2A. In contrast to the wide range of θ_EC, θ_TM generally remained within 0° to 20° with rare excursions to 30° observed in a small fraction of trajectories. Each 1-µs trajectory is divided into four consecutive 0.25-µs intervals, and data points from each interval are distinguished by four different colors (Figure 2—figure supplements 1–4). The variations of θ_EC and θ_TM over time show that large conformational changes predominantly occurred during the first 0.5 µs, followed by convergence of the θ_EC and θ_TM distributions during the second 0.5 µs in most trajectories. The temporal evolution of θ_EC and θ_TM is additionally shown in Figure 2—figure supplements 5–8. For the CT-truncated systems, the joint probability densities of θ_EC and θ_TM calculated from the final 0.5 µs of each trajectory are shown in Figure 2B, while those for the full-length systems are shown in Figure 2—figure supplement 9. Although the combined dataset spans a broad range of θ_EC (0°–40°) and θ_TM (0°–30°), each individual trajectory explores only a portion of this space. For instance, CH^ΔCT1 samples large θ_EC with small θ_TM, whereas CH^ΔCT3 samples small θ_EC with medium-to-large θ_TM. Thus, multiple independent simulations are essential to capture a wide region of the conformational space. Most importantly, no consistent correlation was observed between θ_EC and θ_TM, either within single trajectories or across all trajectories combined. Pearson correlation coefficients of θ_EC and θ_TM in single trajectories varies between -0.5 and 0.5, with examples such as CH^ΔCT1 and CH^ΔCT2 showing similar θ_TM but distinct θ_EC values, and CL^ΔCT1 and CL^ΔCT3 showing similar θ_EC but distinct θ_TM values. We also calculated the dynamical cross-correlation maps (Ichiye et al., 1991) of Cα atoms for all systems using CPPTRAJ (Roe et al., 2013). The results indicate only very weak correlations between the ectodomain and the TMD (Figure 2—figure supplements 10–13).

Figure 2.

Despite the considerable tilting of the ectodomain relative to the membrane, its internal structure remains rigid and well-preserved throughout the simulations. The root-mean-square fluctuation (RMSF) and root-mean-square deviation (RMSD) of the ectodomain were calculated after the snapshots from each trajectory were aligned to the initial structure by maximizing the overlap in the ectodomain. The majority of the ectodomain displays low RMSF (< 2 Å), and the RMSD stabilized around 4 Å after an initial rise during the first 250 ns (Figure 3A, Figure 3—figure supplement 1). Higher RMSF values were observed in the residues missing from the cryo-EM structure (PDB ID: 6B0N) (Sarkar et al., 2018), which was used for the ectodomain in model building (these missing residues are highlighted in red in Figure 1A, B), and in part of HR1 (Q551-H564) that forms a flexible loop at the interface between two neighboring protomers. In addition, the RMSF of the MPER was also calculated with the trajectories aligned by the ectodomain. The entire MPER (K665–R683) and the adjacent HR2 segment (L660–D664) at the C-terminus of the ectodomain exhibited elevated fluctuations. The gp120–gp41 model was built based on the NMR structure in which the MPER adopts a bent conformation consisting of two helices joined by a sharp turn. The resulting trimer widens from the HR2 helix to the midpoint of the MPER (F673) and narrows from F673 to the TMD. However, this specific conformation was not maintained throughout the simulations. We measured the inter-chain distances between the Cα atoms of corresponding residues (G644, E654, D664, and F673) on neighboring protomers to characterize the structural variation along the HR2 helix and the MPER (Figure 3B). The inter-chain distances of G644 and E654 maintained narrow distributions centered on their initial values, while those of D664 and F673 exhibited broader distributions, reflecting inward shifts of three protomers. This effect was more pronounced in uncleaved systems. In cleaved systems, the HR2 helix interacted with residues M530-N543 and L619-N625 of neighboring protomers, which helps to stabilize the bent MPER conformation (Figure 3C, D). However, these interactions were not consistently observed across all protomers and all simulations, and therefore we still observed considerable conformational variability in this region. In uncleaved systems, the closed loop at the cleavage site occupied the space between the HR2 helix and the neighboring protomer, disrupting their interactions and facilitating inward shifts of the HR2 helix (Figure 3E, F).

Figure 3.

The energetically unfavorable R696 in the hydrophobic core results in asymmetric, kinked TMD conformations and disrupts membrane integrity

Unlike the predominantly hydrophobic TMDs commonly found in many viral envelop proteins, the gp41 TMD contains multiple charged residues: R683 at the N-terminal boundary, R707 and R709 at the C-terminal boundary, and a central arginine, R696 (Figure 4A). In simulations, R683 consistently interacted with lipid headgroups in the exoplasmic leaflet, while R707 and R709 interacted with those in the cytoplasmic leaflet, together acting as anchors that secure the TMD within the bilayer. In the NMR structure (PDB ID: 7LOI) (Piai et al., 2021), the side chain of R696 forms cation-π interaction with the side chain of F699 and a hydrogen bond with the backbone carbonyl group of L692, but it is oriented outward from the helical bundle. When embedded in the bilayer, however, this configuration is energetically unfavorable, as the positively charged side chain directly contacts hydrophobic lipid tails (Cheng et al., 2012). Early in the equilibration stage, the TMD rapidly rearranged to allow R696 residues to interact with more favorable partners, including negatively charged lipid headgroups from either leaflet, ions and water molecules diffusing into the bilayer center, as well as polar and positively charged groups in the CT when present. Once the interactions between R696 residues and their binding partners (lipid headgroup, ions or CT residues) were established, they remained stable with minimal changes throughout the production stage. Because the limited space at the TMD core can only accommodate at most two inward-facing arginine residues, at least one R696 is forced outward to interact with lipid headgroups or CT residues. The differences in arginine orientation and interacting partners give rise to asymmetric protomer conformations and distinct TMD tilts (Figure 4A–F). When R696 points outward, its interactions with lipid head groups or CT residues can be strong enough to destabilize the local helix, introducing a kink into the TMD. Representative snapshots from different trajectories illustrate these asymmetric kinked conformations (Figure 4—figure supplements 1–8). Beyond local deformation of the TMD, R696–lipid interactions perturb bilayer organization, inducing the translocation of lipid headgroups and water molecules toward the bilayer center (Figure 4G). We calculated the interaction frequencies of each TMD residue with lipid headgroups, lipid tails, and water molecules. The results show that the membrane disruption is persistent since many non-terminal TMD residues maintain frequent contacts with water and lipid headgroups throughout the simulation (Figure 4H). Previously, Kim et al. reported that the inter-chain interactions between protonated R696 gradually diminished over a short simulation time (23 ns), leading to increased crossing angles and reduced bundle length (Kim et al., 2009). Gangupomu et. al and Baker et. al observed that R696 snorkeled toward either exoplasmic or endoplasmic headgroups in simulations of the TMD monomer, resulting in TMD tilting and membrane thinning due to water penetration and lipid headgroups interacting with R696 (Gangupomu et al., 2010; Baker et al., 2014; Baker et al., 2014). These observations are consistent with our finding. Hollingsworth et. al also reported membrane thinning; however, they attributed this effect to interfacial interactions of R683 and R707 with both leaflets and proposed that R696 only interacted with water and ions permeating into the center of the TMD timer (Hollingsworth et al., 2018).

Figure 4.

To explore whether R696 exhibits a preference for the exoplasmic versus cytoplasmic leaflet, or interacts with either randomly, we generated two initial structures (high and low) with the TMD positioned at two distinct positions, separated by 4 Å along the membrane normal (see Methods). In simulations initiated from the “high” TMD configuration, R696 residues in three protomers interacted with the lipid headgroups in either leaflet (Figure 4B, C, Figure 4—figure supplements 1, 3) In contrast, in simulations initiated from the “low” TMD configuration, R696 residues interacted exclusively with the headgroups in the cytoplasmic leaflet (Figure 4E and F, Figure 4—figure supplements 2, 4). In the full-length systems, the plate-shaped CT occupies substantial space in the cytoplasmic leaflet, displacing lipids during model construction. Because the CT is not thick enough to fully span the cytoplasmic leaflet, an empty gap remained between the CT and the exoplasmic leaflet (Figure 1D). Over time, lipids in the exoplasmic leaflet shifted downward while the CT residues moved upward to fill this space, causing local bilayer thinning. Under these conditions, upward-oriented R696 can still interact with the headgroups in the exoplasmic leaflet, but downward-oriented R696 primarily contacted CT residues, with rare cases of lipids migrating upward from the cytoplasmic leaflet and approaching the protein. A comparison of all full-length systems (Figure 4—figure supplements 5–8) shows that R696 preferentially adopted downward orientations in the simulations initiated from the “low” TMD configuration (Figure 4—figure supplements 6, 8). In these cases, upward-shifted CT residues interact with the C-terminal half of the TMD, rather than solely with R696, resulting in deeper burial of the TMD in the membrane.

MPER adopts diverse conformations, and its exposure depends on both MPER and TMD conformations

Starting from the initial helix-turn-helix conformation consisting of two separate helical segments, the N-terminal half (MPER-N) and the C-terminal half (MPER-C), MPER underwent rapid rearrangements, and a wide variety of conformations were sampled across all trajectories from all systems. In the initial structure, the trimeric MPER was positioned perpendicular to the membrane, with the helical MPER-N tilted inward and MPER-C tilted outward (Figure 4A). Such conformation and orientation were maintained in some trajectories such as CL^ΔCT3 (the third trajectory of the cleaved, CT-truncated system with the low TMD position, Figure 4—figure supplement 2C). In other trajectories, such as CL^CT1, the helix-turn-helix MPER in one protomer shifted into a horizontal orientation parallel to the membrane surface (Figure 4—figure supplement 6A). In UL^ΔCT1, the entire MPER adopted a more vertical arrangement, with both MPER-N and MPER-C tilted outward (Figure 4E, Figure 4—figure supplement 4A). We also observed in UH^ΔCT3 and UL^ΔCT3 that the HR2 helix in the ectodomain, MPER, and TMD merged into a continuous long helix (Figure 4C, F, Figure 4—figure supplements 3C, 4C). In addition, loss of helical structure within the MPER was common, particularly in the MPER-C region, which often transitioned to a random coil. Therefore, the distinct MPER conformations reported in various experimental studies can be all valid as each captures a possible state within the highly flexible conformational landscape of the MPER. It is noteworthy that the MPER in three protomers can adopt different conformations and orientations, resulting in asymmetric local structures.

Because the MPER is a target of multiple bNAbs, we next examined how its exposure is affected by TMD conformation. As described above, interactions between R696 and lipid headgroups affect the burial depth of the TMD, thereby influencing the positioning of the adjacent MPER. To quantify MPER exposure, we measured the vertical distance from the Cɑ of F673, approximately the midpoint of MPER, to the highest point of the neighboring lipid headgroups, denoted by d_F673. Positive values indicate that F673 lies above the bilayer surface, and negative values indicate its membrane burial. In the initial ‘low’ and ‘high’ TMD configurations, d_F673 was 6.1 Å and 9.1 Å, respectively, but across simulations it spanned a wide range from -15 Å to 20 Å (Figure 5A, B, Figure 5—figure supplement 1). Two examples illustrate this variability. In the first example, all three R696 residues interacted with the cytoplasmic leaflet, drawing the TMD deeper into the membrane. Consequently, the entire MPER-C and most of MPER-N were buried in the membrane, with one F673 positioned 11.3 Å below the membrane surface, thus bringing the ectodomain in close proximity to the membrane (Figure 5C, E). In the second example, two R696 residues interacted with the exoplasmic leaflet, while one interacted with the cytoplasmic leaflet. In the protomer with upward-oriented R696, the MPER-N, MPER-C, and TMD formed a continuous helix. As a result, most of the MPER extended outside the bilayer, with d_F673 reaching 18.0 Å and the ectodomain displaced farther from the membrane (Figure 5D, F). The d_F673 distributions across all cleaved CT-truncated systems suggest that d_F673 tends to be smaller when the simulations started from the “low” TMD configuration (Figure 5B). In both “high” and “low” configurations, the mean d_F673 calculated from the simulation trajectories is smaller than its value calculated from the initial structure, due to both protein conformational change and lipid diffusion. In the “low” configuration, d_F673 decreased by 5.0 Å (from 6.1 Å to 1.1 Å), a slightly larger reduction than in the “high” configuration where d_F673 dropped by 3.5 Å (from 9.1 Å to 5.6 Å), indicating that deeper MPER burial arises from a larger-scale protein motion rather than merely from a lower initial placement.

Figure 5.

Ectodomain epitopes are conditionally accessible, whereas MPER epitopes are virtually inaccessible in the closed prefusion state

To access the accessibility of epitopes on different parts of the trimeric Env protein in the prefusion state, we quantified how frequent their epitopes were exposed without steric clashes from neighboring protein residues, glycans or membrane lipids to six selected antibodies: PGT128 targeting the V3 loop and the N332 glycan (PDB ID: 5JSA) (Kong et al., 2016); PG9 targeting the V1/V2 loop (PDB ID: 3U2S) (McLellan et al., 2011); VRC01 targeting the CD4 binding site (PDB ID: 4LST) (Zhou et al., 2013); 35O22 targeting the gp120–gp41 interface (PDB ID: 4TVP) (Pancera et al., 2014); and 10E8 (PDB ID: 6VPX) (Rantalainen et al., 2020) and 4E10 (PDB ID: 1TZG) (Cardoso et al., 2005) targeting the MPER (Figure 6A, Supplementary file 1–Supplementary Tables 2–5).

Figure 6.

The epitope of PGT128 consists of the V3 loop and N332 glycan, and it is widely recognized that PGT128 binding is mediated by the protein-protein interaction with V3 loop and the protein-glycan interaction with N332 glycan, facilitated by its extra-long HCDR3 loop penetrating the glycan shield (Figure 6—figure supplement 1). Due to the absence of defined secondary structures, glycans exhibit greater conformational flexibility than the protein, although the crowded surface of gp120–gp41 may partially restrict the motion of glycans. N332 glycan must adopt specific conformations to enable specific interactions with PGT128, whereas most conformations block antibody approach. Additional glycans on N137, N156, and N301 can further occlude the site (Figure 6B). Epitope accessibility to PGT128 varied widely across protomers and trajectories, with many cases exceeding 35%. When considering all trajectories from all systems collectively, approximately half of them exhibited at least one protomer with >35% accessibility (Supplementary file 1–Supplementary Table 2). For instance, in trajectory CH^ΔCT1, two protomers showed accessibility >40%, whereas in CL^CT2, all three protomers are <5%. The second antibody, PG9, targets the V1/V2 apex, where binding can be hindered by six glycans, N156, N160, and N185E on the same protomer, and N160, N185E, and N185H on the neighboring protomer (Figure 6—figure supplement 2). In most trajectories, the epitopes on all three protomers were occluded in >95% of the snapshots, with a few exceptions (CH^ΔCT1, CH^ΔCT2, and CL^ΔCT3) where accessibility was non-negligible (Supplementary file 1–Supplementary Table 3). The third antibody, VRC01, targets the CD4 binding site where six glycans, N185H, N197, N276, N363, and N462 on the same protomer and N301 on the neighboring protomer, can interfere with antibody binding (Figure 6—figure supplement 3). Similar to PGT128, the VRC01 epitope is moderately to highly accessible on at least one protomer in many trajectories, whereas in some trajectories, the epitopes on all three protomers are nearly completely shielded (Supplementary file 1–Supplementary Table 4).

35O22 targets the gp120–gp41 interface that is not intrinsically membrane-proximal. However, its binding orientation combined with ectodomain tilting can lead to clashes of 35O22 with the membrane (Figure 6—figure supplement 4). Comparing accessibility with and without including the steric effects of membrane lipids shows that the membrane can hinder 35O22 binding, particularly when it approaches from the direction of ectodomain tilt (Figure 6C). For example, in system UH^ΔCT2, the accessibility frequencies of the three protomers were 49%, 5%, and 25%, when only the shielding of three glycans (N88, N234, and N618) was considered. When the steric effects of both glycans and membrane lipids were included, the first protomer’s frequency dropped from 49% to 1%, while the other two remain essentially unchanged (Supplementary file 1–Supplementary Table 5).

The last two antibodies, 10E8 and 4E10, target the MPER. In the PDB structures of 10E8, the epitope is a 17-residue peptide corresponding to residues 671-687 in this work. In the PDB structure of 4E10, the epitope is a 12-residue peptide corresponding to residues 669-680. Superposition of these antibodies onto the initial simulation structure revealed extensive clashes with the protein and glycans linked to N611, N618, and N637 on the neighboring protomer, as well as slight steric clashes with the membrane (Figure 6—figure supplements 5, 6). Using the strict criteria (≤10 heavy-atom clashes with protein/glycan and ≤20 with lipids), no snapshot was identified in which the epitope of either antibody was accessible. Even with relaxed criteria (≤20 heavy-atom clashes with protein/glycan or ≤40 with lipids), accessible cases remain rare, which require either large ectodomain tilts to create space on the opposite side or the MPER transitioning into continuous helices to displace the ectodomain from the membrane (Figure 6D, Figure 6—figure supplements 7–10). Given the extremely low frequency of such specific conformations, MPER epitopes are effectively inaccessible in the prefusion trimer, suggesting that MPER-targeting antibodies such as 10E8 and 4E10 act at later stages of viral entry.

Discussion

The hydrophobic and flexible nature of the MPER has made it difficult to crystallize, and most structural insights have come from NMR studies. However, structural information on the MPER and TMD has been inconsistent in the literature. For example, one study of the MPER in DPC micelles suggested that the MPER adopted a distorted helical structure lying parallel to the membrane surface (Sun et al., 2008), while NMR studies using gp41 fragments containing part of the MPER and TMD embedded in bicelles or bilayers reached conflicting conclusions. Some proposed a kink (Apellaniz et al., 2015) or turn (Kwon et al., 2018) between the MPER and TMD, whereas others suggested that the MPER and TMD formed a continuous helix aligned perpendicular to the membrane (Dev et al., 2016; Chiliveri et al., 2018). In the NMR structure (PDB ID: 7LOI) used to build our model structure, a kink appears between the MPER and TMD, along with a turn in the middle of the MPER (Fu et al., 2018; Piai et al., 2021). Crystal structures of MPER peptides bound to various antibodies show that in addition to the helix-turn-helix conformation, MPER can also adopt alternative structures in which the MPER-C remains helical while part of the MPER-N becomes unstructured (Williams et al., 2017). Such variability underscores the intrinsic structural flexibility of the MPER, with each experimental method capturing only a subset of possible states under specific conditions. In this study, we employed MD simulations to sample a broader conformational landscape and gain a more comprehensive view of the structural heterogeneity of the MPER. The results reveal that the MPER can adopt all of the conformations reported experimentally. This is consistent with the all-atom simulations of MPER–TMD–CT and MPER–TMD in an asymmetric membrane conducted by Majumder et al., which likewise show multiple different conformational states of the MPER and TMD (Majumder et al., 2025). While MPER plasticity has been linked to its role in virus-host membrane fusion because it enables the ectodomain and TMD to adopt distinct orientations during large-scale structural rearrangements, our results show that this flexibility is already present in the prefusion state. At this stage, the MPER functions as a flexible hinge that facilitates ectodomain tilting, which can be critical for the spatial alignment of the CD4-binding site with the host receptor for efficient receptor engagement.

Comparison of the full-length and CT-truncated systems shows that the primary difference arises from changes in the lipid bilayer, particularly in the exoplasmic leaflet, whereas differences in protein conformation and dynamics are less evident. Previous experimental studies have reported that mutations of the TMD residue and CT truncation can substantially affect antigenicity of ectodomain (Edwards et al., 2002; Chen et al., 2015; Dev et al., 2016). However, the ectodomain remains relatively rigid in our simulations for both full-length and CT-truncated systems. It is unclear whether this behavior reflects insufficient conformational sampling or artifacts associated with the model structures. Structural information for the CT is very limited, and the NMR structure (PDB ID: 7LOI) was the only available CT structure at the time the simulation systems were constructed. As a result, the extent to which this structure represents the native CT conformation remains uncertain. Additional experimental structural characterization of the CT will be important for achieving a more complete understanding of its functional role. Similar limitations apply to other modeled regions where structural information is incomplete, including missing loops in the ectodomain, the cleavage site and heptad repeat 2 where two PDB structures (IDs: 6B0N and 7LOI) were merged. These regions introduce additional uncertainty, and the extent to which they influence the interpretation of our results remains an open question. Our analysis of the ectodomain and TMD coupling indicates that the motions of these two domains are largely independent. This observation does not contradict experimental studies demonstrating functional coupling between the TMD, CT, and the antigenic profiles of Env (Chen et al., 2015; Dev et al., 2016). Munro et al. proposed that unliganded Env is intrinsically dynamic, transitioning among three distinct prefusion conformations: a closed ground state (predominant), a transient state, and a CD4-/co-receptor-stabilized state. Both laboratory-adapted and clinically isolated strains can spontaneously transition among these three states, although their relative occupancies differ (Munro et al., 2014). It is therefore possible that TMD mutations or CT truncation also alter the equilibrium distribution among three states, thereby affecting the epitope exposure, particularly for epitopes that are occluded in the closed ground state while exposed in the CD4-/co-receptor-stabilized state. However, transition among three states occur on millisecond-to-second timescales. Our simulations on microsecond timescales primarily capture conformational variations within the closed ground state and suggest that the MPER acts as a hinge, providing substantial flexibility that enables the ectodomain and TMD to move independently while Env remains in the closed ground state. Enhanced sampling methods could be applied to more thoroughly explore the conformational landscape, including not only variations within the closed ground state but also transitions among the closed ground, transient and CD4-/co-receptor-stabilized states. In addition to the limitations inherent to protein structure modeling, the choice of lipid composition remains an open question. In this work, we selected an asymmetric mammalian plasma membrane because it is one of the 18 complex biomembrane systems we previously studied (Pogozheva et al., 2022), and among them, it provides the closest available approximation to the HIV membrane. Nevertheless, experimental studies have reported differences in lipid composition between HIV virions and the host plasma membrane (Aloia et al., 1993; Brugger et al., 2006; Huarte et al., 2016; Mucksch et al., 2019; Tomishige et al., 2023). Although we do not anticipate that our main conclusions regarding Env domain motions and MPER flexibility would change substantially, evaluating the influence of lipid composition represents an important direction for future work.

To probe antibody recognition, we assessed epitope accessibility using snapshots extracted from simulation trajectories. Unlike static experimental structures, this approach incorporates protein dynamics as well as the influence of glycans and membrane lipids, yielding a frequency-based measure indicating how often each epitope is accessible to its corresponding antibody, rather than a binary accessible/occluded classification. It is worth noting that viral glycans are highly variable in both site occupancy and the specific glycoform present at each glycosylation site, as indicated by mass spectrometry data (Behrens et al., 2016; Cao et al., 2017). For structure modeling, however, we assumed full occupancy at all glycosylation sites and selected a single representative glycoform per site. As a result, the estimated glycan shield may not fully reflect the in vivo situation. In addition, this analysis is based on geometric and steric criteria without accounting for potential conformational adaptations of gp120–gp41, antibodies, or the membrane; therefore, the calculated frequency of antibody accessibility should be interpreted as an approximation rather than a definitive indicator of binding competence. Our results indicate that epitopes on the ectodomain, though heavily shielded by glycans, can still become transiently accessible in the closed prefusion state. The frequency provides a quantitative measurement of how vulnerable each epitope is to antibody binding. When estimating steric clashes between the membrane and antibodies, we adopted relaxed cutoffs to account for lipid fluidity, but the membrane due to fast lipid dynamics may in reality accommodate even greater steric overlap. Although the MPER can transiently extend out of the membrane when gp120–gp41 remains in the prefusion state, MPER-targeting antibodies still encounter substantial steric hindrance from the bulky gp120 subunit, surrounding glycans, and the membrane throughout nearly all simulations. Therefore, these antibodies are unlikely to adopt the proper orientations required to bind their epitopes. This is consistent with experiment studies indicating that MPER-targeting antibodies bind effectively only after the gp120–gp41 trimer undergoes major conformational rearrangements toward a fusion-intermediate or post-fusion state (Frey et al., 2008; Alam et al., 2009; Chen et al., 2014; Lee et al., 2016). Shehata et al. also built a model of full-length gp120–gp41 trimer embedded in a lipid bilayer and performed all-atom simulations, in which a tilting motion of the ectodomain was observed. Based on the analysis of accessible surface area using different probe radii, they reported that antibody epitopes on the ectodomain are largely shielded by glycans, while the MPER epitope is mainly occluded by the membrane with tilt angles above 30° required to achieve greater MPER exposure (Shehata et al., 2025).

Overall, the data presented here demonstrate that structural modeling integrated with molecular dynamic simulations can be applied to complex biomolecular systems, such as the gp120–gp41 trimer embedded in a membrane bilayer, which enables detailed characterization of protein dynamics at the molecular level and facilitates the investigation of antigen-antibody interactions, thereby offering potential guidance for rational vaccine design.

Materials and methods

System building

Modeling of full-length gp120–gp41 trimer

The full-length gp120–gp41 trimeric model was generated by combining the crystal structure of the ectodomain with the NMR structure containing the MPER, TMD, and CT (Figure 1A, Figure 1—figure supplement 1A). Considering structure resolution and the number of missing residues, we selected the crystal structure of a soluble Env trimer exhibiting a closed-form, native-like prefusion conformation (PDB ID: 6B0N) (Sarkar et al., 2018), which contains the entire ectodomain including the V1-V5 loops, fusion peptide (FP), heptad repeat 1 (HR1), and heptad repeat 2 (HR2) (Figure 1B). In contrast to the ectodomain well represented in the PDB with numerous crystal and cryo-EM structures, structural data for the MPER, TMD and CT are scarce. For these regions, we chose the NMR structure (PDB ID: 7LOI) (Piai et al., 2021) containing the MPER, TMD, and CT. Both 6B0N and 7LOI are trimeric structures and the C-terminus of 6B0N includes five residues (L660–D664) overlapping with the N-terminus of 7LOI. Measurement of the Cα–Cα distances for D664 in the three protomers revealed that the inter-protomer spacing in 7LOI (16 Å) was much shorter than in 6B0N (39 Å) (Supplementary file 1–Supplementary Figure 1A). To reconcile this difference, we performed short MD simulations on 6B0N with residues E32–Q640 fixed and a distance restraint applied to D664 to reduce the Cα–Cα distance to 33 Å. Similarly, for 7LOI, residues I675–L856 were fixed, and a distance restraint was applied to D664 to increase the Cα–Cα distance to 33 Å (Supplementary file 1–Supplementary Figure 1B). The two adjusted structures were then combined by aligning the overlapping residues (Supplementary file 1–Supplementary Figure 1C). The missing loops in 6B0N (T63, D149, E185A–N185I, S401–G409) and the unresolved region between the TM and CT in 7LOI (F717–G738) were grafted from a modeled structure generated by I-TASSER (Yang et al., 2015). 6B0N is a native flexibly linked (NFL) trimer in which the furin cleavage site ⁵⁰⁸REKR⁵¹¹ is substituted with a 10-residue linker ⁵⁰⁸GGGGSGGGGS⁵¹¹, but the underlined residues are unresolved in the structure. To construct the cleaved gp120–gp41 model, we removed the initial “GG” and final “S” from 6B0N. For the uncleaved gp160 model, we closed the cleaved loop by adding ⁵⁰⁸REKR⁵¹¹ back to the structure. If we simply connect G507 and A512 without adjusting the conformation of flanking loops on both sides of the cleavage site, two neighboring protomers become entangled, yielding a knot-like fold (Supplementary file 1–Supplementary Figure 2A, B). In 6B0N, the helical HR2 region lies at the C-terminus of each protomer, allowing it to thread through a loop formed by the neighboring protomer (Supplementary file 1–Supplementary Figure 2C). However, in the context of the full-length gp120–gp41, it is implausible that the three protomers are entangled as it does not conform to a realistic folding pathway. Therefore, we performed short MD simulations with restraints to adjust the conformations of the flanking loops and the HR2 helix before closing the loop with REKR fragment (Supplementary file 1–Supplementary Figure 2D, E). Two mutations (S764C and S837C) were introduced in the CT to restore the palmitoylation sites, and lipid tails oriented towards the hydrophobic core of the bilayer were then attached to the palmitoylation sites using the PDB Manipulation module in CHARMM-GUI (Jo et al., 2008; Jo et al., 2014; Park et al., 2023) (Figure 1D). In addition to the full-length model, we built the CT-truncated model by removing the modeled residues F717–G738 that are originally unresolved in the PDB structure 7LOI, and the residues E739–L856 that form the large CT plate.

Glycosylation

N-linked glycans were modeled using Glycan Reader & Modeler (Jo et al., 2011; Park et al., 2017; Park et al., 2019) in CHARMM-GUI. The reported mass spectrometry (MS) experiments have revealed the probability of different types of N-linked glycans at each glycosylation site (Behrens et al., 2016; Cao et al., 2017). For each site, we selected the one with the highest probability. As the MS data provide only the glycan type (high-mannose, hybrid, and complex) and composition (number of each monosaccharide unit), a representative isomer was chosen when multiple isomers corresponded to the same composition. For example, HexNAc(2)Hex(9), i.e., the Man₉ glycan, has a single isomer, whereas HexNAc(2)Hex(8), i.e., the Man₈ glycan, has multiple isomers depending on which non-reducing terminal mannose is trimmed from Man₉. For complex N-linked glycans, one composition can correspond to multiple isomers differing in the number of branches on the α1-3 and α1-6 arms and in the placement of non-reducing terminal neuraminic acid (Neu5Ac). In the absence of linkage-specific information, one isomer was arbitrarily selected in cases of multiple possibilities. The selected glycan sequences for the 27 glycosylation sites in each promoter are summarized in Supplementary file 1–Supplementary Table 6. The select glycan sequences were represented in the Glycan Reader Sequence format (Jo et al., 2011; Park et al., 2017) and added to the corresponding glycosylation sites using the Glycan Reader & Modeler graphical interface.

Membrane

The full-length and CT-truncated models were embedded into an asymmetric lipid bilayer with the lipid composition corresponding to a mammalian plasma membrane (van Meer et al., 2008; Sampaio et al., 2011; Ingolfsson et al., 2014; Pogozheva et al., 2022), which consists of phosphatidylcholine (PC), phosphatidylethanolamine (PE), phosphatidylinositol (PI), phosphatidylserine (PS), phosphatidic acid (PA), sphingomyelin (SM), cholesterol (CHOL), and glucosylceramide (GlcCer). In the CT-truncated system, the exoplasmic and cytoplasmic leaflets contain similar numbers of lipids, despite differences in composition. In the full-length system, the cytoplasmic leaflet contains approximately 100 fewer lipids than the exoplasmic leaflet due to the space occupied by the CT. The lipid composition is summarized in Supplementary file 1–Supplementary Table 7. Membrane systems were constructed using CHARMM-GUI Membrane Builder, which provides a user-friendly graphical interface for selecting lipid types and defining their numbers in each leaflet (Jo et al., 2007; Jo et al., 2009; Wu et al., 2014; Lee et al., 2016; Lee et al., 2019). Although it is recognized that the charged residues R683, R707, and R709, flanking the N- and C-termini of the TMD, are anchored within the lipid headgroups, there remains some flexibility to shift the TMD slightly up or down within the membrane. Assuming the lipid bilayer is aligned parallel to the xy-plane, the protein was positioned at two distinct depths along the z-axis, differing by approximately 4 Å, referred to as the “high” and “low” TMD configurations. To examine whether this subtle difference in the initial configurations leads to distinct conformational changes during simulations, we initiated independent simulations from each configuration. The glycoprotein-membrane system was solvated in a box of approximately 210 × 210 × 260 ų using the TIP3P water model (Jorgensen et al., 1983), and KCl was added at a concentration of 0.15 M to neutralize the system. The final simulation input files were automatically generated by CHARMM-GUI.

Simulation details

The combination of cleaved vs. uncleaved, full-length vs. CT-truncated, and high vs. low TMD positions in the membrane results in eight distinct configurations, and we performed three independent 1-μs all-atom MD simulations for each configuration. The CHARMM36(m) force field (Guvench et al., 2009; Klauda et al., 2010; Raman et al., 2010; Guvench et al., 2011; Huang et al., 2017) was used for proteins, carbohydrates, and lipids. The total number of atoms is approximately 1,100,000 (including ∼300,000 water molecules and ∼1,700 ions), with slight variations between systems. van der Waals interactions were smoothly switched off over 10–12 Å using a force-based switching function (Steinbach et al., 2004), and long-range electrostatic interactions were calculated using the particle-mesh Ewald method (Essmann et al., 1995) with a mesh size of ∼1 Å. All simulations were performed using GROMACS (Van Der Spoel et al., 2005). Bond lengths and angles involving hydrogens were constrained using the LINCS algorithm (Hess et al., 1997). The system was equilibrated in the canonical (NVT) ensemble at 310.15 K for 2 × 1.25 ns with a 1-fs time step using the Berendsen thermostat (Berendsen et al., 1984) (coupling constant 𝜏_t = 1 ps), and then switched to the isothermal–isobaric (NPT) ensemble for 1.25 ns with a 1-fs time step, followed by 3 × 5 ns with a 2-fs time step using the Berendsen thermostat (𝜏_t = 1 ps) and barostat (reference pressure = 1 bar, coupling constant 𝜏_p = 5 ps, compressibility = 4.5×10⁻⁵ bar⁻¹). Positional and dihedral restraints were applied to proteins, glycans, and lipids, with force constants progressively reduced over successive intervals (Supplementary file 1–Supplementary Table 8). In the production run, temperature was maintained using the Nosé–Hoover thermostat (Hoover, 1985; Nosé, 2006) and the pressure coupling was applied using the semi-isotropic Parrinello–Rahman barostat (Parrinello et al., 1981; Nosé et al., 2006). A 4-fs time step was used with the hydrogen mass repartitioning technique (Hopkins et al., 2015) and all restraint potentials were removed. The Python library MDTraj (McGibbon et al., 2015) was used to analyze the simulation trajectories.

Assessment of antibody epitope accessibility

We selected six antibodies based on epitope location and the availability of PDB structures containing both antibody and epitope: PGT128 targeting the V3 loop and the N332 glycan (PDB ID: 5JSA); PG9 targeting the V1/V2 loop (PDB ID: 3U2S); VRC01 targeting the CD4 binding site (PDB ID: 4LST); 35O22 targeting the gp120–gp41 interface (PDB ID: 4TVP); 10E8 (PDB ID: 6VPX) and 4E10 (PDB ID: 1TZG) targeting the MPER. For each antibody–epitope complex, TM-align (Zhang et al., 2005) was used to identify the optimal alignment between the epitope from the PDB structure and each protomer of the trimeric protein, and the resulting rotation matrix was applied to place the antibody relative to each protomer. Steric clashes were then assessed by counting any heavy atom of protein, glycan, or lipid within 2 Å of a heavy atom of the antibody. An epitope was considered occluded by protein and glycans if more than 10 heavy-atom clashes occurred with these components. To account for membrane flexibility, we adopted a more permissive cutoff for lipids: the epitope was classified as occluded by the membrane if more than 20 lipid heavy-atom clashes were detected. To calculate the frequency of epitope accessibility, snapshots were extracted from simulation trajectories, and structural alignment and assessment of steric clashes were performed on each snapshot.

Figure supplements

Figure 2—figure supplement 1.

Figure 2—figure supplement 2.

Figure 2—figure supplement 3.

Figure 2—figure supplement 4.

Figure 2—figure supplement 5.

Figure 2—figure supplement 6.

Figure 2—figure supplement 7.

Figure 2—figure supplement 8.

Figure 2—figure supplement 9.

Figure 2—figure supplement 10.

Figure 2—figure supplement 11.

Figure 2—figure supplement 12.

Figure 2—figure supplement 13.

Figure 3—figure supplement 1.

Figure 4—figure supplement 1.

Figure 4—figure supplement 2.

Figure 4—figure supplement 3.

Figure 4—figure supplement 4.

Figure 4—figure supplement 5.

Figure 4—figure supplement 6.

Figure 4—figure supplement 7.

Figure 4—figure supplement 8.

Figure 5—figure supplement 1.

Figure 6—figure supplement 1.

Figure 6—figure supplement 2.

Figure 6—figure supplement 3.

Figure 6—figure supplement 4.

Figure 6—figure supplement 5.

Figure 6—figure supplement 6.

Figure 6—figure supplement 7.

Figure 6—figure supplement 8.

Figure 6—figure supplement 9.

Figure 6—figure supplement 10.

Data availability

The modeled structures, force field parameters, and GROMACS input files are available on Zenodo [https://doi.org/10.5281/zenodo.17259909]. The code used in this article can be found at https://doi.org/10.5281/zenodo.17259951. The simulation trajectories can be found at https://doi.org/10.5281/zenodo.18853902, https://doi.org/10.5281/zenodo.18854615, and https://doi.org/10.5281/zenodo.18854639.

Additional files

Supplementary Figures and Tables

Additional information

Funding

HHS | National Institutes of Health (NIH) (R35-GM153458)

• Wonpil Im