The simulation systems of the HIV-1 PR-inhibitor complexes were prepared using the crystal structures from PDB (Nelfinavir: 1OHR; Amprenavir: 1HPV). The protein complexes were solvated in rectangular TIP3P44 water box using the tleap program of AMBER 10.45 The solute atoms were separated from nearest walls of the water box by at least 12 Å. Sodium chloride counterions were added to make the solution charge neutral. The simulation systems of the solvated inhibitors were prepared in a similar way. The solvated enzyme–inhibitor boxes contain about 34000 atoms. The solvated inhibitor systems contain about 4200 atoms.
The Amber ff99SB parameter set46 (link) was used to model the protein in aqueous solutions. The force field parameters of the ligand were obtained using the Antechamber program47 (link) of AMBER 10. The AM1-BCC method48 was used to assign the atomic charges for the ligands. In the present study, all the MD simulations were performed using the GROMACS program version 4.04.49 The AMBER 10 generated parameters/topology files were converted to the GROMACS format using a Perl script developed by Pande’s group.50 Electrostatic interactions were computed using the particle-mesh Ewald (PME) method,51 with a real space cutoff of 11 Å and a grid spacing of 1.0 Å. A switching function between 9 Å and 10 Å was used for van der Waals interactions. SHAKE52 was used to constrain bond lengths involving hydrogen atoms. A stochastic Langevin dynamics integrator with a friction constant of 0.4 ps−1 and a time step of 2 fs was used to integrate the equations of motion and to provide constant temperature control. The following protocol has been used to minimize and equilibrate the solvated system: the solvent alone was first minimized for 500 steps using steepest descent method followed by 500 steps of conjugated gradient method. Following the minimization steps, the system was equilibrated at the target temperature for 200 ps. The equilibration was performed in the NPT ensemble using Berendsen’s weak coupling method for constant pressure control. Finally, the equilibrated system was simulated in the production MD in the NVT ensemble.
To obtain well equilibrated initial configurations for the decoupling simulations, the solvated enzyme–inhibitor complex and the solvated ligand were equilibrated for 10 and 5 ns, respectively. For ligand–protein decoupling calculation, 11 ns simulation was performed on the solvated enzyme–inhibitor complex at each λ value. The first halves of the trajectory were treated as equilibration and the rest of trajectory were used in the thermodynamic integration. For ligand-water decoupling, 1.2 ns simulation was carried out on the solvated ligands at each λ; the data collected in the last 1 ns were used in the calculation. A total of >3 μs free energy simulation data were collected in the current work. The plots of Hamiltonian derivative 〈∂U/∂λ〉λ obtained from the decoupling simulations were shown in the Supporting Information (SI).
The Amber ff99SB parameter set46 (link) was used to model the protein in aqueous solutions. The force field parameters of the ligand were obtained using the Antechamber program47 (link) of AMBER 10. The AM1-BCC method48 was used to assign the atomic charges for the ligands. In the present study, all the MD simulations were performed using the GROMACS program version 4.04.49 The AMBER 10 generated parameters/topology files were converted to the GROMACS format using a Perl script developed by Pande’s group.50 Electrostatic interactions were computed using the particle-mesh Ewald (PME) method,51 with a real space cutoff of 11 Å and a grid spacing of 1.0 Å. A switching function between 9 Å and 10 Å was used for van der Waals interactions. SHAKE52 was used to constrain bond lengths involving hydrogen atoms. A stochastic Langevin dynamics integrator with a friction constant of 0.4 ps−1 and a time step of 2 fs was used to integrate the equations of motion and to provide constant temperature control. The following protocol has been used to minimize and equilibrate the solvated system: the solvent alone was first minimized for 500 steps using steepest descent method followed by 500 steps of conjugated gradient method. Following the minimization steps, the system was equilibrated at the target temperature for 200 ps. The equilibration was performed in the NPT ensemble using Berendsen’s weak coupling method for constant pressure control. Finally, the equilibrated system was simulated in the production MD in the NVT ensemble.
To obtain well equilibrated initial configurations for the decoupling simulations, the solvated enzyme–inhibitor complex and the solvated ligand were equilibrated for 10 and 5 ns, respectively. For ligand–protein decoupling calculation, 11 ns simulation was performed on the solvated enzyme–inhibitor complex at each λ value. The first halves of the trajectory were treated as equilibration and the rest of trajectory were used in the thermodynamic integration. For ligand-water decoupling, 1.2 ns simulation was carried out on the solvated ligands at each λ; the data collected in the last 1 ns were used in the calculation. A total of >3 μs free energy simulation data were collected in the current work. The plots of Hamiltonian derivative 〈∂U/∂λ〉λ obtained from the decoupling simulations were shown in the Supporting Information (SI).
