Recent Advances in Protein Ligand Interactions

Send Orders for Reprints to [email protected] Current Computer-Aided Drug Design, 2013, 9, 000-000

1

Recent Advances in Protein
Ligand Interactions:
olecular Dynamics Simulations and Binding Free Energy Kshatresh Dutta Dubey, Rakesh Kumar Tiwari and Rajendra Prasad Ojha* Biophysics Unit, Department of Physics, DDU Gorakhpur University, India Abstract: Computational techniques are one of the most emerging topics in structural and molecular biology. Molecular dynamics (MD) simulations are used not only to explore the conformational aspects of biological systems but also to have significant scope in protein
ligand interactions. Then the binding free energy calculations are readily applied to the simulated systems in order to predict the binding affinities. The thermodynamic properties are directly related to protein
ligand interactions which are dependent upon a few specific parameters. In the present review, we highlight some facts related to protein
ligand complexes, by starting with a survey of MD simulations and binding free energy calculations and ending with some successful implementations of these computational techniques.

Keywords: FEP, LIE, ligand binding, MD simulations, MM-PBSA, QM/MM, thermodynamical parameter. 1. INTRODUCTION The accurate determination of absolute binding affinities for complex systems is an issue of central importance in computational chemistry and biology. As to protein
ligand complexes, the binding free energy is crucial to identify new candidates that potentially bind to target receptors [1]. Although the microscopic bimolecular interactions have been relatively well understood, the computational schemes towards the accurate binding free energies remain a challenge. A number of schemes have been proposed at various levels of sophistication; for example, Molecular Mechanics/Poisson
Boltzmann Surface Area (MM/PBSA) and Molecular Mechanics/Generalized Born Surface Area (MM/GBSA) [2-4] that have been used extensively, combine molecular mechanics and implicit solvent models, which have been reviewed elsewhere [5-8]. The linear interaction energy (LIE) method [9] is a similar approach that averages interaction energies from molecular dynamics (MD) simulations to estimate the absolute binding free energies. In contrast to the rigorous methods such as Free Energy Perturbation (FEP) and Thermodynamic Integration (TI) [10], MM/PB(GB)SA and LIE are computationally cheap. The main limitations of FEP and TI methods are the slow convergence and the exhaustive conformational sampling in order to obtain a proper averaging ensemble. In this review, we aim to address the recent developments in MD simulations and binding free energy calculations as well as their applications on protein
ligand complexes. We also discuss the future aspects of MD simulations and binding free energies.

time step (often at the order of fs), where forces on each atom are computed at every time step by using force fields. The force fields are usually potential energy functions including bonded and non
bonded potential terms. The velocity and coordinates of the systems using Newton's law of motion are regularly updated during the MD simulations and the trajectory is obtained meanwhile. The first MD simulations for a biological macromolecule were carried out by Karplus and coworkers 35 years ago [11]. A very short simulation time (9.2 ps) is used for bovine pancreatic trypsin inhibitors. The experimental observations such as hydrogen bonding exchange have already been explained before the use of MD simulations [12]. Two years after the first MD simulations, the role of thermal (B) factor in internal motion of proteins that was calculated during X
ray crystallographic refinement, is investigated [13, 14]. As a consequence of these studies, the estimation of mean square fluctuations versus residue number becomes an ingredient of almost all the MD analyses. Subsequently, a wide range of motional phenomena for proteins and nucleic acids are subjected to MD simulations, focusing on the conformational flexibility and interpretation of experimental results, such as the analysis of fluorescence depolarization [15], dynamics of NMR parameters [16] and effect of solvent and temperature on protein stability [17, 18]. During these time periods, simulated annealing is widely used for x-ray structure refinement [19] and NMR structure determinations [20]. Meanwhile, several applications are found for demonstrating the internal motion in biological functions such as the hinge bending mode for opening and closing the binding sites [21]. 3. IMPORTANCE OF MD SIMULATIONS

2. BEGINNING OF MD SIMULATIONS In MD simulations, the atoms of the macromolecular systems and surrounding water molecules move with a short *Address correspondence to this author at the Biophysics Unit, Department of Physics, DDU Gorakhpur University, India; Tel: +91-551-2202167/ 2332398; Fax: +91-551-2330767/2340459; E-mail: [email protected] 1573-4099/13 $58.00+.00

The MD simulations, which provide a methodology for detailed microscopical modeling at an atomistic scale, are a powerful and widely used tool in chemistry, physics, and material science. This technique is a scheme of studying the natural time evolution for molecular systems that allows prediction of the static and dynamic properties directly from the underlying interactions. Dynamic simulations monitor © 2013 Bentham Science Publishers

2 Current Computer-Aided Drug Design, 2013, Vol. 9, No. 4

the time
dependent processes in molecular systems in order to study their structural, dynamic, and thermodynamic properties by numerically solving an equation of motion, which is the formulation of the rules that govern the motion. That is, MD simulations provide information about the time dependence and magnitude of fluctuations in both position and velocity. Another significant aspect is that, although the potentials used in the MD simulations are approximate, they are completely under the control of users, and by changing or removing some specific circumstances their roles and influences can be examined. The application of MD simulations can be roughly classified into three types [22]: 1) MD simulations are used mainly for the conformational sampling and refinement or determination of the results obtained from X
ray, NMR and other experimental techniques. 2) MD simulations are used to describe the properties of equilibrated systems. In this case, the thermodynamic property, root mean square deviation, thermal fluctuation, motion of center of mass and correlation factor are generally elaborated. These can be assumed as a second stage of MD applications because the systems of interest should have already been correctly refined and the conformational sampling should have been properly carried out. 3) The motion and evolution with the simulation time are of primary interest. An accurate development of the systems is a prerequisite, besides the sufficient sampling of the configurational space. 4. PRESENT AND FUTURE OF MD SIMULATIONS Currently, the MD simulations are being used for investigation of nearly all types of macromolecules. The number of MD publications per year reaches thousands. Therefore, it is not possible to make a complete review here, and so we will focus on our own works and other related contributions. We would like to start with the revolutionary contribution from M. Karplus who is the first to use computational MD simulations. Berendsen [23] and Go [24] also make significant contributions by extending MD methods and developing normal mode dynamics, respectively. To run the MD simulations on the time scales for chemical processes, a number of algorithms have been developed to maintain the system temperature, including velocity rescaling like Berendsen [23], temperature relaxation like Nose
Hoover thermostats [25] and velocity modification such as Andersen [26] and Langevin thermostats [27]. For the Anderson thermocoupling, the particle velocities are periodically reassigned to the pseudo
random values so that the resulting momenta follow a Maxwell distribution around the desired temperature. For Langevin, velocities of the particles during the MD simulations are modified with pseudo
random forces as if they are undergoing stochastic collisions with imaginary particles whose momenta follow a Maxwell distribution at the desired temperature. A good review of thermostats can be found elsewhere [28]. The rescaling of velocities may result in a serious problem, often known as 'flying ice cube effect’. It is an artifact where simulated systems acquire high linear momentum and show extremely damped internal motion, thus freezing the conformational flexibility and causing the

Dubey et al.

system to be rigid as the ice cube floating in space [29]. This artifact is a consequence of the MD algorithms. The ice cube effect arises due to the repeated rescaling of the velocities of particles in simulated systems. The problem can be largely avoided by periodically removing the motion of center of mass and alternatively using the Nose
Hoover thermostats. Recently, a serious problem associated with the implementation of Andersen and Langevin dynamics algorithm has been reported [30]. It is possible for a repeating sequence of pseudo
random numbers to enter into the long MD simulations, which can quickly denature the biomolecules and can be easily detectable if the sequence repeats within a short time interval. With advance of the computer power, we are able to observe the Central Dogma of molecular biology — DNA to RNA to protein, using MD simulations. Atomistic MD simulations have already demonstrated the de novo folding of small (up to 80 residues) protein domains. The effect of force fields on the folding quality has been mentioned [31], and the continuous force field improvement enables further progress. It is expected that the folding of more typical single
domain proteins (~ 300 residues, resulting in ~ 105 atoms including explicit water molecules with the simulation time at the millisecond order) will be viable within the next 10 years, and the folding of large, multi
domain proteins — for instance,

galactosidase (1024 residues), comprised of four strands with each of five domain subunits — is likely within the 25 years. Besides the delivery of folded structures, the MD simulations help us to understand the basic folding mechanisms. Through the folding studies of numerous proteins, these MD simulations will provide the data for developing a model at a conceptual level that describes the general and unambiguous features of protein folding mechanisms [32], and folding simulations may not be needed then. With an ever increase of potential drug candidates (almost one billion drug
like compounds), the MD simulations in drug discovery have an utmost importance. As MD simulations enable us to know the folding mechanisms, they can be implemented for structural determination of target molecules whose experimental structures are not available to us. The MD simulations can be combined with other computational chemistry tools, to identify novel drug binding sites; The MD simulations can also be carried out on the target molecules alone or with presence of suitable candidate fragments. Moreover, results obtained from subsequent binding free energy calculations can be used for the scoring of candidate fragments. 5. PROTEIN
LIGAND BINDING Discovery and development of a new efficient drug is definitely a long
time and expensive process because the optimal potency and minimal side effect should be considered at the same time. With advancement of high
throughout screening and combinatorial techniques, the data of bioactivity and structures of drugs increase dramatically, and analysis of such huge database becomes especially difficult. A new drug has to go through many filtering criterions such as toxicity, selectivity and binding affinity. The explosion of database needs urgently the integration of chemical information with molecular modeling

Recent Advances in Protein-Ligand Interactions

techniques. After the discovery of a new drug, a series of related lead compounds are expected to be explored. A lead series comprises a lot of compounds that usually share some common features, and shows some activity variation when the structure is modified. Discovering a lead series by experimental techniques is apparently more complicated than that of a single molecule. Drug
receptor binding is one of the most active fields for biochemists and bioinformatics. Receptors are the macromolecules involved in chemical signaling between and within cells. Molecules (drugs, DNA, peptides etc.) that bind to receptors are called ligands. Although targeted for a receptor, most drugs have selectivity for different receptors. Selectivity is the degree to which a drug acts on a given site relative to other sites. 6. MAJOR CONTRIBUTIONS OF PROTEIN
LIGAND INTERACTIONS The selective binding of a low
molecular
weight ligand to a protein is determined by the structural and energetic recognition, and the binding affinity can be derived from the experimentally measured binding constant Ki,
G =
RTlnKi =
H
T
S The experimentally determined binding constant Ki is typically in the range 10-2
10-12 M, corresponding to a Gibbs free energy of binding
G between -10 and -70 kJ/mol in aqueous solution [33,34]. There is now a large body of experimental data available on the three
dimensional (3
D) structures and binding affinities of protein
ligand complexes, which clearly indicate that there are several features in all complexes of tightly binding ligands: 1.

