Natural Polymorphisms D60E and I62V Stabilize a Closed Conformation in HIV-1 Protease in the Absence of an Inhibitor or Substrate

HIV infection remains a global health issue plagued by drug resistance and virological failure. Natural polymorphisms (NPs) contained within several African and Brazilian protease (PR) variants have been shown to induce a conformational landscape of more closed conformations compared to the sequence of subtype B prevalent in North America and Western Europe. Here we demonstrate through experimental pulsed EPR distance measurements and molecular dynamic (MD) simulations that the two common NPs D60E and I62V found within subtypes F and H can induce a closed conformation when introduced into HIV-1PR subtype B. Specifically, D60E alters the conformation in subtype B through the formation of a salt bridge with residue K43 contained within the nexus between the flap and hinge region of the HIV-1 PR fold. On the other hand, I62V modulates the packing of the hydrophobic cluster of the cantilever and fulcrum, also resulting in a more closed conformation.


Introduction
HIV-1 infection occurs worldwide and still suffers from drug resistance and virological failure [1].Treatment strategies for HIV-1 infection include a mixture of classes of drugs that target essential components of the HIV-1 lifecycle and are referred to as highly active antiretroviral therapy (HAART) [2].Although this approach is seeing great success and many people today live extended lifetimes, drug resistances emerge to many of the classes of drugs in HAART [3][4][5][6][7][8].
Of particular interest to us was the observation that the closed-conformation of HIV-1 PR, which is typically induced by substate/inhibitor binding, became dominant over the semi-open conformation for subtypes F and H in the absence of inhibitors [80].We also observed the dominance of the closed conformation with single amino acid substitutions A71V [25] and L63P [83] contained within subtype B.Here we report insights, from both pulsed EPR and computational studies, into the effects that natural polymorphisms D60E and I62V have on the PI-naive conformational landscape of subtype B HIV-1 PR.These residues were chosen based upon examination of the amino acid sequences of subtypes F and H, revealing their common NPs: D60E and I62V [80]. Figure 2 shows the sequence of the F and H variants and highlights the location of these substitutions in the cantilever of HIV-1 PR.Of particular interest to us was the observation that the closed-conformation of HIV-1 PR, which is typically induced by substate/inhibitor binding, became dominant over the semi-open conformation for subtypes F and H in the absence of inhibitors [80].We also observed the dominance of the closed conformation with single amino acid substitutions A71V [25] and L63P [83] contained within subtype B.Here we report insights, from both pulsed EPR and computational studies, into the effects that natural polymorphisms D60E and I62V have on the PI-naive conformational landscape of subtype B HIV-1 PR.These residues were chosen based upon examination of the amino acid sequences of subtypes F and H, revealing their common NPs: D60E and I62V [80]. Figure 2 shows the sequence of the F and H variants and highlights the location of these substitutions in the cantilever of HIV-1 PR.

Cloning and Site-Directed Mutagenesis
Site-directed mutagenesis reactions using the Quick-change kit were performed on a pET-23a plasmid encoding the LAI subtype B sequence shown in Figure 2 to generate plasmids encoding the D60E, I62V, and D60E/I62V variants.The subtype B sequence used was described previously and was codon-optimized and cloned between the NdeI and BamHI restriction sites [18,76,84].All constructs for experimental studies contained the following substitutions: D25N to inactive HIV-1PR to aid in protein stability by limiting autoproteolysis [14,78,81], Q7K, L33I, and L63I as stabilizing mutations that are commonly used within the subtype B sequence of an active enzyme to aid in preventing autolysis of HIV-1PR [84,85], C67A and C95A for the removal of native cysteine residues [18,86], and K55C for site-specific spin labeling for DEER experiments, which has been shown not to alter enzyme activity [18,87].The fidelity of each HIV-1 PR gene was confirmed by Sanger DNA sequencing (ICBR Genomics Facility, University of Florida).

Cloning and Site-Directed Mutagenesis
Site-directed mutagenesis reactions using the Quick-change kit were performed on a pET-23a plasmid encoding the LAI subtype B sequence shown in Figure 2 to generate plasmids encoding the D60E, I62V, and D60E/I62V variants.The subtype B sequence used was described previously and was codon-optimized and cloned between the NdeI and BamHI restriction sites [18,76,84].All constructs for experimental studies contained the following substitutions: D25N to inactive HIV-1PR to aid in protein stability by limiting autoproteolysis [14,78,81], Q7K, L33I, and L63I as stabilizing mutations that are commonly used within the subtype B sequence of an active enzyme to aid in preventing autolysis of HIV-1PR [84,85], C67A and C95A for the removal of native cysteine residues [18,86], and K55C for site-specific spin labeling for DEER experiments, which has been shown not to alter enzyme activity [18,87].The fidelity of each HIV-1 PR gene was confirmed by Sanger DNA sequencing (ICBR Genomics Facility, University of Florida).

Expression, Protein Purification, and Spin Labeling
Protein expression, purification, and spin labeling were carried out using protocols previously described [39,[78][79][80]83,85], with the pH of the inclusion body resuspension buffer adjusted to one unit above the isoelectric point (pI) of each construct, which for all three variants studied here is 9.39.The spin label (1-oxyl-2,2,5,5-tetramethyl-d3-pyrrolidine-3-methyl)methanethiosulfonate (MTSL) (Santa Cruz Biotechnology, Santa Cruz, CA, USA) was freshly dissolved into ethanol (~1 mg into 100 µL) and an appropriate volume added to achieve a 4:1 SL:P molar excess during the labeling reaction, which was carried out at 4 • C (to prevent protease precipitation) in the dark overnight.The protein precipitate was removed by centrifugation (17,210× g, 25 • C).The excess free spin label was removed, and protein was buffered to 2 mM NaOAc, pH 5.0, using a HiPrep 26/10 desalting column (GE Healthcare, Chicago, IL, USA).

