A broad spectrum anti-HIV inhibitor significantly disturbs V1/V2 domain rearrangements of HIV-1 gp120 and inhibits virus entry.

Abstract Inhibition of human immunodeficiency virus (HIV) entry into target human cells is considered as a critical strategy for preventing HIV infection. Conformational shifts of the HIV-1 envelope glycoprotein (gp120) facilitates the attachment of the virus to target cells, therefore gp120 remains an attractive target for antiretroviral therapy development. Compound 18A has been recently identified as a broad-spectrum anti-HIV inhibitor. It was proposed that 18A disrupts rearrangements of V1/V2 region in gp120; however, the precise mechanism by which 18A interferes with the inherent motion of V1/V2 domain remains obscure. In this report, we elaborate on the binding mode of compound 18A to the closed conformation of a soluble cleaved gp120 and further examine the dynamic motion of V1/V2 region in both gp120 and the gp120–18A complex via all-atom molecular dynamics simulations. In this work, comparative molecular dynamic analyses revealed that 18A makes contact with Leu179, Ile194, Ile424, Met426 W427, E370 and Met475 in the main hydrophobic cavity of the unliganded gp120 and disrupts the restructuring of V1/V2 domain observed in apo gp120. The unwinding of α1 and slight inversion of β2 in gp120 leads to the shift of VI/V2 domain away from the V3 N-terminal regions and toward the outer domain. Stronger contacts between Trp425 and Trp112 rings may contribute to the reduced flexibility of α1 observed upon 18A binding thereby inhibiting the shifts of the V1/V2 region. Binding of 18A to gp120: (1) decreases the overall flexibility of the protein and (2) inhibits the formation a gp120 conformation that closely ressembles a CD4-bound-like conformation. Information gained from this report not only elaborates on important dynamic features of gp120, but will also assist with the future designs of potent gp120 inhibitors as anti-HIV.


Introduction
The human immunodeficiency virus (HIV-1) is a retrovirus that causes Acquired Immune Deficiency Syndrome and approximately 36 million people have died since the beginning of the epidemic (1,2). Despite the use of highly active anti-retroviral treatments to improve the prognosis of patients, the emergence of drug resistant strains of HIV creates an urgent need for the development of new drugs. HIV virus entry into target cells is the initial step of the infection process; therefore, inhibiting this process is a critical strategy to prevent HIV infection. The binding of the HIV-1 envelope glycoprotein (gp120) to the target cell surface receptor CD4 activates conformational changes to gp120 crucial for HIV-1 entry into cells (3,4). The HIV-1 envelope is a trimer consisting of three gp120 glycoproteins non-covalently bonded to three gp41 glycoproteins. The unliganded gp120 glycoprotein exists in a closed state consisting of two domains: the inner and outer domains ( Figure 1). In the CD4-bound (open state), there is an additional bridging sheet domain. Binding of gp120 to the CD4 induces structural changes on the trimer including (i) the relocation of gp120 V1/V2 and V3 regions, (ii) exposure of the co-receptor-binding site of gp120 and (iii) formation of the heptad repeat 1 coil of gp41 (5,6) leading to viral and target cell membrane fusion. Therefore, gp120 is an attractive target and the identification of virus and target cell fusion inhibitors remains an interesting avenue for the development of antiretroviral therapy.
A newly discovered small molecule, 18A ( Figure 2) was identified by high-throughput screening and reported to inhibit the target cell entry of diverse HIV-1 strains including BMS-806-mutants by targeting gp120 and blocking the function of the envelope (16). Compound 18A was identified as the most potent compound amongst tested compounds sharing an acyl hydrazone pharmacophore. Compound 18A does not disrupt gp120 binding to CD4 and CCR5 and its mechanism of action involves the inhibition of two main structural changes in gp120 induced by CD4 including: (i) rearrangement of the gp120 V1/V2 region and (ii) formation and exposure of the gp41 HR1 region (16). The potential binding site of 18A has been mapped to the b20-b21 and V1/V2 regions after testing for sensitivity of different mutants, furthermore it was suggested that compound 18A binds preferentially to the unliganded gp120 (16). Despite such reports, the binding mode of 18A is not yet known. Kwon et al. previously reported the structures of the unliganded gp120 core that resembled the CD4-bound structures contradicting the traditional view that CD4 induces conformational shifts in gp120 (17). Direct imaging of conformational dynamics within the envelope revealed that unliganded HIV-1 envelope transitions between three conformations whose relative occupancies were remodeled by CD4 and antibody binding, also disputing the previous view of CD4 induced structural changes (18). Several previous molecular dynamics simulation studies have demonstrated that the unliganded gp120 exhibits more conformational fluctuations and less stability than the CD4-bound gp120 (19)(20)(21)(22). Molecular dynamics of gp120 has been examined in computational studies on gp120 cores without the inclusion of V1/V2, V3 and V4, N-and C-termini regions or on gp120 with modeled loops (19)(20)(21)(22). Until recently, the V1/V2 region has been referred to as a loop, but recent crystal structures have unveiled the V1/V2 region to consist of a four-strand anti-parallel beta-sheet fold (7). The V1/V2 domain has been shown to modulate viral tropism, and also shield gp120 against neutralizing antibodies (23). Recently, resolved BG505.SOSIP gp140s structures are the only structures of recombinant soluble trimers that closely resemble the native gp120 structure. The recent elucidation of a 4.7-Å BG505.SOSIP.664 trimer crystal ( Figure 1) has revealed novel features of gp120 antigenicity and architecture (7). Experimental evidence has demonstrated that these SOSIP trimers have conformations similar to gp120 structure and that the V1/V2 regions are well ordered (14). However, information on the structural dynamics V1/V2 domain in the context of the full gp120 structure is still lacking in literature.
The interesting discovery of broad anti-HIV activity and proposed mechanism of action of 18A has prompted us to investigate its binding mode to gp120 as well as the structural dynamics implications of its binding. To this end and in order to provide more insight into the mechanism of disruption of V1/V2 region rearrangements by 18A, we performed comparative molecular dynamic analyses on the unbound as well as 18A-bound protein.
Information gained from this report will not only shed light on the dynamic features of gp120, but will also assist in the design of optimized HIV-1 inhibitors.

