The MM/GBSA or MM/PBSA calculations were applied to six different protein systems, including α-thrombin (7 ligands), avidin (7 ligands), cytochrome C peroxidase (18 ligands), neuraminidase (8 ligands), P450cam (12 ligands) and penicillopepsin (7 ligands). The experimental binding data and the PDB entries for the six proteins are listed in Table S1 in the supporting materials. The chemical structures of the ligands are shown in Figure S1 in the supporting materials. The protonated states for all ligands are shown in Figure 1 in the Supporting Materials.
For ligands bound to α-thrombin, cytochrome C peroxidase, neuraminidase and penicillopepsin, MD simulations were performed based on the crystal structures of the complexes. The starting structures of the six avidin analogues (b2–b7) were generated based on the avidin-biotin complex (PDB entry: 1avd33 (link)). The biotin molecule in the crystal structure was manually mutated to the other ligands. It has been shown that the neutral form of the guanidinium group in b2 and b5 biotin analogues is dominant when it is bound to the protein.34 (link) Therefore, the neutral form of the guanidinium group was used in our simulations. The crystal structures of the nine P450cam ligands were used for MD simulations. Starting structures of the other three P450 ligands (e3, e5 and e6) were obtained by manually modifying the ligand (e1) in the crystal structure of 2cpp35 (link) with the conformation of the protein unaltered. The preparation of the models was accomplished in the SYBYL molecular simulation package.36
In the cytochrome C peroxidase complexes, the lone-pair electrons of the epsilon nitrogen in His175 form resonant bonds with the iron ion and the hydrogen atom is located at the delta nitrogen of His175. In the P450cam complexes, lone-pair electrons of the sulfur atom in Cys357 form resonant bonds with the iron ion and this cysteine residue is thus deprotonated. All the crystal water molecules were kept in the simulations.
The atomic partial charges of all ligands were derived by semiempirical AM1 geometry optimization and subsequent single-point Hartree-Fock (HF)/6-31G* calculations of the electrostatic potential, to which the charges were fitted using the RESP technique.37 The reason why we chose AM1 for optimization, not usually used HF/6-31G(d), is to reduce computational cost.38 (link) The optimization and the electrostatic potential calculations were conducted by Gaussian03.39 Partial charges and force field parameters of the inhibitors were generated automatically using the antechamber program in AMBER9.0.40 (link)
In molecular mechanics (MM) minimizations and MD simulations, the AMBER03 force field was used for proteins41 (link) and the general AMBER force field (gaff) was used for ligands.42 (link) The force field parameters developed by Giammona were used for the heme groups in the cytochrome C peroxidase and the P450cam systems.43 To neutralize the systems, counter ions of Cl− or Na+ were placed in grids that had the largest positive or negative Coulombic potential around the protein. The whole system was immersed in a rectangular box of TIP3P water molecules. The water box was extended 9 Å from solute atoms in all three dimensions.