DEER Experiments and Data Analysis
HIV-1 PR samples were concentrated to 140 µM and buffer exchanged to 20 mM d 3 -NaOAc/D 2 O, pH 5.0, using centrifugal membrane concentrators (Millipore, Billerica, MA, USA).The protein sample was then transferred to a 4 mm quartz EPR tube (Norell, Marion, NC, USA) and flash-frozen in liquid nitrogen.The sample was then immediately inserted into the ER 4118X-MD-5 dielectric ring resonator.DEER experiments were carried out at 65 K using a four-pulse DEER sequence [69,72].The DEER modulation curves were background-corrected and converted to distance distribution profiles via Tikhonov regularization using DeerAnalysis [88].To obtain the fractional occupancy of conformational states of HIV-1PR, we deconstruct the DEER distance profiles into a linear superposition of Gaussian-shaped populations using DEERconstruct (a program developed in our laboratory) [89,90] to represent the conformational states of curled/tucked, closed, semi-open, and wide-open states.

Molecular Dynamics Simulations
MD simulations were performed using the HiPerGator2.0supercomputer at the University of Florida.The starting coordinates of HIV-1PR LAI were taken from PDBID: 1HHP and converted to a dimer [17].The amino acid substitutions described above and outlined in Figure 1 that are not present in the 1HHP sequence that are included in simulations are the following: Q7K, D25N, L33I, C67A, and C95A for subtype B. The amino acid substitutions were incorporated using the tleap module of the AMBER suite of programs.Starting structures for HIV-1PR D60E, I62V, and D60E/I62V were generated from the above subtype B sequence.The starting structure for HIV-1PR bound to Darunavir for simulations of a closed conformation was taken from PDBID: 3BVB [14].Simulations were set up using the procedure previously published [91][92][93][94][95].All molecular dynamics (MD) simulations employed the ff99SB forcefield and the AMBER suite of programs [96,97].Prior to starting simulations, each protein system was run through the H++ protonation state server [98] to model the charged amino acids based on their protonation states [99].Each system was then solvated in a rectilinear periodic box of explicit SPC/E water molecules set at pH 6.8 [100], with the box extending 8 Å from the solute.Counterions were added to net-neutralize the systems.For all simulations, the proteins went through a five-stage minimization followed by a two-stage equilibration prior to running production.In the first four stages of minimization, the following were allowed to be energy minimized in the following order: all water molecules and counterions, all hydrogen atoms, side chain groups, and backbone amide groups.In the last stage of minimization, the entire system, including proteins, water molecules, and counterions, was allowed to be energy minimized through a progressive release of constraints.The systems were slowly heated to 300 K over 100 ps for the canonical ensemble.In the last stage of equilibration, 1 ns of MD was performed at 300 K for an isothermal and isobaric (NPT) ensemble.The simulations were then run at 300 K in the NPT ensemble.A Langevin thermostat with a 2 ps −1 collision frequency was employed to maintain the temperature of the system [101,102].The SHAKE algorithm [103] was used to restrain hydrogen bond lengths, and the particle mesh Ewald method was used to calculate long-range electrostatic interactions.A time step of 2 fs was utilized.A total of 1.6 µs of MD data were collected for each system.MD analyses were performed using the ptraj tool of AMBERTools [96].

D60E and I62V Shift the HIV-1 PR Conformational Sampling to Closed-like States
To characterize the impacts of these NPs on the conformational landscape of the unbound enzyme, DEER analysis was performed for WT, the two single and the double substituted HIV-1 PR constructs, to obtain the K55C-MTSL interspin distance distribution probability profiles given in Figure 3.Because HIV-1PR is a homodimer, the incorporation of a single spin-label site generates a spin-pair for distance measurements that can be used to describe the conformational landscape of HIV-1PR [18,27,39,89].Overall, D60E is found to alter the conformational sampling more dramatically than the I62V substitution.The most probable K55C interspin distance in subtype B is 36.2 ± 0.2 Å, which corresponds to a predominant semi-open conformation (Figure 2A) [76].D60E shifts the most probable interspin distance to 34.0 ± 0.2 Å, like the effects we observed previously for other single amino acid substitutions, L63P [95] and A71V [25] (Table 1).In contrast, a less dramatic shift in the most probable distance is observed for I62V, with a most probable interspin distance of 35.2 ± 0.2 Å, which occurs between the semi-open and closed distances.When both substitutions are included, the conformational landscape contains a marked bimodal shape with two major peaks at 34.2 ± 0.2 Å and 36.5 ± 0.2 Å.
collision frequency was employed to maintain the temperature of the system [101,102].The SHAKE algorithm [103] was used to restrain hydrogen bond lengths, and the particle mesh Ewald method was used to calculate long-range electrostatic interactions.A time step of 2 fs was utilized.A total of 1.6 µs of MD data were collected for each system.MD analyses were performed using the ptraj tool of AMBERTools [96].