Ligand and protein preparation
Compound 18A was constructed and energy-minimized using MMFF94S force-field in Avogadro software (24). The crystal structure, PDB: 4TVP (7), of a cleaved soluble BG505 SOSIP unliganded HIV-1 trimer obtained from the RCSB PDB protein structure database (25) was used as starting coordinates. This structure was selected because it contains additional residues that are missing in the other trimer structures. The gp41 chain, antibodies, water, ions, carbohydrates were deleted from the gp120 monomer using Chimera molecular modeling tool (26).
gp120-18A complex system preparation Based on experimental data, two active site pockets were proposed for drug binding with gp120; these are b20-b21 and V1/V2 regions (16) (Figure 3 and Table 1). Compound 18A was docked into both experimentally identified active sites pockets. Docking parameters are shown in Table 1.
AutoDock Tools graphical interface (27) was used to define the grid box for the selected pockets obtained from submitting 4TVP pdb to the CASTp algorithm (28). Docking calculations using the Lamarckian genetic algorithm (29) were performed with AutoDock Vina. For each docking experiment, 10 runs were performed. The gp120-18A complex with the highest binding affinity and lowest RMSD from  docking (S1 and S2) was submitted to MD simulations. LigPlot (30) was used to construct the ligand interactions. In order to assess the reliability of our docking protocol, the same docking approach/parameters were used for NBD-557-gp120 X-ray crystal structure (PBD accession code: 4DVR) (31). The docked complex was then superimposed against the original complex and RMSD was computed (S3 and S4).

Molecular dynamic simulations
Classical all atom explicit solvation molecular dynamics simulation of two systems, apo gp120 and gp120-18A, was performed using GPU version of PMEMD engine integrated with Amber14 (32). The LEAP program of AMBER14 was used to generate topologies for the systems. Both systems were solvated in a rectangular box TIP3P water molecules with a 10 Å buffer and 10 Cl À ions were added to neutralize the systems resulting to a pH of 7, to correlate with the pH used during the high-throughput screening experiment to obtain 18A (16). A total of 26 843 water molecules were added to the gp120 system to form a rectangular box of 88.6 Â 89.7 Â 133.4 Å 3 , resulting in 87 641 atoms in total. A total of 26 843 water molecules were added to the gp120-18A system to form a rectangular box of 88.6 Â 89.7 Â 133.4 Å 3 , resulting in 87 668 atoms in total. The systems were energy minimized by the steepest descent method with weak constraints for 1000 steps. Systems were subsequently heated up from 0 to 300 K and equilibrated under 1 atm pressure over 50 ps in the NPT and NVT ensembles. Finally, 100 ns production runs were performed on both systems under periodic boundary conditions. The particle mesh Ewald method was used to treat long-range electrostatic interactions at 2 fs integration step. The SHAKE algorithm (33) was used to constrain bond lengths and the Langevin algorithm (34) was applied to couple systems to a 300 K external temperature bath.
Trajectory analysis including the root mean square deviation (RMSD), root mean square fluctuations (RMSF), radius of gyration, solvent accessible surface area, inter atomic distances, dynamical cross-correlation matrices was carried out using PTRAJ and CPPTRAJ (33) modules implemented in AMBER14. The graphical user interface of UCSF Chimera package (26) was used to visualize structures. Plots of data were obtained using the GUI of Microcal Origin data analysis software version 6 (www.originlab.com).