There is a high level of steric complementarities between protein and ligand. This is lively described as a paradigm of lock
and
key. However, it is now considered that lock
and
key idea is obsolete and is not appropriate to describe most biological complexes. The protein
ligand interactions are more like a hand
and
glove association, where both of protein and ligand show flexible and adjust themselves to complement each other. The mechanism is known as an induced fit where protein and ligand modify their shape and mould their complementarities so as to increase favorable contacts and reduce adverse interactions, thus maximizing the binding free energy.

2.

There are usually high complementarities of the surface properties between protein and ligand. The lipophilic parts of ligand are frequently found to be in contact with the lipophilic parts of protein. The polar groups are usually paired with suitable polar protein groups to form hydrogen bonds or ionic interactions. The experimentally determined hydrogen bond geometries display small scatter; in other words, the hydrogen bond geometry is strongly preserved. With few exceptions, there are no repulsive interactions between ligand and protein.

Generally, direct interactions between protein and ligand are important for the binding. Hydrogen bonds, ionic, hydrophobic and cation
pi interactions are examples of such

Current Computer-Aided Drug Design, 2013, Vol. 9, No. 4

3

direct interactions. Structural data on unfavorable protein
ligand interactions are sparse, partially because structures of weakly bound ligands are more difficult to obtain and usually considered less interesting by structural biologists. However, these data are vital for the development of scoring functions. Some clues can be drawn from these data: unpaired buried polar groups at the protein
ligand interface are strongly averse to the binding. Few buried CO and NH groups in the folded proteins fail to form hydrogen bonds [35]. Therefore, in the ligand design process we have to ensure that the polar functional groups, from either protein or ligand, will find suitable counterparts if to be buried upon the ligand binding. Another situation that leads to a decrease to the binding affinity is an imperfect steric fit, leading to holes at the lipophilic part at the protein
ligand interface. The enthalpic and entropic components of the binding affinity can be determined experimentally; e.g., by isothermal titration calorimetry (ITC). Unfortunately, these data are still sparse and are difficult to interpret [36]. The data available indicate that there is compensation between enthalpic and entropic contributions [37-39]. The binding can be enthalpy
driven (e.g., streptavidin
biotin, G = -76.5 kJ/mol, H= -134 kJ/mol) or entropy
driven (e.g., streptavidin
HABA, G = -22.0 kJ/mol, H=7.1 kJ/mol) [40]. 6.1. Hydrogen Bonds in Protein-Ligand Interactions Hydrogen bonding plays a significant role in a variety of chemical and biological processes, including ligand binding and enzyme catalysis [41, 42]. Consideration of hydrogen
bonding properties in drug design is important because of their strong influence on binding specificity, transport, distribution, metabolization, and excretion of ligands. Their ubiquity and flexibility make hydrogen bonds one of the most important physical interactions in biological systems [43]. Because hydrogen atoms comprise approximately half of the atoms of biological macromolecules and two
thirds of the atoms of the solvating water molecules, hydrogen atoms are found among almost every pair of non
covalent bonded heavy atoms. The necessary condition of a hydrogen bonding is that a proton lies between the electron clouds of two atoms and modifies their interaction in a manner that is not explicable in terms of the van der Waals (dispersion
repulsion) effect. Hydrogen bonding competes with van der Waals interactions in number. van der Waals interactions occur unavoidably and with similar strengths between all atoms, their contribution to selectivity of interactions largely lies in the shape selection caused by the repulsive component of these interactions. Consequently, from both evolutionary and design perspectives, modification of local hydrogen
bonding potential is the principal mechanism available for favorably enhancing the interactions between molecule pairs. The popular notion of “hydrophobic” or “lipophilic” forces being important is merely a result of a non
atomic perspective. The hydrophobic forces, while being a simplifying concept, are a complex phenomenon resulting from redistribution and change in the strength of water
water hydrogen bonds as solvent is released upon approach of two apolar chemical groups. This may seem to be a somewhat enthusiastic view of the importance of hydrogen bonding, and yet it is a consequence of regarding all forces from an atomic

4 Current Computer-Aided Drug Design, 2013, Vol. 9, No. 4

perspective and extrapolates a developing trend in the past few decades, which has seen an increasing diversity of hydrogen
mediated interactions that are considered significant in modulation of biological systems. An example of hydrogen bonding in protein
ligand interactions including solvent is shown in Fig. (1).

Dubey et al.

thermodynamic properties and relate them to various structural and functional characteristics of protein and found that the high
energy hydration sites often exist near the protein motifs typically characterized as hydrophilic, such as backbone amide group. In addition, the distribution of high
energy hydration sites on the protein surface can be used to identify the location of binding sites and that binding sites of druggable targets tend to have a larger density of thermodynamically unstable hydration site [46]. Some similar studies on the role of bound water molecules can be found elsewhere [47-54].

Fig. (1). The surface representation of a protein
ligand complex. The ligand in the binding pocket is represented by the ball and stick model. The colour of surface ranges from dodger blue for the most polar residues to white and from orange to orange red for most hydrophobic residues (for the interpretation of references to colors in this figure, the reader is referred to the web version of this article).

6.2. Role of Solvent in Protein
Ligand Interactions The binding of protein
ligand occurs in a salt water environment, which strongly influence the energetics of the systems. Because water has a dielectric constant of about 80 whereas the dielectric constant of vacuum is 1, it leads to a new one
body solvation energy for each atomic charge that arises from the favourable interactions between the charges and high dielectric environment [44]. In addition, water screens the charge
charge interactions of fully hydrated atoms by approximately 80 fold. However, atoms in a protein
ligand interface are sequestered from the solvent and therefore interact with an effective dielectric constant less than 80. Water gives rise to the hydrophobic effect, the tendency of water molecules to drive nonpolar solute together. This promotes the association of nonpolar surfaces of the ligand and the protein. The hydrophobic effect is often accounted for by an additional solvation energy term that is proportional to molecular surface area, with a positive coefficient. The effect is to add positive (unfavorable) solvation energy to conformations with more surface area and thus favour binding, which reduces the surface area. An example of the binding of solvent with protein is shown in Fig. (2). In addition to the hydrophobic effect and other bulk properties of solvent, water also plays an important role in interactions with protein
ligand complexes. A recent study on 28 low micro-molar inhibitors reveals that the implementation of statistical thermodynamic analysis enables a successful explanation of some complex trends in structure
activity relationship (SAR) [45]. Another study indicates that water also plays an essential role in determining the structure and function of many biological systems. In this work, the authors have characterized some

Fig. (2). The interaction of receptor
ligand and solvent
receptor through hydrogen bonds represented by spiral springs (for the interpretation of references to colors in this figure, the reader is referred to the web version of this article).

6.3. Effect of pH on Protein
Ligand Bonding The change in pH plays a crucial role in the binding of protein
ligand complexes [55]. The change in pH controls the electric charges of acidic and basics groups of both protein and ligand. When protein binds ligand, hydrogen ions may be liberated or consumed, and the apparent equilibrium constant K for ligand dissociation is a function of pH. From a thermodynamic point of view, the acid titration curves of protein, ligand, and protein
ligand complex are related to the pH dependent K through the hydrogen ion binding polynomials P (partition functions) of the three constituents. However, the hydrogen ion binding polynomials of protein and protein
ligand complex can be factored into a binding polynomial for the binding site and a binding polynomial for the rest of protein. When ligand does not have acid dissociations in the range of pH where the unoccupied site or occupied site is dissociated, the pH dependent K is given as the ratio of the binding polynomial at the unoccupied site to that at the occupied site [56]. A conformational study on Antheraea polyphemus PBP1 reveals that recombinant Antheraea polyphemus PBP1 (ApolPBP1) picks up hydrophobic molecule(s) endogenous to the Escherichia coli expression host that keeps the protein in the "open" (bound) conformation at high pH but switches to the "closed" (free) conformation at low pH. The study of several protein
ligand binding shows that the change in pH triggers the open state of proteins to the closed state and thus significantly affects the protein
ligand binding [57]. By development of the constant pH MD simulations, the role pH

Recent Advances in Protein-Ligand Interactions

for the system has become computationally efficient. Unlike the classical MD simulations, where protonation states of all titratable groups are kept fixed, the constant pH MD simulations allow the change of protonation states for the titratable groups during the MD simulations. The sampling of protonation states in protein is usually implemented by use of the PB [58-60] or GB model [61], combined with the Monte Carlo sampling. The constant pH MD simulations account for the aqueous environment in a realistic way and therefore their main contribution is likely to lie in the pH
dependent conformational change and stability [62]. The accuracy of the constant pH MD simulations has been assessed by reproducing the experimental pKa values, and these results are generally encouraging [63, 64]. 6.4. Entropic Effect The conformational entropy is an important component of the free energy change upon ligand binding to its target protein. Consequently, the development of computational techniques for reliable estimation of conformational entropy currently receives an increased level of attention in the context of computational drug design. Detailed knowledge of the thermodynamics of ligand binding is crucial for design of potential drugs. In this context, the change of Gibbs free energy,
Gbind, upon binding of ligand to its target determines the ligand binding affinity, which is usually directly related to its efficacy as a biologically active compound. Taking into account that
Gbind is a sum of two terms (
Gbind =
Hbind
T
Sbind), ligand design can benefit from the optimization of both enthalpic (
Hbind) and entropic (T
Sbind) contributions [65]. However, in practical implementations, protein
ligand interactions are usually optimized through the enthalpic component, while entropy is primarily attributed to the hydrophobic effect and solvation, which are believed to be the major contributors [66, 67]. Nonetheless, entropic contributions upon binding also include the changes of conformational entropies from both protein and ligand; i.e.,
Sbind =
SconfP +
SconfL +
Ssol +
SRTP +
SRTL, where
SconfP and
SconfL are the changes in conformational entropies of protein and ligand, respectively,
Ssol is the change of solvent entropy, while
SRTP and
SRTL correspond to entropy changes due to rotational and translational degrees of freedom for protein and ligand, respectively. Although the potential importance of conformational entropy has been emphasized several decades before [68-70], experimental evidence has been provided only recently indicating that it can dramatically influence the binding free energy of protein
ligand association [71-74]. For example, in studies of calmodulin interactions with fragments, Wand and coworkers [71, 74] show that
SconfP and
Ssol can, in fact, be of the comparable magnitude. Similarly, the consideration of changes in ligand conformational entropy (
SconfL) upon binding has been found to be important in design of conformationally restricted binders [75-77]. A combination of MD simulations with NMR experiments is a particularly powerful approach in addressing the challenges inherent in conformational entropy estimation [78
82]. Traditionally, modeling of protein flexibility has been considered in drug design mainly in the context of improving standard docking techniques for

Current Computer-Aided Drug Design, 2013, Vol. 9, No. 4

5