D60E and I62V Shift the HIV-1 PR Conformational Sampling to Closed-like States
To characterize the impacts of these NPs on the conformational landscape of the unbound enzyme, DEER analysis was performed for WT, the two single and the double substituted HIV-1 PR constructs, to obtain the K55C-MTSL interspin distance distribution probability profiles given in Figure 3.Because HIV-1PR is a homodimer, the incorporation of a single spin-label site generates a spin-pair for distance measurements that can be used to describe the conformational landscape of HIV-1PR [18,27,39,89].Overall, D60E is found to alter the conformational sampling more dramatically than the I62V substitution.The most probable K55C interspin distance in subtype B is 36.2 ± 0.2 Å, which corresponds to a predominant semi-open conformation (Figure 2A) [76].D60E shifts the most probable interspin distance to 34.0 ± 0.2 Å, like the effects we observed previously for other single amino acid substitutions, L63P [95] and A71V [25] (Table 1).In contrast, a less dramatic shift in the most probable distance is observed for I62V, with a most probable interspin distance of 35.2 ± 0.2 Å, which occurs between the semi-open and closed distances.When both substitutions are included, the conformational landscape contains a marked bimodal shape with two major peaks at 34.2 ± 0.2 Å and 36.5 ± 0.2 Å.   is the same as reported within ref. [76].Full data analysis is provided in supporting information (Figures S1-S7).± 5 a data taken from ref. [76]; b data taken from ref. [80]; c data taken from ref. [25]; d data taken from ref. [95].
Viruses 2024, 16, 236 6 of 18 The DEER conformational landscape represents the summation of the four previously mentioned HIV-1 PR conformations, specifically the inhibitor-bound closed state, the semi-open state, the wide-open state, and the curled open/tucked flap states [18,27,39,89].Through several investigations, we identified the distances between the K55 spin pairs of each of the four conformational states as 25-30 Å for a flap-curled open/tucked conformation, 33 Å for the closed-conformation, 36 Å for the semi-open conformation, and >40 Å for the wide-open conformations [25,27,79,80,95].The fractional occupancy of a given conformation is obtained by fitting the DEER distance profile to a linear combination of Gaussian-shaped populations that corresponds to the four above-described conformations of HIV-1 PR [89,90].The Gaussian distribution diagrams for the DEER distance profiles in Figure 3A are given in Figures S1-S7.The relative percentages of each Gaussian population, a.k.a.population analysis, are shown in Figure 3B and Table 1, with the relative differences in population to subtype B shown in Figure 3C.Analysis reveals that all constructs have greater occupancy of the closed conformation (65%, 49%, and 59% for D60E, I62V, and D60E/I62V, respectively) compared to the subtype B conformational landscape and demonstrate that the closed conformation is the dominant conformation for each construct.
We previously studied the impact of other cantilever mutations.The mutation L/I63P substitution also shifts the conformation to the closed state (fractional occupancy 57%), increases backbone rigidity, and induces a closed structure.The flaps demonstrated the handedness of an inhibitor-induced closed structure [PDB ID: 5T84] when crystallized in the absence of any substrate analog or inhibitor [95].Additionally, the single substitution A71V also dramatically shifts the conformational sampling to a predominantly closed state (67% fractional occupancy) [25].The fractional occupancy of the wide-open state did not change for D60E and D60E/I62V relative to B, whereas no population was observed for wide-open conformation in I62V.Furthermore, DEER distance profiles for constructs had a statistically relevant population (within the SNR of the data, Supporting Information) below 30 Å (10-12%) that we assigned to a curled/tucked conformation [81,82,95], displaying a higher population than observed in subtype B (4 ± 4%).Although the increase in the fractional occupancy of the curled/tucked conformation was less than that seen for the closed conformation, results reveal that D60E and I62V substitutions lead to an increase in flap curling.Additionally, upon complexation with DRV, these populations are retained within the "bound" distance profiles (Figure S1-S7).

I62V Participates in the Hydrophobic Siding Mechanism
As secondary mutations arise in response to binding with PIs, they are generally located in the core of HIV-1 PR away from the catalytic pocket, yet they may still indirectly facilitate or contribute to conformational changes [30,80,82,104,105].The core of HIV-1 PR, i.e., the fulcrum and cantilever, contains 19 hydrophobic residues outside of the catalytic pocket [105][106][107][108].Of these 19 hydrophobic residues, nine are isoleucine, which has more possible conformations than other hydrophobic amino acids.When a mutation in the core arises, it may disrupt the hydrogen-bonding network or the van der Waal interaction network within the core.With many isoleucine residues in the core, the side chain of a residue can change its conformation, causing the hydrophobic surface to slide from one hydrophobic residue to another by exchanging van der Waal contact, which in turn facilitates conformational change at minimal energetic expense [106].
The residues participating in the hydrophobic sliding mechanism are shown in Figure 4. Structurally, sites 63, 71, 62, 64, 66, 89, and 93 interact with one another within the hydrophobic core cluster [106].Subtypes F and H variants both have four NPs in the hydrophobic core, one of which is I62V.Mutation of the hydrophobic core residues is more likely to affect the dynamic properties of HIV-1 PR and inhibitor binding but not substrate binding [106,108].Thus, it is reasonable that I63V substitutions can modulate the hydrophobic packing/sliding in this cluster, resulting in an induced shift in the stability of the conformational landscape.The combination of D60E and I62V generates an aver-aged impact on the shift to the closed conformation.The offsetting effects of the double substitution are discussed further below.
hydrophobic core cluster [106].Subtypes F and H variants both have four NPs in the hydrophobic core, one of which is I62V.Mutation of the hydrophobic core residues is more likely to affect the dynamic properties of HIV-1 PR and inhibitor binding but not substrate binding [106,108].Thus, it is reasonable that I63V substitutions can modulate the hydrophobic packing/sliding in this cluster, resulting in an induced shift in the stability of the conformational landscape.The combination of D60E and I62V generates an averaged impact on the shift to the closed conformation.The offsetting effects of the double substitution are discussed further below.

Root-Mean-Square Deviation during Simulations
MD simulations from HIV-1PR subtype B (labeled as WT), D60E, I62V, and D60E/I62V were performed to further study the effects these amino acid substitutions have on the dynamics and conformational flexibility of HIV-1 PR. Figure 5 shows the calculated root-mean-square deviation (RMSD) of the backbone in the flap regions (residues 43-59) of the four constructs over the course of the simulation.The RMSD pattern for WT and D60E/I62V is similar, with values near 2 Å for the first 800 ns and a jump to 4.5-5 afterwards.Although similar, the D60E/I62V has RMSD values that fluctuate more frequently than do those for WT.Specifically, there is a jump in RMSD from 2 Å to 5.5 Å around 600-700 ns, followed by a drop in RMSD of 1.5 Å around 700-800 ns.RMSD of the D60E/I62V then increases to ~5 Å, like the behavior observed in WT, but then has additional fluctuations between 2-5 Å over the 800-1300 ns range.Between 1300 and 1550 ns, the RMSD for D60E/I62V decreased from 5 Å to 3.5 Å, like the behavior seen for D60E/I62V in this period.Finally, the RMSD value increased to 5 Å in the last 50 ns.The RMSD of both single mutants, D60E and I62V, fluctuates largely from 1-6 Å in the first 500-600 ns and stabilizes to a value of ~ 3.5 Å during the last 900 ns, which is different from WT.