Principle component analysis
Dominant and collective modes of gp120 and gp120-18A were obtained from the MD trajectories by performing principle component analysis (PCA) using the CPPTRAJ module of AMBER14. Trajectories were aligned against the fully minimized structures and PCA was performed for C-a atoms on 5000 snapshots each at 50 ps intervals. The first two principal components were calculated and the covariance matrices were generated. The first two principal components correspond to the first two Eigen vectors of the covariance matrix. The PCA scatter plots were then created using GUI of Microcal Origin data analysis software version 6. Structural diagrams and porcupine plots were generated using the ProDy interface (35) of normal mode wizard of VMD (36). Figure 3. Cartoon illustration of the closed conformation of gp120 with amino acid residues experimentally identified (16) to be linked to high sensitivity (blue) and low sensitivity (magenta) to 18A. Two pockets (1) and (2) containing the residues are highlighted in circles.

Results and discussion
Binding mode of 18A with gp120 Compound 18A was docked into two pockets of the closed conformation of gp120 (PDB code: 4TVP) as explained in the ''Methods'' section. The best docked conformation (highest binding affinity and lowest RMSD) clearly indicates that 18A interacts mostly with residues (with a cut-off distance of 5 Å from 18A) lining the largest hydrophobic cavity, a highly conserved cavity, that forms the CD4 binding pocket in the open conformation of gp120 ( Figure 4, Table 2 and S5). In our docked model, 18A binds at the tip of the hydrophobic cavity (Figure 4). The ligand-residue interactions plot (S3) shows that the benzothiadiazole group is surrounded mostly by hydrophobic residues Leu175, Ile194, Ile423 and Leu369. The hydrazone region of 18A, which is associated with gp120 inhibition, forms a hydrogen bond with Glu370 and Asp368. The phenolic ring interacts at the bottom of the cavity forming a hydrogen bond with Met475 ( Table 2). The interaction of 18A with residues Leu179, Lys421 Ile424, Ile194, Met426 and Met475 indicates that our docked model is consistent with experimental data (16). Asp368 and Met426 have been identified as two hotspots in the CD-4 binding pocket for the effective inhibition of gp120 functioning (8), therefore the strong hydrogen bond interaction of 18A (Table 2) with the two hot spot residues validates its superior inhibitory properties.

Stability of systems during MD simulations
The gp120 crystal structure in its closed unliganded form and gp120-18A complex were subjected to 100 ns MD simulations. All-atom RMSD of both systems ( Figure 5) was $4 Å , indicating great structural changes during simulations. The radius of gyration and SASA profiles of gp120 and gp120-18A (S6) shows a more compact and stable gp120-18A compared to gp120 during the simulations, indicating that 18A stabilized the closed conformation of gp120. It is hypothesized that molecules that stabilize the closed gp120 conformation should be potential inhibitors (18). The average RMSD of the 18A was 1.283 Å indicating that 18A remained stable in the active site during the simulations, indicating that most structural fluctuations are due protein residues.  An average RMSF of 1.8 and 1.7 Å for gp120 and gp120-18A, respectively, indicates the more flexible nature of gp120 compared to the gp120-18A ( Figure 5B). It is observed that during the gp120 simulation b2 and b3 shift closer to b20/b21 loop, while the a1 begins to slightly unwind and also moves from starting position, pivoting away from a5. The b20/b21 loop begins to refold and shifts closer to a1 in gp120 simulations ( Figure 6). Spontaneous folding of the b20/b21 loop has been reported to be the main force in the restructuring of the unliganded core gp120 (19). It was observed that the interaction of 18A with Asn425, Gln428 restricted the refolding and movement of b20/b21 in gp120-18A simulations (see Figure 8). The presence of 18A enhanced fluctuations mainly in b22, b23, b26, b4, a0 regions that are disordered in gp120 structure, while there is a drastic reduction in fluctuations of a1 and b2 ( Figure 5).