predicting and optimizing properties of protein
ligand complexes [83-85]. However, different efforts have recently been extended in the direction of including different explicit representations of conformational entropy even in docking algorithms [86-88]. What is more, the computationally more demanding, MD
based approaches for entropy estimation are being developed in ligand design tasks [89, 90]. We are witnesses that this trend will keep growing and that the majority of all future computational drug design efforts, regardless of their level of complexity, will take into account the conformational entropy. There are several computational techniques that provide conformational entropy estimations based on the MD simulations [91-94]. The Cluster variation method (CVM) is among such methods where Monte Carlo probability is used without thermodynamic integration (TI) [95]. Quasi
harmonic approximation is another long
standing approach, considering the probability density function (PDF) of internal coordinates as a multivariate Gaussian distribution [96, 97]; however, this simple Gaussian model does not yield accurate entropy estimations when the true multivariate distribution of the spatial fluctuations differs significantly from the Gaussian distribution, especially when the PDF is multi
modal [98, 99]. The nonparametric nearest
neighbor (NN) method of Hnizdo et al. [99], and mutual
information expansion (MIE) by Killian et al. [100], are two other methods for the calculation of configurational entropy. The combination of these methods is also proposed by the same group to furnish an efficient and accurate extraction of configurational entropy from molecular simulations to a given order of correlations among the internal degrees of freedom [101]. Recently, two novel methods using minimally coupled subspace approach (MCSA) [102] and using statistical properties of the solvent
averaged protein potential energy surface [103] have been developed. The first overcomes the limitations such as the quasi
harmonic approximation which neglects the nonlinear and higher
order correlations as well as the multi
minima characteristics of protein energy landscapes, whereas the second is implemented by combining molecular simulations and integral
equation theory of liquids. This method does not assume the Gaussian distribution of protein configurations and can be applied to unfolded or misfolded states in which an average protein structure has not been well defined. Although most of them are still out of the arena of practical drug design, a golden standard for developing novel computationally less demanding methods for such applications has been provided, in addition to furthering our understanding to protein-ligand interactions [73, 78-84]. 6.5. Effect of Protein Flexibility on Ligand Interactions Proteins are in constant motion among various conformational states with similar energies and these motions are often ignored in drug design. However, protein flexibility is fundamental to understand the ways in which drug exerts the biological effects, binding site location, binding orientations etc., thus allowing an increased affinity between the drug and its receptor. A review on protein flexibility and its implementations on drug discovery can be found elsewhere [104]. A variety of biological processes are based on ligand
receptor interactions. When ligand binds to

6 Current Computer-Aided Drug Design, 2013, Vol. 9, No. 4

Dubey et al.

receptor, the optimal ligand-solvent and receptor-solvent interactions are replaced by optimal intermolecular interactions among ligand, receptor and solvent. In this process the flexibility of receptor plays an important role because configurational changes in the receptor may give rise to enthalpic and entropic contributions to the binding free enthalpy and/or Gibbs energy [105, 106]. The degree of these contributions depends on the flexibility of the receptor and is generally difficult to be evaluated [107]. It has been illustrated that in computer
aided studies of biomolecular complexation a receptor can be represented as a rigid body if no structural rearrangement is necessary to accommodate different ligands [108]. In the case of flexible receptors such as DNA where ligand and receptor are in close contact with each other, neglect of the receptor flexibility affects the relative free enthalpies of binding as well as the structures of complexes. The differences in the free enthalpy of ligand binding to flexible versus rigid receptors do not necessarily arise only from non
optimal ligand
receptor interactions, but also ligand
solvent interactions and the loss of configurational entropy upon restraint may also contribute. 7. COMPUTATIONAL APPROACHES PROTEIN
LIGAND BINDING

FOR

where the superscript u indicates the unbound conformation of ligand Unbound protein free energy: The protein is treated by classical mechanics during the simulation for unbound protein. Its free energy is composed of total intramolecular energy Egas [(Evdw) + (Ecoul) + (Ebond)], solvation free energy (Esolv) and entropic contributions (TSideal), G (P) = Egas( Pu) + Gsolv(Pu)
TSideal(Pu) Complex free energy: It is the sum of energies corresponding to the protein and ligand in the bound state. The free energy of complex is composed into the following terms, G(C) = Egas (Pb) + Egas(Lb) + Gsolv(Cb)
TSideal Here Egas is decomposed to Ebond, Ecoul and Evdw with superscript ‘b’ denoting the bound form. The binding free energy for the non
covalent interactions of two associated molecules can be written as, Gbind = G (L+P C) = G(C) - G (L) - G (P) or Gbind = Eint + Gsolv
T
Sideal Here, E = Eint( Pb)
Eint( Pu); int

In the previous sections, some important issues are discussed for protein
ligand interactions. There are several computational techniques, for example, molecular docking, 2
D and 3
D QSAR, MD simulations, MM/PBSA, LIE and FEP that are widely used by computational biologists. Here we focus on the MM/PBSA, LIE and FEP methods. 7.1. MM/PBSA Methods The MM/PBSA method along with its GB variant (MM/GBSA) uses MD simulations of the free ligand, free protein, and their complex as the basis for calculating their average potential and solvation free energy. In this approach G is written in terms of gas phase contribution, energy difference due to translational and rotational motion, desolvation energy and entropic contributions, G = Ggas + Gtran/rot + Gsol
TSideal The first term is gas phase contribution which is the sum of the electrostatic energy, van der Waals energy and internal energy; i.e. Ggas = Evdw + Eele + Eint. The second term is the energy due to translational/rotational motion, and in classical mechanics equals 3RT whereas is generally omitted in MM/PBSA. The third term is the solvation energy which is further composed into nonpolar solvation and polar solvation energies. The last term is the entropic contribution which depends upon the degree of freedom of the translational, rotational and vibrational motion of the system. Theoretically, the above formula is applicable for the ligand, protein and complex. The resultant binding free energy is expressed as [109], G = Gcomp
(Glig+Gprotein) Unbound ligand free energy: The free energy for the ligand is given by G(L) = Egas(Lu) + Gsol(Lu)
TSideal(Lu)

Gsolv = Gsolv(Cb)
Gsol(Lu)
Gsolv(Pu), and
TSideal =
TSideal + TSideal(Lu) + TSideal(Pu) where Eint is the change in protein intramolecular energy. Gsolv [110, 111] is composed of nonpolar and polar contributions. Nonploar solvation energy accounts for the unfavorable cavity formation and favorable van der Waals interaction between solute atoms and solvent [110, 111], Esolv,np = A + b where A stands for solvent accessible surface area (SASA), whereas  and b are empirical constants which may have different values. The polar solvation energy is calculated using the linear Poisson
Boltzmann (PB) equation [112] that relates the charge density, (r) to the electrostatic potential, (r), in a medium with non
uniform dielectric permittivity (r) that equals 1 and 80 for the solute and the solvent, respectively,  (r)  (r) =
4 (r) + 2 (r)  (r) where  is the Debye
Huckel screening parameter taking into account the electrostatic screening effect. TS arises from the change in translational, rotational and vibrational degrees of freedom upon ligand binding. A protocol for free energy calculation is shown in Fig. (3). The entropy term, due to the loss of degrees of freedom upon association, is decomposed into translational Strans, rotational Srot, and vibrational Svib contributions. These terms are calculated using standard equations of statistical mechanics [113]. Srot is a function of the moments of inertia of the molecule, whereas Strans is a function of the mass and the solute concentration. Strans is the only term in the free energy of an ideal solution that depends on solute concentration, leading to the concentration
dependence of the binding process. The vibrational entropy Svib is calculated with the quantum formula from a normal mode analysis (NMA) [114]. A quasi
harmonic analysis of MD simulations is also possible.

Recent Advances in Protein-Ligand Interactions

However, it has been found that it does not always yield convergent values, even using long MD simulation trajectories, and leads to large deviations from the results obtained with NMA, giving an overall unreasonable entropic contribution [114].

Current Computer-Aided Drug Design, 2013, Vol. 9, No. 4

7

reduction of noise arising from flexible remote regions relative to the binding site, and conformational restraints imposed by the complex geometry. Thus, this one
simulation variant is attractive when Hintra may be reasonably neglected. Comparisons between one and three trajectories results can be found in literature [116, 117]. The recent applications of the MM/PBSA method can be found in some recent studies [118, 119]. For some systems, MM/PBSA fails to predict the binding free energy by an order of magnitude, and a new method MM/QMSA (or MM(QM)/COSMO as it is initially called) from the Cavasotto group, in which the whole protein is described by QM calculations, is used to calculate the binding free energy [120]. This method is already established to obtain the binding free energies of a series of phosphopeptides to the SH2 domain of human LCK [121]. Parameters for the COSMO intrinsic solvation model have also been derived [122]. 7.2. LIE Methods The linear interaction energy (LIE) method is first proposed by Aqvist and coworkers [123] to calculate the binding affinity of protein
ligand complexes. It is also known that the linear response method belongs to the semi
empirical methods. Two simulations are required for estimating the absolute free energies, one for ligand in solution and other for ligand in protein binding site. The snapshots saved from the simulations represent Boltzmann ensembles of conformations and are used to compute the Boltzmann
averaged electrostatic and van der Waals interaction energies of ligand with its environments in the bound and free states. The binding free energy is estimated as,

Fig. (3). The protocol for calculations of binding free energies (for the interpretation of references to colors in this figure, the reader is referred to the web version of this article).

Generally, in MM/PBSA methods we run MD simulation till 200
500 frames in explicit solvent and then the trajectories of the MD simulations are averaged. However its GB variant, MM/GBSA, is often performed on single snapshots especially in industrial uses. The ions and the explicit solvent are striped before the MM/PBSA analysis and periodic boundary conditions are applied. A study in the implicit solvent has been performed by Pearlman and Rizzo [115, 116], and it is found that MM/PBSA yields better results with MD simulations restrained around the X
ray structure, compared to unrestrained simulations. There are two possibilities regarding to the MD simulation for MM/PBSA. In the first approach one should make three trajectories, one for the complex and two for isolated molecules. However, a popular alternative is to perform only one MD simulation for the complex. In this variant, the terms relative to one isolated molecule are calculated after removing the atoms of the other molecule in the frame extracted from MD simulation of the complex. As a consequence, the reorganization energy of the molecules upon association is neglected (Hintra = 0). However, this variant is less computationally demanding and leads to an increase of convergence due to cancellation of errors,

G = (
)
(
) + 

+