Root-Mean-Square Deviation during Simulations
MD simulations from HIV-1PR subtype B (labeled as WT), D60E, I62V, and D60E/I62V were performed to further study the effects these amino acid substitutions have on the dynamics and conformational flexibility of HIV-1 PR. Figure 5 shows the calculated rootmean-square deviation (RMSD) of the backbone in the flap regions (residues 43-59) of the four constructs over the course of the simulation.The RMSD pattern for WT and D60E/I62V is similar, with values near 2 Å for the first 800 ns and a jump to 4.5-5 afterwards.Although similar, the D60E/I62V has RMSD values that fluctuate more frequently than do those for WT.Specifically, there is a jump in RMSD from 2 Å to 5.5 Å around 600-700 ns, followed by a drop in RMSD of 1.5 Å around 700-800 ns.RMSD of the D60E/I62V then increases to ~5 Å, like the behavior observed in WT, but then has additional fluctuations between 2-5 Å over the 800-1300 ns range.Between 1300 and 1550 ns, the RMSD for D60E/I62V decreased from 5 Å to 3.5 Å, like the behavior seen for D60E/I62V in this period.Finally, the RMSD value increased to 5 Å in the last 50 ns.The RMSD of both single mutants, D60E and I62V, fluctuates largely from 1-6 Å in the first 500-600 ns and stabilizes to a value of ~3.5 Å during the last 900 ns, which is different from WT.

D60E Mutation Facilitates a Salt Bridge Interaction
Site D60 within HIV-1PR sits in the cantilever region (residues 60-75), as shown in Figures 1 and 6A.When comparing the crystal structures of PI-naïve subtype B (i.e., WT) and subtype F (which contains the D60E substitution, Figure 6B), one can readily see a salt bridge between D60E-K43 in the subtype F structure (3.0 Å) that is altered in WT (4.0 Å) due to the shorter length of Asp (D) compared to GLU (E).The D60E-K43 bridge connects the cantilever (D60E) to the flap/hinge nexus (site K43) in HIV-1PR.It has been shown in the literature that salt bridge interactions in the hinge region are important and can affect mobility as well as the conformational landscape of HIV-1 PR.One such salt bridge interaction that was previously studied in our laboratories is E35D-R57K, which connects the flap region (R57K) to the elbow/hinge region (E35D), which we showed can modulate the relative population of the curled/tucked flap conformation [81,82].Additionally, Schiffer and colleagues demonstrated the impact of reversible disulfide linkages between the elbow/hinge and flap (L38C/G16C) and the cantilever and flap (E65C/G16C) [108].These crosslinking studies showed that restriction of motion by crosslinking the flap to the elbow/hinge diminished catalytic activity (~150 fold lower than WT), whereas the linkage between the cantilever and flap did not impact catalytic activity (~1× WT) [108].

D60E Mutation Facilitates a Salt Bridge Interaction
Site D60 within HIV-1PR sits in the cantilever region (residues 60-75), as shown in Figures 1 and 6A.When comparing the crystal structures of PI-naïve subtype B (i.e., WT) and subtype F (which contains the D60E substitution, Figure 6B), one can readily see a salt bridge between D60E-K43 in the subtype F structure (3.0 Å) that is altered in WT (4.0 Å) due to the shorter length of Asp (D) compared to GLU (E).The D60E-K43 bridge connects the cantilever (D60E) to the flap/hinge nexus (site K43) in HIV-1PR.It has been shown in the literature that salt bridge interactions in the hinge region are important and can affect mobility as well as the conformational landscape of HIV-1 PR.One such salt bridge interaction that was previously studied in our laboratories is E35D-R57K, which connects the flap region (R57K) to the elbow/hinge region (E35D), which we showed can modulate the relative population of the curled/tucked flap conformation [81,82].Additionally, Schiffer and colleagues demonstrated the impact of reversible disulfide linkages between the elbow/hinge and flap (L38C/G16C) and the cantilever and flap (E65C/G16C) [108].These crosslinking studies showed that restriction of motion by crosslinking the flap to the elbow/hinge diminished catalytic activity (~150 fold lower than WT), whereas the linkage between the cantilever and flap did not impact catalytic activity (~1× WT) [108].
Analysis of the MD simulation runs shows that the relative population of this salt bridge is significantly increased in WT when D60E is included.Analyses of MD trajectories for WT and I62V reveal a broad D60-K43 N-O distance, spanning from 3-10 Å with a major peak at 7.5-8.5Å and a smaller shoulder at 5.2 Å (Figure 6D).The presence of a very small peak at 3.0 Å indicates only a slight tendency for a salt bridge interaction in these two constructs.This interaction is observed within the WT crystal structure (PDB ID 1HHP), but only moderately populated when molecular dynamics are included, which highlights the potential drawbacks of looking exclusively at static crystal structures that may stabilize conformations due to crystal packing.In contrast, simulation runs for constructs that contain the D60E substitution (D60E and D60E/I62V) show a dramatic increase in the probability density at 3.0 Å, indicating the D60E-K43 salt bridge is much more favorable when residue 60 is altered to Glu.Additionally, the 7.5-8.5Å peak seen in WT and I62V has significantly reduced, and the shoulder at 5.2 Å shifts to shorter values of 4.5 Analysis of the MD simulation runs shows that the relative population of this salt bridge is significantly increased in WT when D60E is included.Analyses of MD trajectories for WT and I62V reveal a broad D60-K43 N-O distance, spanning from 3-10 Å with a major peak at 7.5-8.5Å and a smaller shoulder at 5.2 Å (Figure 6D).The presence of a very small peak at 3.0 Å indicates only a slight tendency for a salt bridge interaction in these two constructs.This interaction is observed within the WT crystal structure (PDB ID 1HHP), but only moderately populated when molecular dynamics are included, which highlights the potential drawbacks of looking exclusively at static crystal structures that may stabilize conformations due to crystal packing.In contrast, simulation runs for constructs that contain the D60E substitution (D60E and D60E/I62V) show a dramatic increase in the probability density at 3.0 Å, indicating the D60E-K43 salt bridge is much more favorable when residue 60 is altered to Glu.Additionally, the 7.5-8.5Å peak seen in WT and I62V has significantly reduced, and the shoulder at 5.2 Å shifts to shorter values of 4.5 Å in the presence of D60E.
The intensity of the 3.0 Å peak in the D60E/I62V simulation is lower than that in WT.D60E/I62V also shows a slight shift of the 4.5 Å peak to higher distances and a higher intensity for distances greater than 5.0 Å.In WT, I62 Cα has van der Waals contacts with Y59 Cα and L38 Cα1 at 3.6 Å and 4.1 Å, respectively (Figure 6C).When mutated to Val, the side chain of residue 62 is shortened by one carbon bond, and interactions with L38 and Y59 may be lost, possibly resulting in more flexibility in the hinge region, which may explain the lowered intensity observed for the 3.0 Å peak in I62V simulations.Perhaps the addition of the I62V mutation to the D60E mutation slightly destabilizes the D60E-K43 salt bridge to retain flap flexibility and help restore the semi-open population, as seen in DEER results (Table 1).
the side chain of residue 62 is shortened by one carbon bond, and interactions with L38 and Y59 may be lost, possibly resulting in more flexibility in the hinge region, which may explain the lowered intensity observed for the 3.0 Å peak in I62V simulations.Perhaps the addition of the I62V mutation to the D60E mutation slightly destabilizes the D60E-K43 salt bridge to retain flap flexibility and help restore the semi-open population, as seen in DEER results (Table 1).

