Owing to the dramatically increased number of 3D structures of pharmaceutical targets, structure based approaches play an important role in drug design. Especially, rapid and accurate estimation of binding affinities is essential for efficient screening of drugs. For this purpose, the free energy based screen methods have been developed by some groups. These methods are much faster than the free energy perturbation (FEP) method or thermodynamic integration method, and give more accurate relative free energies than the empirical scoring function methods.
The λ-dynamics method, one of these methods, has been recently developed by Brooks and co-workers. The λ-dynamics method is an extension of the FEP method. It differs from the FEP method in the following aspects: (1) In the FEP method, a single coupling parameter λ is used to transform one ligand to another. While in the λ-dynamics method, multiple λs, each corresponding to a given ligand, are used. Because of this feature, the binding free energies of multiple ligands are evaluated simultaneously. (2) In the FEP method λ is fixed during the simulation. While in λ-dynamics, the λs evolve according to their equations of motion via use of an extended system. The free energy difference can be obtained from the probability of the ligand i having the dominant state (i.e. λ<SUB>i</SUB>=1, {λ<SUB>m≠i</SUB>=0}). In this method, the total computation time is not expected to increase with the total number of ligands because only the few favorable binders are able to compete for the λ=1 state. This is in coutrast to the conventional free energy calculation methods.
The λ-dynamics method was used to study the binding of 10 heterocycle derivatives to an artificial cavity created inside cytochrome c peroxidase and the host guest system of β- cyclodextrin-monosubstituted benzene derivatives. Straightforward application of λ-dynamics using a multiple topology approach resulted in trapping in local minima. To extend the λ- dynamics method to multiple topology model, a new restraining potential, which keeps the ligands in lower-energy states, is added to the λ-dynamics Hamiltonian. Relatively short time λ- dynamics simulations with the restraining potential successfully identified the best binders as compared with both experimental data and FEP calculations. Using the iterative procedure with biasing potentials to enhance convergence, λ-dynamics method successfully yielded reasonable estimates of the binding affinity of all ligands. Furthermore, long time λ-dynamics simulations revealed that better ligands tend to have small statistical errors, which is appropriate for screening out the plausible ligands from all candidates. This methodology also provides a means to explore the binding orientations and conformations of the ligands inside the binding pocket much better than does conventional MD. A λ-dynamics simulation starting from random initial orientations, in which some ligands take significantly different orientations as compared with those from the X-ray structure successfully sample the X-ray crystallographic orientations in all ligands, even though conventional MD starting from the same initial structures remain trapped in the local minima from which they start. Such an efficient sampling of ligand orientation and conformation is expected to diminish the limitation that an initial ligand structure must be close to its true bound orientation in order to get a reasonable estimate of binding free energy.
The incorporation of the generalized Born (GB) approach into free energy simulation methods (e.g. FEP or λ-dynamics) using thermodynamic cycles was introduced and applied to the trypsin-benzamidine derivatives system and seven mono-substituted benzen derivatives bound to β-cyclodextrin. Since the GB model is fully analytical continuum solvation representation, derivatives of the energy with respect to individual atoms are available and allow the effect of solvation to be efficiently included in molecular dynamics, The GB energies for the intermediate states, in which more than one ligand obtained competitive λ values, were defined by two ways. Both GB coupling schemes have been incorporated into the program CHARMM. The free energy simulations using both definitions of the GB implicit solvent model gave consistent binding free energy differences (ΔΔG) as compared to those using an explicit solvation model. Non- electrostatic solvation energy contributions, which are not included in the GB energies and approximately related to the solvent-accessible surface area, were successfully included using umbrella sampling techniques. Furthermore, a variant of the λ-dynamics approach, Chemical Monte Carlo / molecular dynamics method (CMC/MD), was implemented in CHARMM and compared with FEP and λ-dynamics methods. In the CMC/MD method, the Metropolis Monte Carlo criterion is used to evolve the λ-space and molecular dynamics is used to evolve the atomic coordinates. Free energy differences calculated using the GB energy agreed well among FEP, λ-dynamics, and CMC/MD. The λ dependent partial charge model was also introduced for the incorporation of hybrid topology model into the λ-dynamics method. In this model, the invariable ligand atoms are represented by the single topology and their partial charges are altered according to the movement of the coupling parameters as they would be the same as those of the end points. The hybrid topology λ-dynamics/GB simulations successfully converged without any restraining potential, however, the sampling configurational space was restricted as compared with that of multiple topology model.
They have also applied the λ-dynamics method for the stability analysis of the DNA- binding domain of the Myb transcriptional regulator. In this case, seven different side-chain mutants simultaneously compete to make the protein stabilize, whereas, multiple ligands compete in the previous λ-dynamics simulations. A single short (300 ps) λ-dynamics simulation successfully identified the best stabilized mutant. Furthemore, a series of λ-dynamics trajectories generated by the iterative techniques based on Weighted Histogram Analysis Method (WHAM) showed excellent correlation with data obtained from the conventional FEP simulations or experiments. The additional biasing potentials along λ coordinates successfully increased the ratio of the end points without wasting time by sampling unphysical intermediate states in the λ- dynamics simulations.
Furthermore, a hybrid Monte Carlo and Langevin dynamics method (MC/LD) was introduced to study the binding orientations of toluene in β-cyclodextrin. In this method, the guest atoms are replicated. One of the replicas is sampled with the full force of the host, while the rest of the replicas are sampled with scaled "ghost force" to find the other local minima. In constant time intervals of Langevin dynamics (LD) simulation, Monte Carlo method is applied to choose a new guest among all replicas, which helps to jump the barrier quickly. The trajectory of the MC/LD simulation successfully gave the broader set of binding orientations as compared with the conventional MD trajectory. The potential of mean force calculated using its trajectory gave the good agreement with that calculated by intensive computational use of umbrella sampling method and WHAM. Therefore, MC/LD can succeed in both exploring the free energy surface much more efficiently and yielding the canonical ensemble of the system.
In conclusion, they have developed and applied the free energy based screening methods such as λ-dynamics and CMC/MD. These methods will be used either to rapidly identify ligands with favorable binding free energy or to estimate specific clange in free energy using the iterative procedure with WHAM. Since they screen the binding free energy of the ligands instead of interaction energy, they provide accurate assessment of binding affinity. The restraining potential is very effective and important for the application of the λ-dynamics method to a multiplc topology model. It successfully enhanced both sampling of λ-space and binding configurations of the ligands. The combination of the GB model with λ-dynamics or CMC/MD has a great potential in the application for drug lead optimization. These combinations may fill the gap between the empirical methods using a single minimized complex structure and the theoretically rigorous methods like FEP or thermodynamic integration. The hybrid topology λ-dynamics representation with λ dependent partial charge model will be the promising representation when one investigates free energy changes for an ensemble of slightly varying ligands. Moreover, the MC/LD method can be applied efficiently to explore the free energy surface which will be useful in many purposes such as the protein folding studies, the loop search, or the conformational search of side chains.