Conformational behavior of unbound versus 18A-bound protein
MD simulations confirmed that the conformational behavior for gp120 and gp120-18A complex is quite distinct. The MD trajectory clearly showed gp120 adopts two conformations: open and closed ( Figure 6). Superposition of snapshots from gp120 simulations onto the starting conformation revealed that at 21 ns gp120 adopts a conformation with some the CD4bound structural changes, this is in great accordance with experimental findings (7). The structural changes observed included a change of perpendicular interactions of Trp112 and Trp427, inversion of Gln428, a1 changes, substantial movement of V1/V2 loops toward the V4 regions ( Figure 7) and opening of the main hydrophobic cavity (Figure 6). This observation substantiates the previous finding that gp120 is intrinsically dynamic and transitions between conformations, one of which is modeled by CD4 (18), contradicting the old paradigm that CD4 binds to gp120 and induces conformational changes. The full bridging sheet domain was not observed in this conformation, however, the observed drastic movement of b2 toward b21 and b3 toward b21 suggests that the bridging sheet was beginning to form but restricted by the V1/V2 and V3 domains as previously reported by Kwong et al. (17). Notably, this CD4-bound-like structural changes were not observed in snapshots from the gp120-18A simulation supporting the hypothesis that 18A inhibits HIV-1 entry by inhibiting some CD4-induced structural changes. On the other hand, gp120-18A complex shows that the binding cavity remained closed throughout the simulations (Figure 8).

V1/V2 domain (residues 126-196) consists of four antiparallel strands (A-D) forming a ''
Greek key'' domain as shown on the crystal structure ( Figure 9B). Compound 18A is reported to inhibit the function of gp120 by inhibiting the rearrangements of the V1/V2 regions induced by CD4 binding (16), however, the mechanism in not clear. The RMSF profile ( Figure 9A) for the V1/V2 regions (especially residues Leu165, Arg166, Asp167, Asn188) indicates that higher fluctuations of these regions in gp120 relative to gp120-18A complex. This sustains the hypothesis that the binding of 18A to gp120 alters the movements of V1/V2 loops.
Analysis of the MD trajectory for the gp120 showed that the unwinding of a1 and subsequently the slight inversion and shifts of b2 drives the rearrangements of V1/V2. During the simulations, V1/V2 strands A-D, and V2 loops move away from the V3 loop N-terminal regions toward the inner domain and outer domain, respectively. The hairpin connecting B and C strands dips down toward the V3 hairpin tip ( Figure 9B). This movement is more drastic in gp120 simulations relative to gp120-18A complex. In order to monitor the conformational events leading to the V1/V2 rearrangements, interatomic distances of Ca atoms between interacting residues: (i) Ser110 ... Gln114 and (ii) Thr123 ... Leu125 ( Figure 10C) driving the unwinding in a1 and slight inversion of b2, respectively, were monitored. The interatomic distance profile of Ser110 ... Gln114 ( Figure 10A) indicates that a1 remains stable and unwound for the gp120-18A simulation, while for the gp120 simulation the distance fluctuates from $8.5 Å at 20 ns up to $9 Å at 60 ns and remains high throughout the simulation. Closer examination of trajectories reveal that in the gp120-18A system, 18A interacts strongly with Trp427 on b20/b21 while Trp427 and Trp112 on a1 maintain strong contact (S7), thereby restricting the unwinding and shifts of a1. It was also found that the shift in a1 leads to a shift and slight inversion of b2 as evident by distance profile of Thr123 ... Leu125 ( Figure 10B). This distance was $6.5 Å for the gp120-18A system and remains stable throughout simulations compared to $7.5 Å for the gp120, indicating a high shift and change in topology of b2. The inversion of b2 drives the motion of the V1 loop and subsequently V2 loop as they pivot away from V3 N-terminal regions toward the outer domain ( Figure 10D). These findings clearly prove that 18A disrupts the movement of V1/V2 away from the V3 consistent with experimental evidence.
To further explore the motions of V1/V2 loops during MD simulations, we measured distances between C a atoms of interacting and conserved residues: (i) Leu179 ... Arg419 and (ii) Thr163 Gly312 ( Figure 11C). The Leu179 ... Arg419 distance profile ( Figure 11A) reveals that in case of gp120 the interatomic distances drop to $5 Å during at $2 0ns      indicating a stronger interaction between the V2 loop and the outer domain -this further confirms that the V1/V2 region pivots away from the V3 N-terminal portion to interact closely with the outer domain. For gp120-18A complex, the almost steady value of $9.5 Å suggests that the V1/V2 regions remain further away from the outer domain and closer to V3 N-terminal portion. Such findings clearly verify that the presence of 18A led to disruption in the rearrangements of V1/V2. Closer inspection of trajectories indicates that the hairpin between strands C and B of V1/V2 reorients downwards toward the V3 tip in both simulations. To further justify this, we monitored Thr163 ... Gly312 interatomic distances to elaborate on the nature of this motion ( Figure 11B). During gp120 simulations, the hairpin between strands C and B shifts further downwards to interact strongly with the V3 tip (average distance $5 Å ), while in gp120-18A simulations an average distance of $7.5 Å indicates weaker interactions between the hairpin and V3 tip ( Figure 11D).
Correlational movements in gp120 and gp120-18A Dynamic cross-correlation maps for gp120 and gp120-18A complex were plotted to provide a comparative correlative motions landscape in both protein systems ( Figure 12). As evident from Figure 12, more concerted motions in gp120 is observed when compared to gp120-18A complex. A drastic reduction in residue correlation is observed especially for the V1/V2 regions in the presence of 18A, substantiating the experimental findings (16). Significant correlations exist in the V1/V2 regions between the motions of a1/b2 and V1/V2 regions indicating the role of a1/b2 residues in the conformational flexibility of V1/V2. The correlation profiles also reveal fluctuations of V3 closely coupled to the V1/V2 regions.