Distance Measurements for Conformational Assignment in Simulations
To better define conformations probed via the flap K55R1-K55R1′ distance distribution profiles seen with DEER spectroscopy, we evaluated the K55C-K55C′ Cα distance together with the I50-I50′ Cα distance of the flap tips to generate distribution profiles to assess and describe conformations sampled during the MD trajectories.The locations of these residues in the flaps of HIV-1PR are shown in Figure 7A.With regards to the K55-K55′ Cα (Figure 7B) distance probability profile, WT protease displays three distances, with a major peak at 25 Å (assigned as semi-open) and two minor peaks at 19 Å and 21 Å,

Distance Measurements for Conformational Assignment in Simulations
To better define conformations probed via the flap K55R1-K55R1 ′ distance distribution profiles seen with DEER spectroscopy, we evaluated the K55C-K55C ′ Cα distance together with the I50-I50 ′ Cα distance of the flap tips to generate distribution profiles to assess and describe conformations sampled during the MD trajectories.The locations of these residues in the flaps of HIV-1PR are shown in Figure 7A.With regards to the K55-K55 ′ Cα (Figure 7B) distance probability profile, WT protease displays three distances, with a major peak at 25 Å (assigned as semi-open) and two minor peaks at 19 Å and 21 Å, which we label a curling or tucking of flaps that represent an open-like state [25,79].When each single amino acid substitution is incorporated, there is a dramatic change observed in the distance distribution profile compared to WT.Specifically, a new major peak appears at 23 Å (assigned as closed) that was not seen in WT, with diminished or no intensity for peaks at 25 Å, 21 Å, and 19 Å.When D60E and I62V are combined, the 25 Å peak regains similar intensity as observed in WT but also retains some probability density at 23 Å and 19-20 Å.We assign these K55-K55 ′ Cα peaks to the four conformations of HIV-1 PR seen in which we label a curling or tucking of flaps that represent an open-like state [25,79].When each single amino acid substitution is incorporated, there is a dramatic change observed in the distance distribution profile compared to WT.Specifically, a new major peak appears at 23 Å (assigned as closed) that was not seen in WT, with diminished or no intensity for peaks at 25 Å, 21 Å, and 19 Å.When D60E and I62V are combined, the 25 Å peak regains similar intensity as observed in WT but also retains some probability density at 23 Å and 19-20 Å.We assign these K55-K55′ Cα peaks to the four conformations of HIV-1 PR seen in DEER data analyses and X-ray structures as:  Furthermore, I50-I50′ Cα distance distribution profiles (Figure 7C) show a similar effect.Specifically, the single amino acid substitutions markedly alter the distribution profile of WT, with diminished probability at 5 Å (the major peak observed for WT) and the growth of new probability peaks at 9.4 Å and 8.5 Å for D60E and I62V, respectively, that are absent in WT data.The combination of the two substitutions again results in an averaged distance distribution profile where the probability at 5 Å grows towards WT values but where D60E/I62V retains the new distance of ~9 Å.
The conformations related to I50-I50′ Cα distances are assigned based on a survey of the I50-I50′ Cα distance from various crystal structures of HIV-1 PR [16,34,35,49,109,110] (Table S1), and our MD trajectories are discussed more below, including that of WT protease bound to Darunavir (Figure S8), which is known to stabilize a closed structure [14].When protease is in an inhibitor-bound closed conformation, the flap tips cross over one another.As the enzyme opens to a semi-open state, the distance between sites I50 and I50′ in the tips is closer to one another than when the inhibitor is bound in the closed state, resulting in a shortened distance between the I50-I50′ Cα atoms.The flap and flap tips Furthermore, I50-I50 ′ Cα distance distribution profiles (Figure 7C) show a similar effect.Specifically, the single amino acid substitutions markedly alter the distribution profile of WT, with diminished probability at 5 Å (the major peak observed for WT) and the growth of new probability peaks at 9.4 Å and 8.5 Å for D60E and I62V, respectively, that are absent in WT data.The combination of the two substitutions again results in an averaged distance distribution profile where the probability at 5 Å grows towards WT values but where D60E/I62V retains the new distance of ~9 Å.
The conformations related to I50-I50 ′ Cα distances are assigned based on a survey of the I50-I50 ′ Cα distance from various crystal structures of HIV-1 PR [16,34,35,49,109,110] (Table S1), and our MD trajectories are discussed more below, including that of WT protease bound to Darunavir (Figure S8), which is known to stabilize a closed structure [14].When protease is in an inhibitor-bound closed conformation, the flap tips cross over one another.As the enzyme opens to a semi-open state, the distance between sites I50 and I50 ′ in the tips is closer to one another than when the inhibitor is bound in the closed state, resulting in a shortened distance between the I50-I50 ′ Cα atoms.The flap and flap tips move furthest from each other in the wide-open conformation.Therefore, we assign these I50-I50 ′ Cα peaks to the HIV-1 PR conformations: closed (8.5 Å), semi-open (5 Å), and wide-open (13-16 Å).When the flaps curl or tuck open, the I50-I50 ′ Cα distance increases; therefore, the population observed for >11 Å is assigned as a combination of the open-like states.We also propose (see double distance plots in Figure 8) that the flap tip distance of 6 Å for I50-I50 ′ Cα also corresponds to curled/tucked conformation.We computationally determined the closed state MD K55-K55 ′ Cα and I50-I50 ′ Cα distance distribution profiles from analysis of the trajectories of WT-bound to DRV, in agreement with ~8 Å and ~6 Å for the I50-I50 ′ Cα distance distribution of the closed and curled/tucked conformations, respectively.For the DRV-bound HIV-1PR simulations, the K55-K55 ′ Cα is centered at 22.8 Å (Figure S8).This agrees with the K55-K55 ′ Cα data: WT assembles different conformations.When the D60E mutation is introduced, the conformation of the PR is shifted to closed.The I62V mutation also shifted conformation to closed, but to a lesser extent.When both mutations are introduced, the conformations seen in WT are regained while also retaining a portion of the closed conformation.
move furthest from each other in the wide-open conformation.Therefore, we assign these I50-I50′ Cα peaks to the HIV-1 PR conformations: closed (8.5 Å), semi-open (5 Å), and wide-open (13-16 Å).When the flaps curl or tuck open, the I50-I50′ Cα distance increases; therefore, the population observed for >11 Å is assigned as a combination of the open-like states.We also propose (see double distance plots in Figure 8) that the flap tip distance of 6 Å for I50-I50′ Cα also corresponds to curled/tucked conformation.We computationally determined the closed state MD K55-K55′ Cα and I50-I50′ Cα distance distribution profiles from analysis of the MD trajectories of WT-bound to DRV, in agreement with ~8 Å and ~6 Å for the I50-I50′ Cα distance distribution of the closed and curled/tucked conformations, respectively.For the DRV-bound HIV-1PR simulations, the K55-K55′ Cα is centered at 22.8 Å (Figure S8).This agrees with the K55-K55′ Cα data: WT assembles different conformations.When the D60E mutation is introduced, the conformation of the PR is shifted to closed.The I62V mutation also shifted conformation to closed, but to a lesser extent.When both mutations are introduced, the conformations seen in WT are regained while also retaining a portion of the closed conformation.