As usual, the angle bracket indicates the Boltzmann average.
and  are two parameters, whereas  is a constant that can be added to get a corrected energy. To determine the binding free energy G one needs to perform two simulations, one of ligand in solvent and the other of ligand in bounded form within protein. The first term describes the electrostatic contributions according to the linear response approximation (LRA) theory [123, 124]. The second term holds for the nonpolar contributions. Its linear relationship with the surrounding van der Waals energy is based on the observation that the solvation energies of non
polar compounds are linearly correlated with the surrounding van der Waals energies [125, 126]. In the initial implementation,  is fixed to 0.5 following the LRA approximation, while
is empirically fitted at 0.16 to reproduce the experimental activity of four structurally related endothiapepsin inhibitors [127].  is kept to 0 to limit the over
parameterization. Although these parameters give satisfying results for protein
ligand systems, it has been found that using the FEP calculations,  can be considered to have a function of the ligand nature, suggested to be 0.50 for ionic molecules and 0.43, 0.37, and 0.33 for neutral compounds with one, two and more hydroxyl groups, respectively [128]. A value of 0.18 is found to be optimal for
. Non
zero values are necessary for  in order to reproduce G for some systems

8 Current Computer-Aided Drug Design, 2013, Vol. 9, No. 4

[129]. More recently, it is suggested that  can be expressed as a function of the buried solvent accessible surface area (SASA) of the ligand that is buried upon complexation, leading to the modified equation, G = (
) + (
) + (
) However, in the study of Aqvist et al. [130, 131] the addition of this term is questioned since buried SASA is correlated with van der Waals term. Some efforts have been done to reduce the computational costs where the electrostatic terms are replaced by the function of Coulombic interaction between ligand and protein and the solvation reaction field energy [132]. This variant is attractive due to the reduction of computational costs. Alam et al. [133] show the reliability of the LIE method. A training set of 76 podophyllotoxin analogues is used to build a binding affinity model for estimating the binding free energy for 36 inhibitors (test set) with diverse structural modifications. They find that the average root mean square error (RMSE) between the experimental and predicted binding free energy values is merely 0.56 kcal/mol, This is comparable to the level of accuracy achieved by the most accurate computational methods, such as the free energy perturbation (FEP) or thermodynamic integration (TI). The squared correlation coefficient between the experimental and SGB
LIE estimations for the binding free energies of the test set compounds are also significant. The similar applications are also found in many other studies [134, 135]. A few important contributions to molecular recognition have been neglected in LIE, such as the conformational rearrangement upon complexation of ligand and receptor, the receptor desolvation energy, and the entropy. However, it has been argued that these terms are implicitly taken into account by the LRA approximation and the adjustable parameters [136]. The quality of the results obtained by LIE methods is somewhat surprising, because the LIE method does not account explicitly for the standard configurational entropy or the internal energy of ligand. This method may be successful in part because it is generally used to compare ligands in a series of chemical structures. 7.3. FEP Methods Among all the known methods for the binding free energy calculations, the free energy perturbation (FEP) methods give the most promising results. However, due to the large computational time, these methods are not applicable for large systems. The change in the free energy between two states of a system (before and after binding) can be written as
RTIn, where U is the change in the energy function between the initial and final states, and the angle bracket indicates a Boltzmann average in the initial state. The change in the energy can be found by the interaction energy term and the Boltzmann average taken by a classical MD run. In reality, such a simulation is extraordinarily difficult to converge unless the initial and final states are simple enough. The FEP methods solve this problem by

Dubey et al.

breaking the change into small steps U (i.e., perturbations) and running separate simulations for each resulting energy function Ui to obtain the stepwise free energy changes. The perturbation methods can also be used to compute the standard binding free energy of protein and ligand. The free energy difference between two states A and B can formally be obtained from the Zwanzig’s formula [137], G = GB - GA = 
1InA where  = 1/KT and the angular bracket denotes a MD generated ensemble average. The main criterion for the above equation to be practically useful is that the configurations sampled on the potential VA should have a reasonable (at least non
vanishing) probability of occurrence also on VB. This means that thermally accessible regions of the two potentials should have a significant degree of overlap; otherwise, the results will be a very slow convergence. The convergence can be assessed by interchanging the labels A and B and change of the signs of G in the above equation, thus making the formula “backwards”. In order to solve the convergence problem, a multi
stage approach is normally adopted. A path between the states A and B is defined by introducing a set of intermediate potential energy functions that are usually constructed as the linear combinations of the initial (A) and final (B) state potentials, Vm = (1
m) VA + mVB where m varies from 0 to 1. In practice this path is discretized into a number of points (m = 1, ... , n), each being represented by a separate potential energy function that corresponds to a given value m. This coupling parameter approach rests upon the fact that the free energy difference is uniquely defined by the initial and final states (i.e. a state function) and can be computed along any reversible path connecting these states. Now the total free energy change can be obtained by summing over the intermediate states along the  variable, G = GB - GA = 
1In m This approach is generally referred to as the free energy perturbation (FEP) method. The FEP methods can also be used to compute the standard binding free energy for protein-ligand complexes. Many successful applications of the FEP methods can be found elsewhere [138]. Although the equation used in the FEP methods is exact, many studies have demonstrated that except in the cases of rather small changes, the equation converges as a function of the amount of data collected that should be far from ideal [138]. 7.4. TI Methods The energy expression of thermodynamic integration (TI) methods can be obtained by taking derivative of the free energy with respect to some continuous parameter  describing a series of intermediate alchemical states, dG⁄d = d⁄d
e -H() dx =  G =


Recent Advances in Protein-Ligand Interactions

where the pathways of intermediates between the states of interest is parameterized between  = 0 and  = 1. When the ending states have different masses, the momenta will be 
dependent as well. The TI trades variance for bias. Averaging over requires fewer uncorrected samples to reach a given level of relative error than averaging over e
H(x) , as long as is well behaved. However, to compute the total free energy from a series of individual simulations, we have to use some types of numerical integration for the integral, which by definition introduces bias for a number of different numerical techniques [139, 140]. For alchemical changes that result in the smooth, monotonic curves for , the TI methods can be quite accurate using a relatively small number of points. However, if the curvature becomes large as can frequently be in the alchemical simulations where Lennard-Jones potentials are turned on or off, the bias introduced by discretization of the integral can then become large [141, 142]. Even in the case of small curvature (i.e., charging of SPC water in water) the reasonably large errors can be caused; i.e., 5
10% of the total free energy with 5 values [143]. For relatively long simulations, the TI methods turn to be very bad approximation and can be avoided except in the case of Jarzynski's relationship. 8. CASE STUDIES In this section, some application cases of MD simulations and subsequent free energy calculations will be discussed. 8.1. MD and Free Energy Calculation of C
Abl Tyrosine Kinase Abl kinase plays an important role in the most fatal human pathogen chronic mylogenous leukemia (CML). There are several studies to clarify the related mechanisms [144-147] and suggest that a few regulatory elements should be responsible for the activity. The first is a centrally located segment termed as activation loop (Asp381-Pro402), whose phosphorylation state determines the activity of kinase. When phosphorylated, it permits access to the peptide substrate binding site and when not phosphorylated it can be disordered and assumed to present various stable conformations. A second regulatory element is the DFG (Asp381-Phe382-Glu383), which is a highly conserved motif located at the N
terminal of the activation loop. In the active kinase it is located inwards of the binding site, whereas in the inactive kinase it is flipped outwards to the binding site. A C
helix (Glu279-Glu292) is the third regulatory element which contains a glutamate residue (Glu286) that is verified to be crucial for the kinase activity. The conformational mechanism of Abl kinase complex with the known inhibitor, imatinib, can be explored by molecular dynamic simulations [148]. The study of Dubey et al. shows that Imatinib (STI) binds to the cleft of this inactivated kinase to form a stable complex, which results in the inhibition of biological activity of kinase [149]. The absolute free energy of 6.04 kcal/mol supports the tight binding of STI with c
Abl tyrosine kinase. A cavity is created by replacing the solvent molecule from the large SA surface area, which is confirmed by the cavity energy of 61.64 kcal/mol. The binding
free energy of 11.87

Current Computer-Aided Drug Design, 2013, Vol. 9, No. 4

9

kcal/mol obtained from the SIE method is in accordance with the experimental binding energy of 10.37 kcal/mol [150]. The conformational flexibility suggests that the activation loop is the most flexible, which facilitates the nucleotide binding and release processes. We have, for the first time, showed the computational flexibility of the regulatory elements of c
Abl protein kinase. The interactions between the STI and Abl core domains are qualified and quantified, whereas other physico
chemical parameters are often neglected, such as the energy of solvation, cavity energy, entropy contributions, which can also play a major role in the inhibition. It is clear that the binding energy calculated by the SIE methods gives very accurate results that are the closest to the experimental value. From the calculated binding
free energies of both methods, we can see that the van der Waals interaction energy represents the main contributions to the binding, which can be explained by a large binding site. It is concluded that the hydration plays a crucial role for the stability of the kinase
STI complex. The water molecules may become the life partner for the stability of the complex. Using MD simulations the authors propose that Glu21: Lys13, Glu52: Lys37, Glu82: Lys144 and Glu157: Lys170 salt bridges need to be conserved for the in-vivo stability of kinases and their optimal activity. The loop regions in kinase are the preferred sites targeted by ATP because of conformational adaptability required for the active
site binding. Therefore, the structural stabilization in the loop regions is critical for kinase functions. The mutational experiments for some of these salt bridges show the reduction in the activity because of loss of the conformational rigidity. These findings highlight the importance of considering kinase enzyme conformation in the rational design of inhibitors for cancer treatment. 8.2. Implementation of Hybrid QM/MM Energy in Binding Free Energy Calculations Although, the above computational methods are widely used and play a vital role in understanding structural functions, they have some limitations. These occur due to the neglect of electronic contributions that exert a substantial impact in the processes such as formation and/or cleavage of covalent bonds and fluctuation of charge during the MD simulations. The quantum mechanical (QM) approach gives an insight to overcome these limitations but it is not practical to be applied on the entire biomolecular systems. Thus, Quantum Mechanics/Molecular Mechanics (QM/MM) hybrid calculations have been used in the recent years in order to increase the computational accuracy [151-154]. In the QM/MM approach, we apply MM force fields to the receptor/protein molecule while the QM potentials for the ligand/core fragments (Fig. 4). The corresponding total energy function is written as, E = EMM + EQM + EQM/MM Here, EMM and EQM are the energies for the QM and the MM regions, which are calculated using the selected QM methods and force field equations, respectively. The last term, EQM/MM, describes the interaction between the QM and MM parts and typically contains the terms for electrostatic, van der Waals, and bonding interactions across the region.

10 Current Computer-Aided Drug Design, 2013, Vol. 9, No. 4

The use of classical MM/PBSA methods gives the binding affinity of protein
ligand complexes but ignores the energy contributions due to the EQM/MM term. In a recent study, the combined QM/MM methodology is incorporated

Dubey et al.

