Molecular Dynamics and Docking of Biphenyl: A Potential Attachment Inhibitor for HIV-1 gp120 Glycoprotein

Purpose: To develop a new drug that inhibits viral attachment and entry for the treatment of HIV/AIDS patients. Methods: Two Protein Databank (PDB) crystal structures of HIV-1 gp120-CD4 complexes, namely, 1RZK and 1G9N, were mutated at amino acid position 43 to a biphenylalanine (biPhe-43) residue. FireDock web server was used for the docking experiments and 5ns molecular dynamics (MD) using Gromacs 4.0 was performed on the protein complexes to verify the docking results based on the Gibbs free binding energies. Results: Molecular docking by FireDock web server showed that biPhe-43 and Trp-43-mutated CD4 inhibited the binding of gp120 more efficiently, -113.8 and -101.7 kJ/mol (SD = 0, n = 3), respectively, than the alternate aromatic wild type amino acid Phe-43 and the mutant His-43 and Tyr-43. FireDock revealed that electrostatic and Van der Waals interactions were mainly involved in the CD4-gp120 binding and helped to stabilize the protein interactions. In a 5ns MD simulation, biPhe-43 and Trp-43 mutated CD4 demonstrated best Gibbs free binding energies (-16271  29 and -16266 ± 18 kJ/mol, respectively) to gp120 in the identification and confirmation of biPhe-43 and Trp-43 mutated CD4 as excellent inhibitors to gp120. Conclusion: The docked energies and probability outcomes by FireDock anticipated that a ligand for an efficient inhibition of HIV gp120 should involve an extended but conformational flexible aromatic group, i.e. a biphenyl.


INTRODUCTION
Human immunodeficiency virus (HIV) that caused the acquired immunodeficiency syndrome (AIDS) is still the most life-threatening disease which had killed over 25 million people and infected more than 60 million since the beginning of the epidemic [1]. The highly active antiretroviral therapy or HAART has improved the survival and prolonged many lives of HIV/AIDS patients. HAART regimens involved antiretroviral agents from three classes: nucleoside reverse transcriptase inhibitors (NRTIs), non-nucleoside reverse transcriptase inhibitors (NNRTIs), and protease inhibitors (PIs). HAART regimens are effective to suppress the virus but not to eradicate HIV infection, besides the long-term toxicities, multiple side effects and resistance which limit their effectiveness [2,3]. Hence, the hunt for new potential inhibitors that inhibit viral entry is required for HIV victims.
To date, HIV therapy has progressed to the clinical testing of alternative entry inhibitor compounds which consists of three categories: (i) attachment inhibitors, (ii) co-receptor binding inhibitors, and (iii) fusion inhibitors consistent with the beginning early stage of HIV infection [4]. The binding of HIV-1 gp120 to CD4 receptor belongs to the first category which is an attachment process that creates a well-protected cavity and conserved among different HIV strains. In addition, this cavity does not contain glycosylation sites. The CD4 phenylalanine is the only residue binds to this cavity; hence, the cavity of HIV-1 gp120 has been designated as the Phe-43 cavity. It is estimated that CD4-gp120 binding via Phe-43 cavity accounts for a significant 23% of the total binding energy of CD4-gp120 [5,6].
Computational investigations have applied molecular docking and molecular dynamics (MD) simulations to comprehend the binding properties of gp120 and its interaction with the CD4 receptor and co-receptor in the immune cell [7]. The study has suggested that the V3 loop of gp120 is a contact point for cell entry of HIV-1 binding and high occupancy salt bridges may be responsible for transient structural stability in preparation for co-receptor binding. A mathematical model that describes the binding of HIV-1 virus to T cells has also been developed to determine the analytical thresholds for the dosage and dosing interval of HIV fusion inhibitor enfuvirtide [8].
This current study presents an extension using MD simulation to validate the protein-protein docking in a previous study [9] on the binding properties between biphenylalanine-43 (biPhe-43) mutated CD4 and gp120. The molecular dynamics simulation for biphenylalanine was accomplished by the integration of its residue topology parameters in Gromacs simulation software.