Alterations in the Landscape by Double Distance Measurement
Even though Figure 7 shows probability distance profiles of flaps and flap tips for conformers sampled throughout the simulation, time information is not captured within Figure 7, which makes it difficult to correlate these distances and assign conformations.Double-distance plots showing both distances, flaps (K55-K55 ′ ), and flap tips (I50-I50 ′ ), in the same snapshot, can better reveal the conformational landscape of HIV-1 PR throughout the simulation (Figure 8). Figure 8A shows WT assembled predominantly from three different conformational populations represented by three different areas labeled 1-3, with area 3 representing the semi-open conformation.There is also a smaller area labeled 4, representing a fourth conformational population with a much lower probability compared to the other three conformations.During the MD runs for WT, the wide-open conformation was not sampled.Instead, access to the active site appeared due to the flap curling/tucking as evidenced by populations 1 and 2 with short K55-K55 ′ flap distances and I50-I50 ′ flap tip distances spanning 6 to 18 Å.
The conformational landscape of the D60E single mutant differs dramatically from WT. D60E sampled larger flap distances (K55-K55 ′ ) above 30 Å, which was not the case for WT.In D60E, conformations 1 and 2 are absent, with the appearance of a new conformational population represented by area 5.This area 5 had a flap distance of 21-25 Å and a flap tip distance of 6-9 Å, which we assigned to the closed conformation as these distances were observed in the MD trajectories of the DRV-bound subtype B (Figure S8).This result is consistent with the experimental DEER data, where the D60E substitution increases the closed conformation percentage in HIV-1 PR.
The double distance plot of I62V shows a similar ensemble pattern compared to WT.The distance pairs, however, spread a wide range of distance combinations within the conformational landscape instead of assembling into more defined areas, as seen in the WT and D60E single mutants.Nevertheless, closed conformation represented by area 5 was also present, as evidenced by data points at 21-24 Å for flap distance and 6-9 Å for flap tip distance (Figure 8A).The blue density in Figure 8B at the position of area 1 showed that the conformation represented by area 1 was decreased to a lesser degree compared to the difference plot of D60E-WT (Figure 8B).In other words, the conformation for area 1 in I62V was decreased to give rise to the closed conformation represented by area 5.
A similar conformational ensemble to the I62V single mutant was observed for the D60E/I62V double mutant.The conformational ensemble of D60E/I62V was spread over a wide range of distance combinations, even more like WT than I62V.In contrast to the double distance plot of both single mutants, area 5 was not clearly observed in the double distance plot of D60E/I62V.The data points for areas 5 (increase in closed conformation) and 1 (decrease in flap curling/tucking) in Figure 8B were even less intense compared to the I62V single mutant.These changes in the conformational landscape of the double mutant seemed to be minimal, and the overall conformational landscape is shown to be very similar to that of WT.

Average Structure during Simulations
The average coordinates of each construct over the course of 1.6 µs simulations were also calculated, and the average structures are shown in Figure 9.In general, the average structures did not have large variation in other regions of the PR except for the flap region, especially flap tips.The average structure of WT is in semi-open conformation, with one flap raised higher than the other flap and slightly curled flap tips.On the other hand, both single mutants adopted the closed conformation on average, with their flap tips rather straight (Figure S9).The average structure of D60E/I62V has more in common with the average structure of WT.On average, the double mutant adopted more of a semi-open conformation as well as curled flap tips, like the average structure of WT.The double mutant had one flap slightly raised higher than that of either single mutant but not as high as the WT.The combination of both mutations resulted in an average structure halfway between WT and mutants.The flap tips of double mutants were also more curled compared to the flap tips of WT.
flap raised higher than the other flap and slightly curled flap tips.On the other hand, both single mutants adopted the closed conformation on average, with their flap tips rather straight (Figure S9).The average structure of D60E/I62V has more in common with the average structure of WT.On average, the double mutant adopted more of a semi-open conformation as well as curled flap tips, like the average structure of WT.The double mutant had one flap slightly raised higher than that of either single mutant but not as high as the WT.The combination of both mutations resulted in an average structure halfway between WT and single mutants.The flap tips of double mutants were also more curled compared to the flap tips of WT.The flap-handedness of WT is the same as in the starting structure and crystal structure of B (PDB ID: 1HHP).An inversion in flap-handedness was observed for both single mutants.Not only did the double mutant adopt similar flap conformation as WT on average, but it also had the same flap-handedness as WT.This agrees with our hypothesis that the I62V mutation serves as a compensatory mutation in the presence of the D60E mutation to regain conformation of WT. Figure 9 also shows that the flap conformation as well as the flap-handedness of WT and the double mutant are similar in their starting structure.On the other hand, the flap conformation and flap-handedness of the single mutants are similar to each other and different from the starting structure.