studies of binding events where quantum mechanical effects play a major role. The QM(MM)/PBSA approach presented here has a major advantage over purely classical mechanics
based (MM/PBSA) methods in that the force field parameters for ligands do not affect the results. The free energy calculation procedure can thus be fully automated, and this permits the atomic
based binding affinity calculations for a wide range of ligands with satisfying accuracy and reasonable CPU costs. 8.3. Stability of Qaudruplexes Studied Simulations and Free Energy Calculations

Fig. (4). A schematic representation of QM/MM calculations. The ligand (in circle) is treated as the QM (quantum mechanical) region while outside receptor is treated as the MM (molecular mechanical) region (for the interpretation of references to colors in this figure, the reader is referred to the web version of this article).

into the free energy calculations of protein
ligand binding using the AMBER software package [155]. It reveals that the binding free energy calculations using the classical MM/PBSA methods are frequently used to calculate the binding free energy of protein
ligand complexes but the results are less reliable due to the neglect of the electronic contributions. The contributions of electronic interactions and polarization effects due to the QM treatment are evident from the large values of EQM. Implementation of the QM/MM techniques for protein
ligand complexes is a fertile field of the computational researches, but the implementation of this technique for the exact calculation of binding affinity is really a troublesome task. To the best of our knowledge, this study represents one of the most successful applications of the AMBER 10 software package to calculate the QM/MM
PBSA binding free energy. The calculated binding free energies are in good agreement with the experimental results obtained for the c
Abl complex with Imatinib [150], and the results obtained from implementation of various QM potential methods show that the DFTB method gives the most accurate binding free energy. The disagreement with experiments probably may arise from the approximation involved in combining a QM/MM treatment for protein
ligand interaction with purely classical treatment of solvation, from approximation inherent to various methods as well as from the lack of a protein polarization term in these potentials. The accuracy of the QM/MM/PBSA approach will be improved by optimizing the implicit solvent model and using higher level QM methods. Our findings show a requirement for considering the electronic contributions to the binding affinity calculations which are frequently neglected during the classical MM/PBSA calculations. The methodology thus developed provides a useful approach to the atomic structure
based ligand
binding calculations and thereby opens up a new field of applications by permitting accurate computational

by

MD

Besides the protein
ligand interaction studies, MD simulations and binding free energy are also widely used for other biological systems like drug
DNA, protein
protein, quadruplexes etc. A study has been performed by Chaubey et al. on the modified quadruplexes to investigate the stability of different modelled telomeric quadruplexes. The structure and stability of telomeres are of great interest as they are closely related with cancer [156-159], aging [160, 161] and genetic stability [162-164]. The single stranded G
rich telomeric sequences at the end of chromosomes can fold into intramolecular quadruplex structures, a DNA secondary structure consisting of stacked G
tetrad planes connected by a network of Hoogsteen hydrogen bonds and stabilized by mono-valent cations, such as Na+ and K+. Human telomeric DNA repeats are highly conserved, which has been suggested to be related to their ability to form DNA G
quadruplexe [165]. The formation and stabilization of the DNA G
quadruplex in the human telomeric sequence have been shown to inhibit the activity of telomerase. Thus, the telomeric DNA quadruplex has been considered to be an attractive target for cancer therapeutic intervention [156, 158, 159, 166-168]. Chaubey et al. [169] have used the molecular modeling methods to build models of G
quadruplex arrays of monomeric structures formed from extended human telomeric DNA. Different reports suggest variable conformations of human telomeric quadruplex in the presence of K+ ions. Seven different types of human telomeric d(TAGGGT)4 modified quadruplex conformations, viz., parallel propeller type (all strands are parallel), basket
type, chair
type and mixed hybrid types, are used for the simulation studies to understand their relative stabilities and interactions between strands. Furthermore, in order to characterize the most energetically favorable quadruplex structures, they have estimated the binding free energy and entropy. Only a few studies have been reported based on the MD simulations of human telomeric quadruplex and the role of LNA with DNA in duplex [170, 171], triplex and quadruplex. The relative stability of the modified d(TAGGGT)4 G
quadruplex configurations is clarified. The MD simulations are performed for seven different modified quadruplex structures. Among all the modified structures, one of the mixed hybrid type structures possessing two parallel and two anti-parallel strands shows the lowest free energy. The accompanying energy calculations suggest that this mixed hybrid type conformation is the most structurally stable. The RMS deviation of MOD2 complex, which is also a mixed hybrid type structure, is smaller, whereas the RMS deviations of other complexes are much higher and vary significantly during the MD simulations. This suggests that

Recent Advances in Protein-Ligand Interactions

Current Computer-Aided Drug Design, 2013, Vol. 9, No. 4 [6]

the modified monomeric DNA sequence d(TAGGGT)4 can form only mixed hybrid type of structure, and the modified sequences may help to inhibit the growth of the tumor and cancerous cell. The tendency of human telomeric G
quadruplex structures to adopt different conformations under the various environmental conditions represents grand challenges in the design of specific drug candidates that can efficiently bind and stabilize these structures.

[10]

9. CONCLUDING REMARKS

[11]

With increase of computational efficiency and robotic techniques, MD simulations that quantify at the molecular level the physical phenomena modeled by statistical simulations have emerged as a frequently used method providing a tangible link between theory and experiment. With twenty years of “hindsight” gained from the methodological development and characterization, a variety of problems of chemical and biological relevance can now be tackled with confidence. Hopefully, we will be able to perform MD simulations for the folding of proteins of large size at micro
to milli
second levels. Using the MD simulations, the structure of the target (or target complex), if not available, can be derived easily from the folding simulations. Alternatively, simulations can be used to prepare homology models as accurate as those of the experimental structures. The free energy methods, in the near future, will turn to be an element that can be avoided in screening pipelines and discriminating between candidates. The implementation of quantum mechanical (QM) methods in the free energy calculations will drastically enhance the accuracy of the binding free energy and we can easily elaborate the polarization and electronic effect in protein
ligand binding as well as protein
protein interactions in the upcoming years.

[7] [8] [9]

[12] [13]

[14] [15]

[16] [17]

[18]

[19]

[20] [21]

CONFLICT OF INTEREST The authors confirm that this article content has no conflict of interest.

[22]

ACKNOWLEDGEMENTS

[24]

Authors acknowledge DST for the FIST scheme to the department. KDD acknowledges CSIR, New Delhi for SRF, while RKT and RPO acknowledge the support from UGC, New Delhi.

[25]

REFERENCES

[28]

[1]

[29]

[2]

[3]

[4]

[5]

Wlodawer, A. Rational approach to aids drug design through structural biology. Annu. Rev. Med., 2002, 53, 595-614. Wang, J.M.; Hou, T.J.; Xu, X.J. Recent advances in free energy calculations with a combination of molecular mechanics and continuum models. Curr. Comput. Aided Drug. Des., 2006, 2, 287-306. Wang, W.; Donini, O.; Reyes, C.M.; Kollman, P.A. Biomolecular simulations: Recent developments in force fields, simulations of enzyme catalysis, protein-ligand, protein-protein, and protein-nucleic acid noncovalent interactions. Annu. Rev. Biophys. Biomol. Struct., 2001, 30, 211-243. Kollman, P.A.; Massova, I.; Reyes, C.; Kuhn, B.; Huo, S.H.; Chong, L.; Lee, M.; Lee, T.; Duan, Y.; Wang, W.; Donini, O.; Cieplak, P.; Srinivasan, J.; Case, D.A.; Cheatham, T.E. Calculating structures and free energies of complex molecules: Combining molecular mechanics and continuum models. Acc. Chem. Res., 2000, 33, 889-897. Gilson, M.K.; Zhou, H.X. Calculations of protein binding affinities. Ann. Rev. Biophys., 2007, 36, 21-42.

[23]

[26] [27]

[30]

[31]

[32] [33]

11

Christ, C.D.; Mark, A.E.; van Gunsteren, W.F. Basic ingredients of free energy calculations: A review. J. Comput. Chem., 2010, 31, 1569-1582. Zhou, H.X.; Gilson, M.K. Theory of free energy and entropy in noncovalent binding. Chem. Rev., 2009, 109, 4092-4107. Shirts, M.R. Best practices in free energy calculations for drug design. Meth. Mol. Biol., 2012, 819, 425-467. Beveridge, D.L.; Dicapua, F.M. Free-energy via Molecular simulation applications to chemical and biomolecular systems. Annu. Rev. Biophys. Biomol. Struct., 1989, 18, 431-492. Aqvist, J.; Medina, C.; Samuelsson, J.E. New method for predicting binding-affinity in computer-aided drug design. Protein Eng., 1994, 7, 385-391. McCammon, J.A.; Gelin, B.R., Karplus, M. Dynamics of folded proteina. Nature, 1977, 267, 585-590. Berger, A; Linderstrom-Lang, K. Deutrium exchange of poly-dl-alanine in aqueous solution. Arch. Biochem. Biophys., 1957, 69, 106-118. Smith, J; Causack, S; Pezzeca, U; Brooks, B.R; Karplus, M. Inelastic neutron scattering analysis of low frequency motion in protein. A normal mode study of the bovine pancreative trypsine inhibitors. J. Chem. Phys., 1984, 85, 3636-3654. Brunger, A.T; Brooks, C.L III; Karplus, M. Active site of ribonuclease. Proc. Nat. Acad. Sci., USA, 1985, 82, 8458-8462. Frauenfelder, H.; Hartmann, H.; Karplus, M.; Kuntz, I.D.Jr.; Kuriyan, J.; Parak, F.; Petsko, G.A.; Ringe, D.; Tilton, R.F.Jr.; Connolly, M.L.; Max, N. Thermal expansion of protein. Biochemistry, 1987, 26, 254261. Brunger, A.T; Kuriyan, J; Karplus, M. Crystallographic R factor refinement by molecular dynamics. Science, 1987, 235, 458-460. Nilson, L; Clore, G.M; Gronenborn, A.M; Brunger, A.T; Karplus M. Structure refinement of oligonucleotides by molecular dynamics with nuclear Overhauser effect interproton distance restraints: application to 5' d(C-G-T-A-C-G)2. J. Mol. Biol., 1986, 188, 455-475, Colonna-Ceseri, F.; Perahia, D.; Karplus, M.; Eklund, H.; Brädén, C.I.; Tapia, O. Interdomain motion in liver alcohol dehydrogenease: structural and energetics of the hinge bending mode. J. Biol. Chem., 1986, 261, 15273-15280. Harvey, S.C.; Prabhakaran, M.; Mao, B.; MacCammon, J.A. Phenylalanine transfer RNA: molecular dynamics simulation. Science, 1984, 223, 1189-1191. Case, D.A; Karplus M. Dynamics to ligand binding to heme proteins. J. Mol. Biol., 1979, 132, 343-368. Brooks, B.R; Karplus, M. Harmonics dynamics of proteins: normal modes and fluctuations in bovine pancreatic trypsine inhibitor. Proc. Nat. Acad. Sci. USA., 1983, 80, 6571-6575. Karplus, M.; McCammon, J.A. Molecular dynamic simulation of biomolecules. Nat. Struct. Biol., 2002, 9, 646-652. Berendsen, H.J.C; Hayword S. Collective protein dynamics in relation to function. Curr. Opin. Struct. Biol., 2000, 10, 165-169. Hayword, S; Kitao, A; Go, N. Harmonicity and anharmonicity in protein dynamics. Prot. Struct. Func. Gent., 1995, 23, 177-186. Nose, S. A unified formulation of the constant temperature moleculardynamics methods. J. Chem. Phys., 1984, 81, 511-519. Andersen, H.C. Molecular dynamics at constant pressure and/or temperature. J. Chem. Phys., 1980, 72, 2384-2393. Izaguirre, J.A; Catarello, D.P; Wozniak, J.M; Skeel, R.D. Langevin stabilization of molecular dynamics. J. Chem. Phys., 2001, 114, 20902098. Hünenberger, P.H; Thermostats algorithms for molecular dynamics simulations. Adv. Polym. Sci., 2005, 173, 105-149. Harvey, S.C; Tan, R.K.Z; Cheatham, T.E. The flying ice cube: Velocity rescaling in molecular dynamics leads to violation of energy equipartition. J. Comp. Chem., 1998, 19, 726-740. Cerutti, D.S; Duke, R; Freddolino, P.L; Fan, H; Lybrand, T.P. Vulnerability in popular molecular dynamics packages concerning langevin and andersen dynamics. J. Chem. Theory. Comput., 2008, 14, 1669-1680. Piana, S; Lindorff-Larsen, K; Shaw, D.E. How robust are protein folding simulations with respect to force field parameterization? Biophys. J., 2011, 100, L47-L49. Duan, Y; Harvey, S.C; Kollman, P.A. Protein folding and beyond. In: Keinan E, Schechter I (eds) Chemistry for the 21st century. Wiley, Weinheim., 2001. Babine, R.E; Bender, S.L. Molecular recognition of protein-ligand complexes:
applications to drug design. Chem. Rev., 1997, 97, 13591472.