Principal component analysis
Principal component analysis (PCA) was used to detect conformational states visited by gp120 and gp120-18A during the simulations. PCA plot ( Figure 13A) shows that gp120 samples a higher distribution of confromational states relative to the gp120-18A system confirming that 18A binding causes reduction in the structural flexibilty of protein as compared to the unbound state. The first principal component indicates the opening and closing of the hydrophobic cavity in gp120. Furthermore, procupine plots ( Figure 13B and C), which reveal the direction of motion of V1/V2 regions in both systems, clearly indicate a more open cavity in gp120 relative to gp120-18A complex.

Conclusions
In this report, a compilation of molecular dynamics analytical tools were applied to explore the dynamic motion of the V1/ V2 domain of gp120 protein and the mechanism by which the recently discovered broad-spectrum anti-HIV inhibitor, 18A, inhibit this protein. Experimental studies reported that 18A distrubs the rearrangment of V1/V2 domain of gp120, however, the molecular description of such conformational events remain obscure. Analyses of molecular dynamic trajectories reveal that the unwinding of a1 and slight inversion of b2 in gp120 induces a shift of VI/V2 domain away from the V3 N-terminal regions but toward the outer domain. Results also showed that 18A interacts with residues Leu179, Ile194, Ile424, Met426, Trp427, Asp368 and Met475 in the main hydrophobic active site cavity and disrupts the restructuring of V1/V2 domain. Stronger contacts between Trp425 and Trp112 rings, as a result of binding with 18A, may contribute to the reduced flexibility of a1, thus inhibiting the shifts of the V1/V2 region in gp120-18A complex. Closer observation of a snapshot obtained at 21 ns of gp120 simulations reveals a conformation that closely ressembles a CD4-bound-like conformation (observed inversion of Gln428, shift of Trp427 and a1 changes) consistent with the previous report by Munro et al. that gp120 transitions between conformations one of which is captured by CD4 (18). Notably, the full formation of the bridging sheet minidomain was not observed in the above conformation, consistent with studies of Kwong et al. that reported the restriction of bridging sheet domain formation by V1/V2 aand V3 regions. This CD4-like confromation was not observed in the gp120-18A simulations substatiating experimental evidence that 18A inhibits CD4 induced structural changes (16).
Information gained from this report not only provides insight into the various conformational states of the gp120 protein and sequential events driving or inhibiting restructuring of V1/V2 regions in gp120 as a result of binding to inhibitor, but also will asssist in the future design of potent HIV therapeutics.