Conclusions
The distance profiles of the single and double mutants in the unbound state confirm that the D60E and I62V mutations shift protein conformation to a higher closed conformation percentage.This increase was also observed in MD simulations.This result explains the most probable closed conformation distance of the F and H variants presented in our previous work [80].MD simulations reveal a salt bridge interaction between D60E and K43 was facilitated in the presence of the D60E mutation to stabilize a closed flap conformation.I62V likely participates in a hydrophobic sliding mechanism, resulting in a change in conformation.Together with other NPs present in F and H variants, D60E and I62V indirectly contribute to conformational change and possibly affect binding to PIs other than DRV [80].

Supplementary Materials:
The following supporting information can be downloaded at: www.mdpi.com/xxx/s1.Figure S1: DEER data processing for D60E-unbound.The flap-handedness of WT is the same as in the starting structure and crystal structure of B (PDB ID: 1HHP).An inversion in flap-handedness was observed for both single mutants.Not only did the double mutant adopt similar flap conformation as WT on average, but it also had the same flap-handedness as WT.This agrees with our hypothesis that the I62V mutation serves as a compensatory mutation in the presence of the D60E mutation to regain conformation of WT. Figure 9 also shows that the flap conformation as well as the flap-handedness of WT and the double mutant are similar in their starting structure.On the other hand, the flap conformation and flap-handedness of the single mutants are similar to each other and different from the starting structure.

Conclusions
The distance profiles of the single and double mutants in the unbound state confirm that the D60E and I62V mutations shift protein conformation to a higher closed conformation percentage.This increase was also observed in MD simulations.This result explains the most probable closed conformation distance of the F and H variants presented in our previous work [80].MD simulations reveal a salt bridge interaction between D60E and K43 was facilitated in the presence of the D60E mutation to stabilize a closed flap conformation.I62V likely participates in a hydrophobic sliding mechanism, resulting in a change in conformation.Together with other NPs present in F and H variants, D60E and I62V indirectly contribute to conformational change and possibly affect binding to PIs other than DRV [80].

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/v16020236/s1.Figure S1: DEER data processing for D60Eunbound.Figure S2: DEER data processing for the D60E-DRV.Figure S3: DEER data processing for I62V-unbound.Author Contributions: T.T.T. performed all experimental and computational aspects of the manuscript.G.E.F. and T.T.T. prepared the drafts and final version of the manuscript.G.E.F.oversaw all aspects of the experimental design, data analysis, and interpretation of results.All authors have read and agreed to the published version of the manuscript.

Figure 1 .
Figure 1.Ribbon diagram of a highly drug-resistant HIV-1 PR, PRS17, in complex with substrate analog CA-P2 [PDB ID 6O5X], showing the flaps closed around the substrate analog/inhibitor in the active site, with active site aspartic acid residues D25 and D25′ rendered in ball and stick format.Several structural regions discussed within are rendered in varying colors, such as the "flap region" (residues 43-59) in green, the "hinge/elbow region" (residues 35-42) in bright pink, the "fulcrum region" (residues 9-24) in yellow, the "cantilever region" (residues 60-75), and the "near-active site region" (residues 80-84) in cyan.

Figure 1 .
Figure 1.Ribbon diagram of a highly drug-resistant HIV-1 PR, PRS17, in complex with substrate analog CA-P2 [PDB ID 6O5X], showing the flaps closed around the substrate analog/inhibitor in the active site, with active site aspartic acid residues D25 and D25 ′ rendered in ball and stick format.Several structural regions discussed within are rendered in varying colors, such as the "flap region" (residues 43-59) in green, the "hinge/elbow region" (residues 35-42) in bright pink, the "fulcrum region" (residues 9-24) in yellow, the "cantilever region" (residues 60-75), and the "near-active site region" (residues 80-84) in cyan.

Figure 2 .
Figure 2. Comparison of HIV-1 PR sequences of subtypes B, F, and H. (A) Ribbon diagram of the closed conformation (PDB ID 2PBX) showing the locations of natural polymorphisms (NPs) and (B) sequence alignment comparisons.All NPs are represented as spheres, with D60E in royal blue, I62V in red, hydrophobic core residues in teal, and other NPs in gray.Hydrophobic core residues are underlined in the PI-naive subtype B sequence.

Figure 2 .
Figure 2. Comparison of HIV-1 PR sequences of subtypes B, F, and H. (A) Ribbon diagram of the closed conformation (PDB ID 2PBX) showing the locations of natural polymorphisms (NPs) and (B) sequence alignment comparisons.All NPs are represented as spheres, with D60E in royal blue, I62V in red, hydrophobic core residues in teal, and other NPs in gray.Hydrophobic core residues are underlined in the PI-naive subtype B sequence.

Figure 3 .
Figure 3. Conformational sampling of HIV-1 PR in the unbound state.(A) Stack plot of DEER distribution probability profiles.Data are vertically offset for clarity and probability normalized to 100%.The dash line and solid line are placed at 33 Å (closed conformation) and 36 Å (semi-open conformation), respectively.(B) Relative populations for each HIV-1PR construct via deconstruction analysis of the DEER distance profiles, where the flaps are assigned four conformations with the following distances between spin-pairs: curled/tucked (24-30 Å), closed (~33 Å), semi-open (~36 Å), and wide-open (40-45 Å). (C) Difference in population analysis relative to subtype B. Data for B is the same as reported within ref. [76].Full data analysis is provided in supporting information (Figures S1-S7).