12 Current Computer-Aided Drug Design, 2013, Vol. 9, No. 4 [34]

[35] [36] [37] [38] [39]

[40] [41]

[42] [43]

[44] [45]

[46]

[47]

[48]

[49] [50]

[51]

[52]

[53] [54]

[55] [56]

[57]

Gohlke, H; Klebe, G. Approaches to the description and prediction of the binding affinity of small-molecule ligands to macromolecular receptors. Angew. Chem. Int. Ed. Engl., 2002, 41, 2644-2676 McDonald, I. K; Thornton, J.M. Satisfying Hydrogen Bonding Potential in Proteins. J. Mol. Biol., 1994, 238, 777-793. Dunitz, J.D. Win some, lose some: enthalpy-entropy compensation in weak intermolecular interactions. Chem. Biol., 1995, 2, 709-12. Gilli, P.; Ferretti, V.; Gilli, G; Brea, P.A. Enthalpy-entropy compensation in drug-receptor binding. J. Phys. Chem., 1994, 98, 15151518. Williams, D.H.; O’Brien, D.P.; Bardsley, B. J. Am. Chem. Soc., 2001, 123, 737-738. Weber, P.C.; Wendoloski, J.J.; Pantoliano, M.W.: Salemme, F.R. Crystallographic and thermodynamic comparison of natural and synthetic ligands bound to streptavidin. J. Am. Chem. Soc., 1992, 114, 3197-3200. Born, M. Volumes and heats of hydration of ions. Z. Phys. 1920, 1, 4548. Kretschmer, R.; Kinzel, D.; and González, L. The role of hydrogen bonds in protein-ligand interactions. DFT calculations in 1,3dihydrobenzimidazole-2 thione derivatives with glycinamide as model HIV RT inhibitors. Int. J. Quantum Chem., 2012, 112, 1786-1795. Panigrahi, S.K.; Desiraju, G.R.; Strong and weak hydrogen bonds in the protein-ligand interface. Proteins, 2007, 67, 128-141. Levy, Y.; Onuchic, J.N. Water mediation in protein folding and molecular recognition. Annu. Rev. Biophys. Biomol. Struct., 2006, 35, 389-415. Kauzmann, W. Some factors in the interpretation of protein denaturation. Adv. Protein Chem., 1959, 14, 1-63. Shah, F.; Gut, J.; Legac. J.; Shivakumar, D.; Sherman, W.; Rosenthal, P.J.; Avery, M.A. Computer-aided drug design of falcipain inhibitors: virtual screening, structure-activity relationships, hydration site thermodynamics, and reactivity analysis. J. Chem. Inf. Model., 2012, 52, 696-710. Beuming, T.; Che, Y.; Abel, R.; Kim, B.; Shanmugasundaram, V.; Sherman, W. Thermodynamic analysis of water molecules at the surface of proteins and applications to binding site prediction and characterization. Proteins, 2012, 80, 871-883. Snyder, P.W.; Mecinovi, J.; Moustakas, D.T.; Thomas, III. S.W.; Harder, M.; Mack, E.T.; Lockett, M.R.; Héroux A, Sherman W, Whitesides GM. Mechanism of the hydrophobic effect in the biomolecular recognition of arylsulfonamides by carbonic anhydrase. PNAS., 2011, 108, 17889-17894. Abel, R.; Salam, N.K.; Shelley, J.; Farid, R.; Friesner, R.A.; and Sherman, W. Contribution of explicit solvent effects to the binding affinity of small-molecule inhibitors in blood coagulation factor serine proteases. Chem. Med. Chem., 2011, 6, 1049-1066. Robinson, D.D.; Sherman, W.; and Farid, R. Understanding kinase selectivity through energetic analysis of binding site waters. Chem. Med. Chem., 2010, 5, 618-627. Higgs, C.; Beuming, T.; Sherman W. Hydration site thermodynamics explain SARs for triazolylpurines analogues binding to the A2A receptor. ACS Med. Chem. Lett., 2010, 1, 160-164. Pearlstein, R.A.; Hu, Q.Y.; Zhou, J.; Yowe, D.; Levell, J.; Dale, B.; Kaushik, V.K.; Daniels, D.; Hanrahan, S.; Sherman, W.; Abel, R. New hypotheses about the structure-function of proprotein convertase subtilisin/kexin type 9: Analysis of the epidermal growth factor-like repeat A docking site using WaterMap. Proteins, 2010, 78, 2571-2586. Beuming, T.; Farid, R.; Sherman, W. High-energy water sites determine peptide binding affinity and specificity of PDZ domains. Protein Sci., 2009, 18, 1609-1619. Abel, R.; Young, T.; Farid, R.; Berne, B.J.; Friesner, R.A; The role of the active-site solvent in the thermodynamics of factor Xa ligand binding. J. Am. Chem. Soc., 2008, 130, 2817-2831. Young, T.; Abel, R.; Kim, B.; Berne, B.J.; Friesner, R.A. Motifs for molecular recognition exploiting hydrophobic enclosure in proteinligand binding. PNAS., 2007, 104: 808-813. Alberty, R.A. Effect of pH on Protein-Ligand Equilibria. J. Phys. Chem. B, 2000, 104, 9929-9934. Katre, U.V.; Mazumder, S.; Prusti, R.K.; Mohanty, S. Ligand binding turns moth pheromone-binding protein into a pH sensor: effect on the Antheraea polyphemus PBP1 conformation. J. Biol. Chem., 2009, 284, 32167-32177. Mongan, J.; Case, D.A.; McCammon, J.A. Constant pH molecular dynamics in generalized Born implicit solvent. J. Comput. Chem., 2004, 25, 2038-2048.

Dubey et al. [58] [59]

[60] [61]

[62] [63]

[64] [65] [66]

[67] [68] [69] [70]

[71] [72]

[73]

[74] [75]

[76] [77]

[78]

[79]

[80] [81]

[82] [83]

Machuqueiro M, Baptista AM. Acidic range titration of HEWL using a constant-pH molecular dynamics method. Proteins, 2008, 72, 289-298. Baptista, A.M.; Teixeira, V.H.; Soares, C.M. Constant-pH molecular dynamics using stochastic titration. J. Chem. Phys., 2002, 117, 41844200. Dlugosz, M.; Antosiewicz, J.M. Constant-pH molecular dynamics simulations: a test case of succinic acid. Chem. Phys., 2004, 302, 161170. Machuqueiro, M.; Baptista, A.M. Constant-pH molecular dynamics with ionic strength effects. Protonation-conformation coupling in decalysine. J. Phys. Chem. B., 2006, 110, 2927-2933. Mongan, J.; Case, D.A.; McCammon, J.A. Constant pH molecular dynamics in generalized Born implicit solvent. J. Comput. Chem., 2004, 25, 2038-2048. Dugosz, M.; Antosiewicz, J.M. Constant-pH molecular dynamics study of protonation - structure relationship in a heptapeptide derived from ovomucoid third domain. Phys. Rev. E. Stat. Nonlin. Soft Matter Phys., 2004, 69, 021915. Chaires, J.B. Calorimetry and thermodynamics in drug design, Annu. Rev. Biophys., 2008, 37, 135-151. Freire, E. Do enthalpy and entropy distinguish first in class from best in class? Drug Discov. Today, 2008, 13, 869-874. Lipinski, C.A. Drug-like properties and the causes of poor solubility and poor permeability. J Pharmacol Toxicol Methods., 2000, 44, 235249. Go, N.; Go, M.; Scheraga, H.A. Molecular theory of the helix-coil transition in polyamino acids, I. Formulation. Proc. Natl. Acad. Sci. USA., 1968, 59, 1030-1037. Karplus, M., Kushick, J.N. Method for estimating the configurational entropy of macromolecules. Macromolecules, 1981, 14, 325-332. Karplus, M.; Ichiye, T.; Pettitt, B.M. Configurational entropy of native proteins. Biophys. J., 1987, 52, 1083-1085. Frederick, K.K.; Marlow, M.S.; Valentine, K.G.; Wand, A.J. Conformational entropy in molecular recognition by proteins. Nature, 2007, 448, 325-329. Tzeng, S.R.; Kalodimos, C.G. Dynamic activation of an allosteric regulatory protein. Nature, 2009, 462, 368-372. Diehl, C.; Engstrom, O.; Delaine, T.; Hakansson, M.; Genheden, S.; Modig, K.; Leffler, H.; Ryde, U.; Nilsson, U.J.; Akke, M. Protein flexibility and conformational entropy in ligand design targeting the carbohydrate recognition domain of galectin-3. J. Am. Chem. Soc., 2010, 132, 14577-14589. Marlow, M.S.; Dogan, J.; Frederick, K.K.; Valentine, K.G.; Wand, A.J. The role of conformational entropy in molecular recognition by calmodulin. Nat. Chem. Biol., 2010, 6, 352-358. Chang, C.E.; Chen, W.; Gilson, M.K. Ligand configurational entropy and protein binding. Proc. Nat. Acad. Sci. USA., 2007, 104, 1534-1539. DeLorbe, J.E.; Clements, J.H.; Teresk, M.G.; Benfield, A.P.; Plake, H.R.; Millspaugh, L.E.; Martin, S.F. Thermodynamic and structural effects of conformational constraints in protein-ligand interactions. Entropic paradoxy associated with ligand preorganization. J. Am. Chem. Soc., 2009, 131, 16758-16770. Mann, A. In: Wermuth CG (ed) The Practice of Medicinal Chemistry, 2nd edn, Academic Press, London., 2003. Showalter, S.A.; Johnson, E.; Rance, M.; Bruschweiler, R. Toward quantitative interpretation of methyl side-chain dynamics from NMR by molecular dynamics simulations, J. Am. Chem. Soc., 2007, 129, 1414614147. Krishnan, M.; Smith, J.C. Response of small-scale, methyl rotors to protein- ligand association: a simulation analysis of calmodulin-peptide binding. J. Am. Chem. Soc., 2009, 131, 10083-10091. Diehl, C.; Genheden, S.; Modig, K.; Ryde, U.; Akke, M. Conformational entropy changes upon lactose binding to the carbohydrate recognition domain of galectin-3. J. Biomol. Nmr., 2009, 45, 157-169. Li, D.W.; Bruschweiler, R.A dictionary for protein side-chain entropies from NMR order parameters. J. Am. Chem. Soc., 2009, 131, 72267227. Trbovic, N.; Cho, J.H.; Abel, R.; Friesner, R.A.; Rance, M.; Palmer, A.G.3rd. Protein side-chain dynamics and residual conformational entropy. J. Am. Chem. Soc., 2009, 131, 615-622. Teague, S.J. Implications of protein flexibility for drug discovery. Nat. Rev. Drug. Discov., 2003, 2, 527-541. Cavasotto, C.N.; Orry, A.J. Ligand docking and structure-based virtual screening in drug discovery. Curr. Top. Med. Chem., 2007, 7, 10061014.