Molecular docking of HIV-1 gp120 and mutated biPhe-43 CD4
Two PDB crystal structures namely 1RZK [10] and 1G9N [11] were selected after careful elimination for those containing faulty bonding from a total of 22 PDB crystal structures as described previously [9]. In order to study the interactions of biphenylalanine-43 (biPhe-43) residue on CD4, the Phe-43 cap was mutated using the Chem3D molecular modeling software. The protein complexes were minimized in Gromacs 4.0 [12] using GROMOS96 43a1 force field after the addition of hydrogen atoms by pdb2gmx module; initially with harmonic constraints by steepest-descent and then followed by conjugate gradient without constraints to energy convergence of 0.01 kJ/mol. The complex was separated into individual proteins and subjected to the PatchDock & FireDock [13][14][15] web server directed and blind bound-docking. Directed docking applied the neighboring amino acids Asp190, Glu 192, Asn239 and Trp241 for 1RZK; and Thr108, Asp190, Glu192, Ile193, Asn239, Met240, Trp241 and Gly287 for 1G9N [9], respectively, around the hydrophobic pocket of gp120. Haddock docking web server [16,17] was also employed as a comparison to FireDock. Additionally, molecular docking between wild type Phe-43, mutated His-43 (neutral and protonated), Trp-43 and Tyr-43 CD4 and gp120 was performed to substantiate the reliability of PatchDock and FireDock and Haddock docking web servers.

Docked conformations analysis
FireDock and Haddock output the global binding energy with electrostatic and Van der Waals energies. In addition, FireDock provides additional hydrogen bond and  stacking energies. The effective inhibition of gp120 refers to the more negative global energy and the respective constituent energies. The docked structures were visualized using the Swiss-PdbViewer [18] to qualify the properly docked structures by lock-and-key principle and the probability from a total of 10 docked structures was calculated. The root mean square deviation (RMSD) with respect to the initial minimized bound protein structures was then computed by the Swiss-PdbViewer to provide as criteria of lock-and-key interaction between CD4 and gp120.

Molecular dynamics of mutated biPhe-43 CD4 and gp120 complex
The initially minimized structures of mutated biPhe-43 CD4 and gp120 of 1RZK and 1G9N were soaked in SPC explicit water solvent [19], resulting to a total of 147,077 atoms in a truncated periodic box measuring 10  10  15 nm 3 . Then, the residue topology parameters for biphenylalanine must be incorporated in Gromacs simulation software. As the protein system was not neutral but showed net positive charges, Clˉ ions were added to neutralize the system, so that the PME potential could be calculated without infinite summation. This input structure was minimized in Gromacs 4.0 by steepest-descent with harmonic constraints followed by conjugate gradient without constraints to energy convergence of 0.01 kJ/mol. Molecular dynamics (MD) was performed in Gromacs for 5ns under constant pressure of one atmosphere with a pressure coupling constant of 1.0 ps and at 310 K with a temperature coupling constant of 0.1 ps. The MD time step was set to 0.001ps (1fs). The nonbonded pair list was updated every 10 steps. A 1.0 nm cutoff was applied to Van der Waals and short-range Coulomb interactions. The output was updated for every ps and frames were saved for every 10 ps. Each 5ns of 'bound' MD needs 30 days computational time and seven days for 'unbound' MD on UMHPC Linux Cluster SGIAltix serial processor. The g_lie module in Gromacs was then used to calculate the Gibbs free energy of binding (G bind ) of Aqvist [20] as shown in Eq 1. The approach is based on the assumption that a receptor-ligand protein-protein interaction can be estimated by determining the difference in solvent (water) interaction energy of the free ligand (E unbound ) and the complex (E bound ). The remaining variable, i.e. the receptor and its water interaction, can be ignored, as it remains constant in the comparison, thus only affects the total energy, but not the difference between the ligands. In the g_lie module, the long-range Lennard-Jones (LJ-14, E VdW =Elj) and Coulomb (C-14, E Coul =Eqq) potentials between the ligand and solvent must be provided in the -Elj and -Eqq options. The root mean square deviation (RMSD), Van der Waals and electrostatic energies were calculated by Gromacs modules g_rms and g_energy, respectively, over 5ns MD trajectories. The MD snapshots and hydrogen bonding after 5ns simulation were displayed by the PyMOL Molecular Graphics System, Version 1.1 Schrödinger, LLC.

Statistical analysis
The docking and simulation data were expressed in mean  standard deviations (SD). The docking experiments by FireDock web server were performed in three replications and the mean and SD were calculated by Microsoft Excel 2007. The mean and SD for molecular dynamics were computed by the g_energy module in Gromacs 4.0 over 5ns MD trajectories.