Figure 3 .
Figure 3. Conformational sampling of HIV-1 PR in the unbound state.(A) Stack plot of DEER distribution probability profiles.Data are vertically offset for clarity and probability normalized to 100%.The dash line and solid line are placed at 33 Å (closed conformation) and 36 Å (semi-open conformation), respectively.(B) Relative populations for each HIV-1PR construct via deconstruction analysis of the DEER distance profiles, where the flaps are assigned four conformations with the following distances between spin-pairs: curled/tucked (24-30 Å), closed (~33 Å), semi-open (~36 Å), and wide-open (40-45 Å). (C) Difference in population analysis relative to subtype B. Data for Bis the same as reported within ref.[76].Full data analysis is provided in supporting information (FiguresS1-S7).

Figure 4 .
Figure 4. Hydrophobic core residues.(A) Ribbon diagram of subtype B HIV-1 PR with hydrophobic core residues shown in spheres.(B) Sequence of subtypes F and H compared to subtype B with hydrophobic core residues highlighted in gray as well as D60E and I62V mutations shown in blue.

Figure 4 .
Figure 4. Hydrophobic core residues.(A) Ribbon diagram of subtype B HIV-1 PR with hydrophobic core residues shown in spheres.(B) Sequence of subtypes F and H compared to subtype B with hydrophobic core residues highlighted in gray as well as D60E and I62V mutations shown in blue.

Figure 6 .
Figure 6.(A) Ribbon diagram (PDBID 2P3C) in gray showing a side view of HIV-1 PR, highlighting interactions of residues 60 and 62 with the local environment.D60E and I62V are rendered in stick format, whereas other residues are shown as sticks and balls.(B) Overlay of crystal structures of subtype B (1HHP, silver) and subtype F (2P3C, orange) showing D60E to K43 distance.Residues D60 and K43 in subtype B are shown as silver sticks and have an N-O distance of 4.0 Å.Residues E60 and K43 in subtype F are shown as orange sticks and have an N-O distance of 3.0 Å. (C) Interaction of I62 with L38 and Y59 in 1HHP.The residues are shown as sticks, with I62 shown in blue.(D) Probability distribution profile for N-O side chain distances of D/E60 and K43 indicative of salt bridge interactions determined from MD trajectories.

Figure 6 .
Figure 6.(A) Ribbon diagram (PDBID 2P3C) in gray showing a side view of HIV-1 PR, highlighting interactions of residues 60 and 62 with the local environment.D60E and I62V are rendered in stick format, whereas other residues are shown as sticks and balls.(B) Overlay of crystal structures of subtype B (1HHP, silver) and subtype F (2P3C, orange) showing D60E to K43 distance.Residues D60 and K43 in subtype B are shown as silver sticks and have an N-O distance of 4.0 Å.Residues E60 and K43 in subtype F are shown as orange sticks and have an N-O distance of 3.0 Å. (C) Interaction of I62 with L38 and Y59 in 1HHP.The residues are shown as sticks, with I62 shown in blue.(D) Probability distribution profile for N-O side chain distances of D/E60 and K43 indicative of salt bridge interactions determined from MD trajectories.

Figure 7 .
Figure 7. (A) Top view of HIV-1 PR (PDB ID 1HHP) with residues I50, I50′, K55, and K55′ shown as sticks and dash lines showing the distances measured.(B) Flap distance (K55-K55′) and (C) flap tip (I50-I50′) distances as probes for conformation in simulations.Grey coloring indicates distances corresponding to a closed flap conformation, red coloring indicates distances corresponding to the semi-open conformation, blue coloring indicates distances corresponding to a curled/tucked open flap conformation, purple coloring indicates a wide-open flap conformation, and the teal coloring indicates distances consistent with both the wide-open and curled/tucked open flap conformations.

Figure 7 .
Figure 7. (A) Top view of HIV-1 PR (PDB ID 1HHP) with residues I50, I50 ′ , K55, and K55 ′ shown as sticks and dash lines showing the distances measured.(B) Flap distance (K55-K55 ′ ) and (C) flap tip (I50-I50 ′ ) distances as probes for conformation in simulations.Grey coloring indicates distances corresponding to a closed flap conformation, red coloring indicates distances corresponding to the semi-open conformation, blue coloring indicates distances corresponding to a curled/tucked open flap conformation, purple coloring indicates a wide-open flap conformation, and the teal coloring indicates distances consistent with both the wide-open and curled/tucked open flap conformations.

Figure 8 .
Figure 8. Conformational landscape of HIV-1 PR during simulations.(A) Double distance plots representing the conformational landscape of WT, D60E, I62V, and D60E/I62V in MD simulations determined from the distances of each conformer within the simulation.Labels 1 and 2 are assigned to a flap curling/tucking open-like state, with label 3 representing the semi-open conformation, label 4 representing an unknown conformational population and label 5 representing a flap closing-as would be induced by inhibitor or substrate binding.(B) Difference plots between the double distance plots of mutants and WT.

Figure 9 .
Figure 9. (Top row) Front view of the average structure of mutants compared to WT during MD simulations.(Bottom row) Top view of the average structure showing the handedness of WT and mutants.
Figure S2: DEER data processing for the D60E-DRV.
Figure S7: DEER distance profiles and populations with

Figure 9 .
Figure 9. (Top row) Front view of the average structure of mutants compared to WT during MD simulations.(Bottom row) Top view of the average structure showing the handedness of WT and mutants.
Figure S7: DEER distance profiles and populations with DRV.
Figure S8: Distance profiles for flaps and flap tips in a closed conformation.Figure S9.Ribbon diagrams comparing the starting equilibrated structure (1HHP) in yellow and (A) average structures obtained from analysis of conformers sampled during the 1.6 µs simulation runs with (B) focusing on details of the flap tips showing that D60E and I62V flip their handedness during the simulations.Table S1: Analysis of flap and flap-tip distances from HIV-1 PR crystal structured deposited within the PDB.

Table 1 .
Summary of most probable nitroxide distances and fractional occupancy for HIV-1PR constructs.

Table 1 .
Summary of most probable nitroxide distances and fractional occupancy for HIV-1PR constructs.