Recent Advances in Protein-Ligand Interactions [84]

[85] [86]

[87] [88]

[89] [90] [91] [92]

[93]

[94] [95] [96] [97]

[98] [99] [100]

[101]

[102] [103]

[104]

[105]

[106] [107] [108]

[109]

B-Rao, C.; Subramanian, J.; Sharma, S.D. Managing protein flexibility in docking and its applications. Drug. Discov. Today., 2009, 14, 394400. Chang, M.W.; Belew, R.K.; Carroll, K.S.; Olson, A.J.; Goodsell, D.S. Empirical entropic contributions in computational docking: evaluation in APS reductase complexes. J. Comput. Chem., 2008, 29, 1753-1761. Huang, S.Y.; Zou, X. Inclusion of solvation and entropy in the knowledge-based scoring function for protein-ligand interactions. J. Chem. Inf. Model., 2010, 50, 262-273. Cosconati, S.; Forli, S.; Perryman, A.L.; Harris, R.; Goodsell, D.S.; Olson, A.J. Virtual screening with AutoDock: theory and practice. Expert. Opin. Drug Discov., 2010, 5, 597-607. Baron, R.; McCammon, J.A. Thermodynamic role of receptor flexibility, entropy, and motional correlation in protein
ligand binding. Chemphyschem., 2008, 9, 983-988. Crespo, A.; Fernandez, A. Induced disorder in protein-ligand complexes as a drug-design strategy. Mol. Pharm., 2008, 5, 430-437. Edholm, O.; Berendsen, H.J.C. Entropy estimation from simulations of non-diffusive systems. Mol. Phys., 1984, 51, 1011-1028. Schlitter, J. Estimation of absolute and relative entropies of macromolecules using the covariance matrix. Chem. Phys. Lett., 1993, 215, 617-621. Andricioaei, I.; Karplus, M. On the calculation of entropy from covariance matrices of the atomic fluctuations. J. Chem. Phys., 2001, 115, 6289. Baron, R.; van Gunsteren, W.F.; Hunenberger, P.H. Estimating the configurational entropy from molecular dynamics simulations: Anharmonicity and correlation corrections to the quasi-harmonic approximation. Trends. Phys. Chem., 2006, 11, 87-122. Teage S.J. Implications of protein flexibility for drug discovery. Nat. Rev. Drug Discov., 2003, 2, 527-541. Tepeschy, P.D.; Astay, M.; Ceder. G. Computation of configurational entropy using Monte Carlo probabilities in cluster-variation method entropy expressions. Model. Simul. Mater. Sci. Eng., 1998, 6, 787-797. Karplus, M.; Kushick, J.N. Method for estimating the configurational entropy of macromolecules. Macromolecules, 1981, 14, 325-332. Levy, R.M.; Karplus, M.; Kushick, J.; Perahia, D. Evaluation of the configurational entropy for proteins: application to molecular dynamics simulations of an a-Helix. Macromolecules, 1984, 17, 1370-1374. Chang, C.E.; Chen, W.; Gilson, M.K. Evaluating the accuracy of the quasiharmonic approximation. J. Chem. Theory. Comput., 2005; 1: 1017. Hnizdo, V.; Darian, E.; Fedorowicz, A.; Demchuk, E.; Li, S.; Singh, H. J. Comput. Chem. 2007; 28: 655. Killian, B.J.; Kravitz, J.Y.; Gilson, M.K. Extraction of configurational entropy from molecular simulations via an expansion approximation. J. Chem. Phys., 2007, 127, 024107. Hnizdo, V, Jun, T.; Killian, B.J.; Gilson, M.K. Efficient calculation of configurational entropy from molecular simulations by combining the mutual-information expansion and nearest-neighbor methods. J. Comput. Chem., 2008, 29, 1605-1614. Chong, S.H.; Ham, S. Configurational entropy of protein: A combined approach based on molecular simulation and integral-equation theory of liquids. Chem. Phys. Lett., 2011, 504, 225-229. Hensen, U.; Lange, O.F.; Grubmuller, H. Estimating absolute configurational entropies of macromolecules: the minimally coupled subspace approach. PLoS ONE, 2010, 5, e9179. Cozzini, P.; Kellogg, G.E.; Spyrakis, F.; Abraham, D.J.; Costantino. G.; Emerson, A.; Fanelli, F.; Gohlke, H.; Kuhn, L.A.; Morris, G.M.; Orozco, M.; Pertinhez, T.A.; Rizzi, M.; Sotriffer, C.A. Target flexibility: an emerging consideration in drug discovery and design. J. Med. Chem., 2008, 51, 6237-6255. Mangoni, R.; Roccatano, D.; Di Nola, A. Docking of flexible ligands to flexible receptors in solution by molecular dynamics simulation. Proteins, 1999, 35, 153-162. Lin, J.H.; Perryman, A.L.; Schames, J.R.; McCammon, J.A. Computational drug design accommodating receptor flexibility: the relaxed complex scheme. J. Am. Chem. Soc., 2002, 124, 5632-5633. Grunberg, R.; Nilges, M.; Leckner J. Flexibility and conformational entropy in protein-protein binding. Structure, 2006, 14, 683-693. Gohlke, H.; Kuhn, L.A.; Case, D.A. Change in protein flexibility upon complex formation: Analysis of Ras-Raf using molecular dynamics and a molecular framework approach, Proteins, 2004, 56, 322-337. Perdih, A.; Bren, U.; Solemajer, T. Binding free energy calculations of N-sulphonyl-glutamic acid inhibitors of MurD ligase. J. Mol. Mod., 2009, 15, 983-996.

Current Computer-Aided Drug Design, 2013, Vol. 9, No. 4 [110]

[111] [112] [113]

[114]

[115] [116]

[117]

[118]

[119] [120]

[121]

[122]

[123] [124]

[125] [126]

[127] [128]

[129]

[130] [131] [132]

[133]

13

Gohlke, H.; Case, D.A. Converging free energy estimates: MM-PB (GB)/SA studies on the protein -protein complex Ras-Raf. J. Comput. Chem., 2004, 25, 238-250. Fogolari, F.; Brigo, A.; Molinari H. Protocols for MM/PBSA molecular dynamics simulations of proteins. Biophys. J., 2003, 85,159-166. Grochowski, P.; Trylska, J. Continuum molecular electrostatics, salt effects, and counterion binding - a review of the Poisson-Boltzmann theory and its modifications. Biopolymers, 2008, 89, 93-113. Schwarzl, S.M.; Tschopp, T.B.; Smith, J.C.; Fischer, S. Can the calculation of ligand binding free energies be improved with continuum solvent electrostatic and ideal-gas entropy corrections? J. Comp. Chem., 2002, 23, 1143-1149. Tidor, B.; Karplus, M. The contribution of vibrational entropy to molecular association: the dimerization of insulin. J. Mol. Biol., 1994, 238, 405-414. Pearlman, D.A. Evaluating the molecular mechanics poissonboltzmann surface area free energy method using a congeneric series of ligands to p38 MAP kinase. J. Med. Chem., 2005, 48, 7796-7807. Rizzo, R.C.; Toba, S.; Kuntz, I.D. A molecular basis for the selectivity of thiadiazole urea inhibitors with Stromelysin-1 and Gelatinase-A from generalized born molecular dynamics simulations. J. Med. Chem., 2004, 47, 3065-3074. Shao, J.; Tanner, S.W.; Thompson, N.; Cheatham, T.E.III. Clustering the molecular dynamics trajectories:1. Characterizing the performance of different clustering algorithms. J. Chem. Theory. Comput., 2007, 3, 2312-2334. Balasubrimaniam, C.; Ojha, R.P.; Maiti, S. Hydration pattern of A4T4 and T4A4 416 DNA: a molecular dynamics study. Biochem. Biophys. Res. Com., 2007, 355, 1081-1086. Agrawal, S.; Ojha. R.P.; Maiti, S. Energetics of the human tel-22 quadruplex telomestatin interaction: a molecular dynamics study. J. Phys. Chem. B., 2008, 112, 6828-6836. Anisimov, V.M.; and Cavasotto, C.N. Quantum mechanical binding free energy calculation for phosphopeptide inhibitors of the Lck SH2 domain. J. Comput. Chem., 2011, 32, 2254-2263. Victor, A.; Ziemys, A.; Kizhake, S.; Yuan, Z.; Natarajan, A.; Cavasotto, C. Computational and experimental studies of the interaction between phospho-peptides and the C-terminal domain of BRCA1. J. Comp. Aided. Drug. Des., 2011, 25, 1071-1084. Anisimov, V.M.; Cavasotto, C.M. Hydration free energies using semiempirical quantum mechanical hamiltonians and a continuum solvent model with multiple atomic-type parameters. J. Phys. Chem. B., 2011, 115, 7896-7905. Åqvist, J.; Medina, C.; Samuelsson, J.-E. A new method for predicting binding affinity in computer-aided drug design. Protein Eng., 1994, 7, 385-391. Lee F.S.; Chu, Z.T.; Bolger, M.B.; Warshel. A. Calculations of antibody-antigen interactions: microscopic and semi-microscopic evaluation of the free energies of binding of phosphorylcholine analogs to McPC603. Protein Eng., 1992, 5, 215-228. Åqvist, J.; Hansson, T. On the validity of electrostatic linear response in polar solvents. J. Phys. Chem., 1996, 100, 9512-9521. Hansson, T.; Marelius, J.; Åqvist, J. Ligand binding affinity prediction by linear interaction energy methods. J. Comput. Aided Mol. Des., 1998, 12, 27-35. Carlson, H.A.; Jorgensen, W.L. An extended linear response method for determining free energies of hydration. J. Phys. Chem., 1995, 99, 10667-10673. Ljungberg, K.B.; Marelius, J.; Musil, D.; Svensson, P.; Norden, B.; Åqvist, J. Computational modelling of inhibitor binding to human thrombin Eur. J. Pharm. Sci., 12, 441-446. Zhou, R.H.; Friesner, R.A.; Ghosh, A.; Rizzo, R.C.; Jorgensen, W.L.; Levy, R.M. New linear interaction method for binding affinity calculations using a continuum solvent model. J. Phys. Chem. B., 2001, 105, 10388-10397. Åqvist, J.; Luzhkov, V.B.; Brandsdal, B.O. Ligand binding affinities from MD simulations. Acc. Chem. Res., 2002, 35, 358-365. Åqvist, J.; Marelius, J. The linear interaction energy method for predicting ligand binding free energies. J. Combin. Chem. High Through Scr., 2001, 4, 613-626. Carlsson, J.; Ander, M.; Nervall, M.; Åqvist, J. Continuum solvation models in the linear interaction energy method. J. Phys. Chem. B., 2006, 110, 12034-12041. Alam, M.A.; Naik, P.K. Applying linear interaction energy method for binding affinity calculations of podophyllotoxin analogues with tubulin