Docking outputs and the inferences
The docked energies for Haddock and FireDock are shown in Table 1. Haddock web server indicated best binding energy for 1RZK Trp-43 but least binding for 1G9N Trp-43 and showed contradictory docking energies among 1RZK and 1G9N mutants. Besides, Haddock could not compute the docking of biPhe-43 because the topology parameters were not available in the database. On the other hand, FireDock consistently showed greatest binding energies for biPhe-43 and Trp-43 for both 1RZK and 1G9N by directed docking. However, as expected, the blind docking results showed a slight discrepancy for biPhe-43 and Trp-43 for 1RZK and 1G9N. Hence, the best scheme is using FireDock directed docking opposed to blind docking. As shown in the parenthesized fractions in Table 1, the docked probability favoured smaller amino acids, i.e. Phe-43 and His-43 but modest docked probability for biPhe-43 due to its structural flexibility whereas Trp-43 showed moderate to lower docked probability attributable to its rigidity. Table 2 shows that the root mean square deviation (RMSD) of around 0.1 to 0.5 Å for both 1RZK and 1G9N wild type and mutants for properly docked conformations that were in conformity to lock-and-key theory. The RMSD results proposed a measure to qualify proper docked conformation as RMSD  1.0 Å. The energy contribution to the global energies in Table 1 was mainly caused by electrostatic and Van der Waals energies (Table 3). Table 4 shows the Gibbs free binding energy computed from 5ns of molecular dynamics (MD) simulation. MD results showed that best binding energy was demonstrated by biPhe-43 and Trp-43. These outcomes were in a comparable inclination with the global energies from directed docking by FireDock ( Table 1). The root mean square deviation (RMSD) for biPhe-43 1RZK and biPhe-43 1G9N solvated protein systems over 5ns of MD is shown in Fig 1. In general, the RMSD for biPhe-43 1RZK was fluctuating slightly but has reached convergence towards structural stability for biPhe-43 1G9N. The RMSD for 1RZK

DISCUSSION
As shown in Table 1, the docking results were in good agreement with an enzyme linked immunosorbent assay (ELISA) study signifying enhanced binding affinity of gp120 for a biphenylalanine (biPhe-43) mutated CD4 [21]. The data in Table 3 enclose the same opinion with a review on mechanisms of action and resistance pathways of HIV entry inhibitors [5] which indicated that electrostatic forces mainly drive the CD4-gp120 binding and the Van der Waals forces help to stabilize the CD4-gp120 interactions. The hydrogen bonding and - interactions (limited to Trp-43) were inferior in stabilizing the CD4-gp120 binding.
The paramount electrostatic and Van der Waals interactions found in biPhe-43 and Trp-43 (Table  3) were in complementary proportion to the global energy for 1RZK and 1G9N as shown in Table 1. The hydrogen bonding did not diverge radically along with the mutated amino acids except for Tyr-43. This was owing to the nonpolar physicochemical property of Tyr-43. MD was proficient to validate the FireDock docking results to infer biPhe-43 and Trp-43 as efficient inhibitor to gp120 in both 1RZK and 1G9N protein models.
The Gibbs free binding energies of 1G9N Phe-43, His-43, Trp-43 and Tyr-43 were not included on account of inconclusive results. Our MD RMSD data were within the range to the protein dynamics of 4-9 Å RMSD investigated for the structure and correlated motions between the contact points for cell entry of HIV-1 gp120 to CD4 co-receptors [7]. The fluctuation in 1RZK was more drastic than 1G9N indicating that 1G9N was more stable throughout the simulation.
The present study serves as a motivation of pursue for an attachment inhibitor against HIV-1 gp120 glycoprotein to complement the currently  available entry inhibitors such as co-receptor binding and fusion inhibitors [4]; and the traditional antiretroviral HAART regimens. The benefit of exploration headed for a lead structure by computational molecular modeling and docking is that it is a truthful and precise rational drug design approach.
The present study serves as a motivation of pursue for an attachment inhibitor against HIV-1 gp120 glycoprotein to complement the currently available entry inhibitors such as co-receptor binding and fusion inhibitors [4]; and the traditional antiretroviral HAART regimens. The benefit of exploration headed for a lead structure by computational molecular modeling and docking is that it is a truthful and precise rational drug design approach. This is because the mechanism and insights into the target binding sites have been revealed unambiguously by means of protein X-ray crystallography technique. However, the future design of HIV-1 attachment inhibitor should be a chemical compound or a small ligand which can be orally prescribed. This is because the presently available fusion inhibitor such as enfuvirtide (fuzeon/T-20) targeting HIV-1 gp41 glycoprotein [22,1] is an extremely expensive 36-mer synthetic oligopeptide and must be administered by subcutaneous injection. Hence, the frequent occurrence of painful injection site reactions has limited the long-term application [23].

CONCLUSION
The search for a lead attachment inhibitor is a pressing need in addition to the presently available anti-HIV/AIDS regimens. Molecular docking and molecular dynamics have shown harmonized and promising results in the identification and verification of biPhe-43 and Trp-43 as the lead candidates en route for the design of new attachment inhibitors against HIV-1 gp120 affection glycoprotein to human Tlymphocyte CD4 receptor. A ligand for an efficient inhibition of HIV-1 gp120 should engage a flexible aromatic group, i.e. a biphenyl.