14 Current Computer-Aided Drug Design, 2013, Vol. 9, No. 4

[134]

[135]

[136] [137] [138]

[139]

[140]

[141] [142] [143]

[144]

[145] [146] [147]

[148]

[149]

[150]

Dubey et al.

using continuum solvent model and prediction of cytotoxic activity. J. Mol. Graph. Mod., 2009, 27, 930-947. Bortolato, A.; Moro, S. In silico binding free energy predictability by using the linear interaction energy (LIE) method:
Bromobenzimidazole CK2 inhibitors as a case study. J. Chem. Inf. Model., 2007, 47, 572-582. Singh, P.; Mhaka, A.M.; Christensen, S.B.; Gray, J.J.; Denmeade, S.R.; Isaacs, J.T. Applying linear interaction energy method for rational design of noncompetitive allosteric inhibitors of the sarco- and endoplasmic reticulum calcium-ATPase. J. Med.Chem., 2005, 48, 3005-3014. Kyani A.; Goliaei. B. Binding free energy calculation of peptides to carbon nanotubes using molecular dynamics with a linear interaction energy approach. J. Mol. Struct., 2009, 913, 63-69. Zwanzig, R. W. High‐temperature equation of state by a perturbation method. I. nonpolar gases J. Chem. Phys., 1954, 22, 1420-1426. Chandani, S.; Lee, C.H.; Loecher, E.L. Free-energy perturbation methods to study structure and energetics of dna adducts:
results for the major N2-dG adduct of benzo[a]pyrene in two conformations and different sequence contexts. Chem. Res. Tox., 2005, 18, 1108-1123. Shirts, M.R.; Pande, V.S. Comparison of efficiency and bias of free energies computed by exponential averaging, the bennett acceptance ratio, and thermodynamic integration. J. Chem. Phys., 2005, 122, 144107. Lu, N.D.; Singh, J.K.; Kofke, D.A. Appropriate methods to combine forward and reverse free-energy perturbation averages. J. Chem. Phys., 2003, 118, 2977-2984. Resat, H.; Mezei, M. Studies on free energy calculations. I. Thermodynamic integration using a polynomial path. J. Chem. Phys., 1993, 99, 6052-6061. Pitera, J.W.; van Gunsteren, W.F. A comparison of non-bonded scaling approaches for free energy calculations. Mol. Simulat., 2002, 28, 45-65. Lopex, M.A.; Kollman, P.A. Application of molecular dynamics and free energy perturbation methods to metalloporphyrin-ligand systems ii: Co and dioxygen binding to myoglobin. Prot. Sci., 1993, 2, 1975-1986. Nagar, B.; Hantschel, O.; Young, M.A.; Scheffzek, K.; Veach, D.; Bornmann, W.; Clarkson, B.; Superti-Furga, G.; Kuriyan, J. Structural basis for the autoinhibition of c-Abl tyrosine kinase. Cell, 2003, 112, 859-871. Nagar, B. c-Abl tyrosine kinase and inhibition by the cancer drug Imatinib. J. Nutr., 2007, 137, 1518S-1523S. Li, W.; Miller, W.T.; Role of activation loop tyrosine in regulation of the insulin-like growth factor I receptor tyrosine kinase. J. Biol. Chem., 2006, 281, 23785-23791. Emrick, M.A.; Lee, T.; Starkey, P.J.; Mumby, M.C.; Resing, K.A.; Ahn, N.G. The gatekeeper residue controles, auto activation of ERK2 via a pathway of intramolecular connectivity. Proc. Nat. Acad. Sci., 2006, 103, 18101-18106. Shan, Y., Seelinger, M.A.; Eastwood, M.P.; Frank, F.; Xu, H.; Jensen, M.O.; Dror, R.O.; Kuriyan, J.; Shaw, D.E. A conserved protonation dependent switch controls drug binding in the Abl kinase. Proc. Natl. Acad. Sci., 2009, 106, 139-144. Dubey, K.D.; Ojha, R.P. Conformational flexibility, binding energy, role of salt bridge and alanine-mutagenesis for c-Abl tyrosine kinase complex. J. Mol. Mod., 2012, 18, 1679-1689. Pricl, S.; Fermeglia, M.; Ferrone, M.; Tamborini, E. T315I mutated Bcr-Abl in chronic myeloid leukemia and Imatinib: Insights from a computational study. Mol. Cancer Ther., 2005, 4, 1167-1174.

Received: June 14, 2012

[151]

[152] [153] [154]

[155]

[156] [157] [158] [159] [160]

[161] [162] [163]

[164] [165]

[166] [167]

[168] [169] [170]

[171]

Lin, H.; Truhlar, D.G. QM/MM: What have we learned, where are we, and where do we go from here? Theor. Chem. Acc., 2007, 117, 185199. Senn, H.M.; Thiel, W. QM/MM studies of enzymes. Curr. Opin. Chem. Biol., 2007, 11, 182-187. Bruice, T.C. Computational approaches: reaction trajectories, structures, and atomic motions. Enzyme reactions and proficiency. Chem. Rev., 2006, 106, 3119-3139. Murphy, R.B.; Philipp, D.M.; Friesner, R.A. A mixed quantum mechanics/molecular mechanics (QM/MM) method for large-scale modeling of chemistry in protein environments. J. Comput. Chem., 2000, 21, 1442-1457. Dubey, K.D.; Ojha, R.P. Binding free energy calculation with QM/MM hybrid methods for Abl kinase inhibitors. J. Biol. Phys., 2011, 37, 6978. Mergny J.L, Helene C. G-quadruplex DNA: a target for drug design. Nat. Med., 1998, 4, 1366-1367. Sun, D.Y, Hurley, L.H. Methods in enzymology, drug nucleic acid interactions. Academic Press, Inc, 2001, 340: 573-592. Hurley, L.H. DNA and its associated processes as targets for cancer therapy. Nat. Rev. Can., 2002, 2, 188-200. Neidle, S.; Parkinson, G. Telomere maintenance as a target for anticancer drug discovery. Nat. Rev. Drug. Discov., 2002, 1, 383-393. Bodnar, A.G.; Ouellette, M.; Frolkis, M.; Holt, S.E.; Chiu, C.P.; Morin, G.B.; Harley, C.B.; Shay, J.W.; Lichtsteiner, S.; Wright WE. Extension of life-span by introduction of telomerase into normal human cells. Science, 1998, 279, 349-352. Harley, C.B. Telomere loss: mitotic clock or genetic time bomb? Mutat. Res., 1991, 256, 271-282. Sun, H.; Karow.; J.K.; Hickson, I.D.; Maizels, N. The Bloom’s syndrome helicase unwinds G4 DNA. J. Biol. Chem., 1998, 273, 27587-27592. Sun, H.; Bennett, R.J.; Maizels, N. The saccharomyces cerevisiae Sgs1 helicase esciently unwinds G-G paired DNAs. Nucleic Acids Res., 1999, 27, 1978-1984. Hackett, J.A.; Feldser, D.M.; Greider, C.W.; Telomere dysfunction increases mutation rate and genomic instability. Cell, 2001, 106, 275286. Salazar, M.; Thompson, B.D.; Kerwin, S.M.; Hurley, L.H. Thermally induced DNA:RNA hybrid to G-quadruplex transitions:possible implications for telomere synthesis by telomerase. Biochemistry, 1996, 35, 16110-16115. Hurley, L.H. Secondary DNA structures as molecular targets for cancer therapeutics. Biochem. Soc. Trans., 2001, 29, 692-696. Hurley, L.H.; Wheelhouse, R.T.; Sun, D.; Kerwin, S.M.; Salazar, M.; Fedoro, O.Y.; Han, F.X.; Han, H.Y.; Izbicka, E. et al. G-quadruplexes as targets for drug design. Pharmacol. Ther., 2000, 85, 141-158. Neidle, S.; Read, M.A. G-quadruplexes as therapeutic targets. Biopolymers, 2000, 56, 195-208. Chaubey, A.K.; Dubey, K.D.; Ojha, R.P. Stability and free energy calculations of LNA modified Human Tel-24 quadruplex structure: A molecular dynamics study. J. Comp. Aid. Mol. Des., 2012, 29, 289-299. Ivanova, A.; Rosch, N. The structure of LNA: DNA hybrids from molecular dynamics simulations: the effect of locked nucleotides. J. Phys. Chem. A, 2007, 111, 9307-9319. Pande, V.; Nilsson, L. Insights into structure, dynamics and hydration of locked nucleic acid (LNA) strand-based duplexes from molecular dynamics simulation. Nucleic Acids Res., 2008, 36, 1508-1516.

Revised: October 29, 2013

Accepted: April 28, 2013