Computational Molecular Bioscience
 Vol.2 No.1(2012), Article ID:18269,17 pages DOI:10.4236/cmb.2012.21002

Role of Ligand Reorganization and Conformational Restraints on the Binding Free Energies of DAPY Non-Nucleoside Inhibitors to HIV Reverse Transcriptase

Emilio Gallicchio

BioMaPS Institute for Quantitative Biology, Department of Chemistry and Chemical Biology, Rutgers, The State University of New Jersey, Piscataway, USA

Email: emilio@biomaps.rutgers.edu

Received January 28, 2012; revised February 21, 2012; accepted March 2, 2012

Keywords: HIV-RT; TMC125; TMC278; Etravirine; Rilpivirine; Binding Free Energy

ABSTRACT

The results of computer simulations of the binding of etravirine (TMC125) and rilpivirine (TMC278) to HIV reverse transcriptase are reported. It is confirmed that consistent binding free energy estimates are obtained with or without the application of torsional restraints when the free energies of imposing the restraints are taken into account. The restraints have a smaller influence on the thermodynamics and apparent kinetics of binding of TMC125 compared to the more flexible TMC278 inhibitor. The concept of the reorganization free energy of binding is useful to understand and categorize these effects. Contrary to expectations, the use of conformational restraints did not consistently enhance convergence of binding free energy estimates due to suppression of binding/unbinding pathways and due to the influence of rotational degrees of freedom not directly controlled by the restraints. Physical insights concerning the thermodynamic driving forces for binding and the role of “jiggling” and “wiggling” motion of the ligands are discussed. Based on these insights we conclude that an ideal inhibitor, if chemically realizable, would possess the electrostatic charge distribution of TMC125, so as to form strong interactions with the receptor, and the larger and more flexible substituents of TMC278, so as to minimize reorganization free energy penalties and the effects of resistance mutations, suitably modified, as in TMC125, so as to disfavor the formation of non-binding competent extended conformations when free in solution.

1. Introduction

The strength of the association between a ligand molecule and its target receptor is measured by the binding constant or, equivalently, by the standard free energy of binding. From a medicinal perspective there is great interest in the development of computational models capable of predicting accurately protein-ligand binding free energies [1,2]. A wide variety of methods have been developed to model the strength of protein-ligand association, spanning the field of QSAR knowledge-based approaches, to structure-based methodologies at various levels of theory ranging from empirical docking & scoring to quantum-mechanical descriptions [2]. This work is concerned a class of computational methodologies aimed at computing free energies of binding by employing high-level molecular mechanics descriptions of molecular interactions and classical descriptions of atomic motion [3-7]. Given a sufficiently accurate model of molecular interactions, these methods have the potential to incorporate details about the energetics and dynamics of the association equilibrium to address subtle aspects of drug development such as drug specificity and resistance.

Despite recent progress in physics-based binding free energy estimation methods, stemming from more accurate force field models, more extensive conformational sampling, improved free energy estimation methods, and faster computers, the conformational reorganization aspect of binding free energy calculations has received relatively little attention. The binding free energy is often the result of a large cancellation between the favorable work of forming receptor-ligand interactions and the unfavorable work to localize and reorganize the conformational ensembles of the ligand and receptor to their bound conformational states. The extent of this effect is likely widespread as for example it has been observed that that in protein-ligand complexes ligands assume energetically strained conformations [8,9]. While drug design is often concerned with strengthening receptor-ligand interactions, the reorganization component can play a fundamental role in regulating binding specificity in cases where variations of binding energies are expected to be small. In such cases optimization of binding affinity can proceed by strategies aimed at preorganizing the ligand for binding, that is by minimizing the unfavorable reorganization. For example this strategy has been used to optimize the presentation of viral epitopes in HIV vaccine candidates [10,11]. Examples exist of cases where optimization of a class of inhibitors was achieved by chemical rigidification of the ligands into their bound conformations [12-14]. Better correlation with experimental affinities have been reported when energetic scores are combined with ligand reorganization free energy estimates [15,16].

TMC125 (etravirine) [17] and TMC278 (rilpivirine) [18] (Figure 2) are two of the newest and most effective drugs for antiviral AIDS therapy. Rilpivirine’s discovery has been announced in 2005 [18] and has been approved for use by the Food and Drug Administration in May 2011. Both drugs act by inhibiting the function of the reverse transcriptase enzyme of the HIV virus (HIV-RT), which is essential for the initial transcription of viral RNA into DNA. TMC125 and TMC278 are members of the diarylpyrimidine (DAPY) class of non-nucleoside inhibitors (NNRTI) of HIV-RT. These bind to an allosteric pocket in the so-called palm domain of HIV-RT (Figure 1), causing a conformational change of the enzyme preventing it from properly processing the viral RNA.

The insurgence of drug-resistant HIV strains has been one of the most serious drawbacks of antiviral drug therapy, responsible for the growing inefficacy of the older

Figure 1. Ribbon representation of portion of the reverse transcriptase (RT) enzyme of the HIV-1 virus. The palm, fingers, thumb, and connection domains of the p66 subunit are shown in red, blue, green, and yellow respectively. The p51 subunit is in dark purple. The crystal structure of HIV-RT in complex with TMC125 is shown (PDB id 3MEC [33]). TMC125 is shown in ball and stick representation within the non-nucleoside binding pocket. The circle represents schematically the region of the receptor that is included in the computational model.

Figure 2. Chemical structures of TMC125 (etravirine) and TMC278 (rilpivirine). Hydrogen atoms are not shown. The nomenclature of the torsional angles are indicated for TMC278. The same nomenclature applies to TMC125, except for the angle which is specific to TMC278.

drugs and causing the rebound of the viral load in patients after prolonged exposure to drugs. DAPY inhibitors have been developed specifically to target wild-type HIV-RT as well as common mutant forms that are known to be resistant to other NNRTI’s, so as to be more broadly applicable for more prolonged therapy regimens. Detailed structural studies [19-21] have characterized the molecular mechanisms of development of resistance and the properties of these drugs that help them evade resistance mutations better than older drugs. The main conclusion is that DAPY compounds possess specific flexibility that allows them to reposition (“wiggle”) into the binding site and to change shape (“jiggle”) to avoid clashes with mutated residues and form alternative contacts to replace those lost to mutations. While TMC125 and TMC278 have some of the most favorable properties among DAPY inhibitors in this respect, TMC278 has been shown to be superior to TMC125 in terms of adaptability to mutant forms of the HIV-RT enzyme partly due to its greater flexibility, and partly because of its unique cyanovinyl functionality which binds into an hydrophobic tunnel lined by conserved residues and provides an additional and reliable anchor in spite of mutations capable of disrupting the action of most other inhibitors.

The intrinsic flexibility of these DAPY compounds make them ideal proving ground for the study of conformational reorganization effects in binding free energy calculations. The conformational propensities of TMC278 have been previously characterized by both modeling and experimental techniques [22-24]. The main conclusion is that in solution TMC278 preferentially assumes extended conformations very different from to the compact “U”- shaped conformation necessary for complexation with HIV-RT. Naturally, a conformational reorganization of TMC278 has to occur at a free energy cost disfavoring binding. In contrast, TMC125 is believed [25] to behave oppositely to TMC278 by incurring little conformational reorganization in going from the solution environment to the protein-bound state. This is because as it preferentially assumes the same “U”-shaped conformation in both conditions. We will show that the present modeling investigation confirms this hypothesis.

In this work we will address two main aspects of conformational reorganization. The first is understanding at a molecular level of how the opposite conformational propensities of the two drugs affect their affinity to HIV-RT. While direct measurements of binding constants are lacking, cell-based essays indicate that TMC278 has slightly higher affinity for wild-type HIV-RT than TMC125, despite that fact that reorganization effect are expected to favor the latter [21]. If so, the assumed larger reorganization free energy penalty for TMC278 must be somehow counterbalanced by stronger protein-ligand interactions so as to achieve higher or equivalent affinity than TMC125, for which reorganization is expected to be less significant.

The second aim of this work is to evaluate the role of conformational sampling [26,27] and ligand conformational reorganization in protein-ligand binding free energy calculation practices. We have recently illustrated [28,29] how an appropriate distribution of intermediate alchemical states is a necessary but not sufficient condition to ensure convergence of binding free energy calculations. Proper equilibration between unbound and bound conformational states, and between alternative binding modes, is also a critical requirement for obtaining reliable and reproducible results. The conformational sampling algorithm employed must be sufficiently powerful to allow frequent conformational transitions so as to rapidly explore conformational space with suitable frequency. While the Hamiltonian replica exchange (RE) algorithms we developed [30] are very useful to speed up equilibration along the alchemical thermodynamic path, the calculations reported here indicate that RE provides limited help in accelerating convergence when the kinetics is dominated by slow conformational reorganization processes. We will in particular examine the use of conformational restraints, devices that are often employed to speed up the convergence of binding free energy calculations [6,28,31,32]. Confirmation that results are independent of the simulation setup, including the nature of imposed conformational restraints, is an important aspect of validation tests of binding free energy protocols. In addition we will investigate the conditions under which conformational restraints accelerate convergence and the molecular mechanisms at the basis of these effects.

2. Methods

2.1. Molecular Systems

The chemical structures of the two DAPY NNRTI inhibitors of HIV-RT investigated are shown in Figure 2. In this figure we indicate the main torsional degrees of freedom of the DAPY compounds based on established nomenclature [21,23]. Non-nucleoside inhibitors bind to an allosteric pocket of HIV-RT in the so-called palm domain of the enzyme located between the thumb and fingers domains (see Figure 1). HIV-RT is a heterodimer of the p66 and p51 subunits; while the non-nucleoside pocket is lined by residues of both subunits, residues of p66 accounts of most of the interactions with the inhibitors. The structures of HIV-RT complexes and the characteristics of the interactions with DAPY inhibitors has been thoroughly described [19,21,33].

The initial structures of the complexes of HIV-RT in complex with TMC125 and TMC278 were obtained by the corresponding crystal structures (PDB id’s 3MEC [33] and 2ZD1 [21], respectively), hydrogen atoms were added using the Maestro program (Schrdinger, Inc.), and protonation states were assigned based on neutral pH. Harmonic restraints with force constant of 1 kcal/mol/Å2 were imposed on all heavy atoms of the receptor, allowing for a range of motion of approximately 3 Å at room temperature. The structural models were then subjected to energy minimization before further processing. The calculations employed a reduced representation of the receptor which includes any residue of the receptor with atoms less than 12 Å from any atom of the ligand, resulting in a region surrounding the ligand of approximately 36 Å in diameter (Figure 1). The reduced model was adopted only for reasons of computational convenience; the 12 Å cutoff on non-bonded interactions and the rigidity of the receptor, make model predictions independent from atomic interactions beyond the modelled region. The model of the receptor includes 114 residues (residues 88 - 112, 171 - 195, 220 - 243, 314 - 323, 347 - 350, and 378 - 385 of the p66 subunit and residues 132 - 142 of the p51 subunit) and 1905 atoms. Residues at chain termini were capped using N-methyl-amine (NMA) and acetyl (ACE) groups, placed based on the coordinates of the backbone atoms of the corresponding deleted residues.

2.2. BEDAM Binding Free Energy Computational Protocol

Absolute binding free energies of the complexes with HIV-RT with TMC125 and TMC278 were computed using the Binding Energy Distribution Analysis Method method (BEDAM) [30]. Briefly, BEDAM computes the standard binding free energy between a receptor and a ligand with implicit solvation by means of the expression

(1)

which follows, without approximations, from a wellestablished statistical mechanics theory of molecular association [34]. In Equation (1), is the standard concentration of ligand molecules (set to M, or equivalently 1668 Å–3), is the volume of the binding site (see below), and is the excess free energy of binding defined as the difference in free energy between the coupled state of the complex, in which the ligand and the receptor are fully interacting, and the decoupled state in which, while the ligand is confined in the binding site, the receptor and the ligand are interacting only with the solvent continuum and not with each other.

Central to the method is the concept of a binding energy function

(2)

defined for each conformation of the complex as the difference between the effective potential energies and of, respectively, the ligand-bound and ligand-separated conformations of the complex without internal conformational rearrangements. BEDAM is based on a λ-dependent hybrid potential energy function of the form

(3)

where

(4)

is the potential energy of the uncoupled state of the complex. It is easy to see from Equations (2)-(4) that and correspond to the effective potential energy functions of the coupled and decoupled states of the complex, respectively. The excess binding free energy is by definition the difference in free energy between the and states. As described below, is obtained from statistical analysis of values of the binding energy [Equation (2)] sampled from these states as well intermediate states with values of between these two extremes.

Rather than simulating each λ state independently, BEDAM employs a Hamiltonian replica exchange (HREM) λ-hopping strategy whereby simulation replicas periodically attempt to exchange λ values through Monte Carlo (MC) λ-swapping moves. λ-exchanges are accepted with the Metropolis probability [30] where is the difference in λ’s being exchanged and is the difference in binding energies between the replicas exchanging them. Replica exchange strategies of this kind yield superior conformational sampling and more rapid convergence rates by allowing conformational transitions to occur at the value of at which they are most likely to occur and to be then propagated to other states [26].

To improve convergence of the free energy near, in this work we employ a modified “soft-core” binding energy function of the form

(5)

where is some large positive value (set in this work as 106 kcal/mol). This modified binding energy function, which is used in place of the actual binding energy function [Equation (2)] wherever it appears, caps the maximum value of the binding energy while leaving unchanged the value of favorable binding energies [29, 35].

In this work we employed multistate Bennett acceptance ratio estimator (MBAR) [36] to estimate the excess free energy of binding, using as input binding energy samples obtained from the HREM simulations. Binding free energies are computed directly from the MBAR dimensionless free energies using the relationship

(6)

The MBAR dimensionless free energies, , are defined as the negative of the logarithm of the -dependent biased partition functions. In this case the dimensionless free energies are estimated by the self-consistent solution of the set of equations [36]

(7)

where is the th binding energy sample from replica, is the number of replicas and is the total number of binding energy samples from replica. For the MBAR analysis we employed the code provided by John Chodera and Michael Shirts (http://alchemistry. org). Statistical uncertainties were obtained from the standard deviation of the binding free energies from the last four 1 ns blocks of binding energy data (see below).

For later use, the reorganization free energy for binding is defined by the expression [37]

(8)

where is the average binding energy at and is the standard binding free energy from Equation (1). The former is computed as described below and is computed by difference using Equation (8).

2.3. Binding Free Energy Calculations with Conformational Restraints

The concept of the reorganization free energy is useful to formalize situations in which, by imposing suitable conformational restraints, auxiliary states are introduced which do not match exactly with the end point states of the complex [32,37]. Consider for example Figure 3, in which the excess free energy of binding is first decomposed in the reorganization free energy of restraining the ensembles of conformations of the receptor and the ligand in solution to chosen macrostates and. For instance, in the application below the macrostate is defined as the ensemble of conformations of the ligand with torsional angles within a specified range. The free energy for this process is related to the population , defined as the probability of finding a conformation belonging to the specified macrostate, in absence of restraints and ligand-receptor interactions ():

(9)

Following this step, we consider the binding free energy, , between the and species, that is the excess binding free energy when the receptor and the ligand are limited to the chosen macrostates. In general the resulting state of the complex, denoted by ()in Figure 3, does not match the full complex state because in the former the receptor and the ligand are limited to their respective macrostates., defined similarly to, is the free energy difference between these two states. In the present application we assume because the restrained macrostate encompasses all of the conformational ensemble of the complex; that is the and species are virtually equivalent.

Based on the above we have

(10)

Equation (10) expresses the notion that a binding free energy calculation can be conducted with or without conformational restraints, and, if the binding free energy, , obtained with restraints is properly combined with the free energy, , of imposing the restraints, either approach should return the same binding free energy. Hence the use of restraints is purely a matter of efficiency; restraints should be considered whenever

Figure 3. Thermodynamic cycle illustrating the restrainand-release decomposition of the excess free energy of binding [Equation (10)]. Although not indicated, the ligand here is assumed to be always sequestered in the binding site. R and L represent the free receptor and ligand, R* and L* represent the receptor and ligand restrained within a conformational macrostate, (RL)* represents the complex in which receptor and ligand are restrained within their macrostates, and RL represents the free complex.

they can lead to faster convergence of binding free energy estimates.

2.4. Computational Details

We employed the OPLS-AA/AGBNP2 effective potential in which the effect of the solvent is represented implicitly by means of the AGBNP2 implicit solvent model [38,39] together with the OPLS-AA [40, 41] force field for covalent and non-bonded interatomic interactions. Parallel molecular dynamics simulations were conducted with the IMPACT program [42]. The simulation temperature was set to 300K. We employed 24 intermediate steps at 0, , , , , , , , , 0.01, 0.05, 0.075, 0.1, 0.175, 0.25, 0.3125, 0.375, 0.4375, 0.5, 0.58, 0.67, 0.75, 0.875 and 1. The binding site volume was defined as any conformation in which the middle nitrogen atom (the one directly in between the two sidechains) of the central pyridine ring of the ligand is within 2 Å of the center of mass of the C atoms of residues 97, 102, 179, and 181 of the p66 subunit of HIV-RT. The guest was sequestered within this binding site volume by means of a flat-bottom harmonic potential. Based on this definition, the volume of the binding site, , is calculated as approximately Å3 and in Equation (1) is kcal/mol. BEDAM parallel -hopping calculations were performed for 8 ns of molecular dynamics per replica (192 ns of molecular dynamics for each of four binding free energy calculations with 24 replicas each, resulting in a total of 768 ns of simulations). The second half of the trajectory (4 ns per replica) was used for data analysis. Binding energies were sampled with a frequency of 1 ps. Statistical uncertainties are evaluated as the standard error of the distributions of estimates obtained from each 1 ns section of the trajectories.

Probability distributions and averages are computed by weighted averages over the trajectories of all of the replicas [36]. Briefly, each sample from replica is assigned a -dependent thermodynamic weight given by

(11)

where the are the solution of the MBAR Equation (7) and the other symbols are as for Equation (7). Thermodynamic averages are formulated as weighted averages over the observations using the weights defined by Equation (11). For example the average binding energy at is formulated as

(12)

and the probability density, , of a structural quantity (a torsional angle for instance) at is expressed as

(13)

where is the value of the structural quantity of sample and is a function defined as 1 if the argument falls within the bin centered at with width and zero otherwise.

BEDAM binding free energy calculations were conducted for TMC125 and TMC278 with and without torsional angle restraining potentials designed to restrict the ligand within the bound ligand macrostate, defined as any conformation with and torsional angles within a interval around (see Figures 1 and 3). Restraining potentials were implemented as flat-bottom harmonic restraints on the and torsional angles centered on with a tolerance of on either side. The population of the bound macrostate at, required to compute from Equation (9), was evaluated by numerical integration of the joint probability density, , of the and angles at within the region. probability densities are displayed in terms of the corresponding potential of mean force (PMF)

(14)

taking the smallest value of the PMF as the zero of energies.

3. Results

3.1. Conformational Analysis of Free and Bound TMC125 and TMC278

Since they result in the largest conformational variations (see Figure 4), we will focus first on the and torsional degrees of freedom of the ligands as defined in Figure 2. The computed conformational distributions for TMC125 and TMC278 free in solution with respect to these torsional angles (see Methods) are shown in Figures 5 and 6, respectively. As it can be seen from Figure 5, TMC125 assumes predominantly U conformations in solution with a minor L1 component. Remarkably, TMC278 displays the opposite behavior by occupying mostly extended (E) and L2 conformations (Figure 6), which are virtually never observed for TMC125. Of the two ligands, TMC278 can be considered the more flexible compound in solution since, although less frequently than TMC125, it also visits U and L1 conformations; so overall TMC278 is capable of adopting all of the main conformations available to DAPY compounds.

By numerically integrating the conformational distributions in Figures 5 and 6 within the U conformational state (indicated by the square central region of the diagrams), we obtain the restraining free energies

defined by Equation (9) and reported in Table 1. The restraining free energy for TMC125 (kcal/mol) is much smaller than the one for TMC278 (kcal/mol) reflecting the relative abundance of U conformations for

Figure 4. The four main conformational states of the DAPY compounds based on the τ3 and τ4 torsional angles, illustrated for TMC278. The U conformational state is centered around the “U”-shaped conformation with τ3 = τ4 = 0˚. The U state corresponds to the bound conformation. The L1 and L2 conformational states are centered around the “L”- shaped conformations at (τ3, τ4) = (180˚, 0) and (τ3, τ4) = (0˚, 180˚), respectively. Finally, the E or extended conformational state is centered around (τ3, τ4) = (180˚, 180˚).

Figure 5. The distribution of conformations of free TMC125 in solution (λ = 0) expressed in terms of the potential of mean force (in kcal/mol) with respect to the τ3 and τ4 torsional angles (see Figure 1). The regions corresponding to the main conformational states (see Figure 3) are indicated.

Figure 6. The distribution of conformations of free TMC278 in solution (λ = 0) expressed in terms of the potential of mean force (in kcal/mol) with respect to the τ3 and τ4 torsional angles (see Figure 1). The regions corresponding to the main conformational states (see Figure 3) are indicated.

the two molecules in solution. TMC125 has a natural propensity of forming U conformations so restraints that force it to adopt that state have little influence on its thermodynamics. Conversely, restraints have a large effect on the thermodynamics of TMC278 because in absence of restraints this preferentially adopts other conformations.

Despite the large differences in conformational propensities when free in solution, the two compounds have similar conformational distributions in the bound state. The interactions with the protein receptor force each ligand to adopt conformations well contained within the U conformational macrostate. The thermodynamic analysis above indicates that TMC278 suffers a substantial reorganization free energy penalty associated with the reorganization of the and degrees of freedom in order to adapt to the receptor. This is confirmed by the significantly smaller reorganization free energy (kcal/mol; Table 2, 5th column) obtained for TMC278 binding with torsional restraints on the and angles, forcing it to remain withing the U state, compared to unrestrained conditions (kcal/mol). In contrast, the computed reorganization free energies for TMC125 with or without torsional restraints are very similar to each other (Table 2, 5th column) indicating minor reorganization of the distributions of and angles.

3.2. Binding Free Energy Estimates

The computed binding free energies for the complexes of wild-type HIV-RT with the TMC125 and TMC278 inhibitors with and without torsional restraints are shown in Table 1. Based on Equations (1) and (9), the binding free energy values, ,

Table 1. Experimental and computed standard binding free energies of the complexes of HIV-RT with the TMC125 and TMC278 ligands with and without torsional restraints.

Table 2. Decomposition of computed binding free energies into average binding energy and reorganization free energy components.

computed with torsional restraints (4th column) are combined with the free energies, , of imposing the restraints (5th column) to obtain the binding free energy estimates (6th column), which can be compared to the experimental measurements. The values obtained from the unrestrained simulations yield directly the actual binding free energy estimates without restraining corrections. The binding free energy estimates, , from the unrestrained and restrained simulations are in agreement within statistical uncertainty for both TMC125 and TMC278, thereby providing numerical validation to the theoretical requirement that binding free energy estimates should not depend on the calculation settings. Indeed, as discussed above, the binding free energy calculations conducted with conformational restraints, when appropriately combined with the free energy of imposing the restraints, should return equivalent results to those of binding free energy calculations without restraints.

The binding free energy estimates (6th column in Table 1) for TMC278 are in good agreement with the experimental measurement (2nd column). The computed binding free energy values for TMC125 are instead 1 to 2 kcal/mol more favorable to binding than the experimental value, causing the affinity ranking of the two compounds to be reversed relative to the experiments. Overall the agreement between calculations and experiments is reasonable considering the approximations in the potential energy model and the fact that we have not attempted to model the influence of receptor reorganization. The present model also neglects the influence of protonation and aggregation equilibria that are known to be important in these systems [23]. Furthermore, it should be noted that direct binding free energy measurements and inhibition constants are not available for these complexes. The experimental values listed in Table 1 have been derived from EC values from cell-based essays [18,43] which are influenced by many characteristics of the drugs (propensity for aggregation, membrane transport, non-specific binding, degradation, etc.) [22] and represent only qualitatively their thermodynamic affinity for the HIV-RT enzyme, which is the only quantity modeled here.

4. Discussion

4.1. Use of Conformational Restraints in Binding Free Energy Calculations

Conformational restraints, often in the form of additional potential terms which restrict the position and conformation of the ligand and the receptor, are commonly used in binding free energy calculations. One use of restraints derives from the need to define a binding site volume which specifies the complexed state and introduces the effect of ligand concentration [34]. This issue and its in fluence on binding free energy estimates has been discussed [30,37,44]. Conformational restraints are, in addition, often used as computational devices to speed up the convergence of binding free energy calculations by reducing the size of conformational space of the complex that needs to be sampled to achieve convergence [6, 28,31] When used in the latter context, it is important that the inclusion of restraints affects only the efficiency of the calculations and not their outcomes. As discussed in Section 2, when appropriate corrections are applied representing the work corresponding to imposing and releasing the restraints, equivalent results should be obtained without restraints or with different restraining schemes.

4.1.1. Independence of Binding Free Energy Estimates

Protein-ligand binding free energy protocols are rarely subjected to validation tests with respect to variations of restraints. In this work we have numerically confirmed that the BEDAM protocol [30] combined with multi-state estimation of conformational free energies [24] returns consistent binding free energy estimates for two complexes of HIV-RT with DAPY inhibitors with or without the presence of ligand conformational restraints. Conformational restraints on key torsional degrees of freedom have been optionally used here to force the ligands to remain within the bound conformational macrostate even when free in solution, and all along the thermodynamic path. Restraints have very minor thermodynamic effects for TMC125, which preferentially adopts bound conformations in solution without the help of restraints. Restraints have instead a large effect on the affinity of TMC278 because this does not often visit bound conformations unless restraints are applied. The computed binding free energy of restrained TMC278 (kcal/mol) is significantly larger in magnitude than the binding free energy of unrestrained TMC278 (kcal/mol) because the latter includes a larger portion of the reorganization free energy penalty for reorganizing the ligand conformational ensemble into the bound macrostate. When, however, the work of imposing the restraints (estimated at kcal/mol) is considered, the two approaches yield equivalent binding free energy estimates within statistical uncertainty.

As it can be clearly seen from the data in Table 1, restraints have very little influence on average binding energy values, which measure the strength of ligand-receptor interactions. This is expected because in this case the bound conformational macrostate is well contained within the range allowed by the restraints. The influence of restraints is strictly limited to reorganization free energy penalties. The computed reorganization free energy of TMC278 is approximately kcal/mol with restraints compared to kcal/mol without them. The larger reorganization by kcal/mol without restraints is clearly due to the reorganization free of the ligand, estimated at approximately kcal/mol, due to the reorganization from the unrestrained ensemble in solution (Figure 6) to the ensemble allowed by the restraints (the U macrostate in Figure 6).

4.1.2. Effects of Restraints on the Rate of Convergence

As discussed above increasing computational efficiency is the main motivation for using conformational restraints in binding free energy calculations. One measure of computational efficiency is the simulation time necessary to achieve a given statistical uncertainty of the free energy estimate or, equivalently, the magnitude of the statistical uncertainty for a given simulation time, smaller being better. One puzzling result of the calculations presented here is that the use of restraints did not lead to smaller statistical uncertainties (see Table 1). The estimated uncertainty for the binding free energy of TMC125 with restraints is significantly larger than for the unrestrained calculation (vs kcal/mol, respectively). The statistical uncertainties for the two calculations for TMC278 binding are very similar (both are approximately kcal/mol) with the uncertainty (kcal/mol, see Table 1) of the restraining free energy, , contributing somewhat to the overall uncertainty of the free energy estimate with restraints.

Some insights on this issue can be gained by examining how binding free energy estimates vary as the simulation progresses (Figure 7). In this figure we plot binding free energy estimates from data from 1 ns sections of the simulation trajectories. Note that these are not cumulative results. We see that for each calculation there is a relatively short equilibration time (the results presented in Table 1 are from the second half of the trajectory); furthermore the free energy estimates for TMC125 show significantly smaller variations than for TMC278. The unrestrained TMC125 calculation shows the smallest variations reflecting its small statistical uncertainty. The restrained TMC125 calculation displays larger variations. There are no obvious differences in the variations of the free energy estimates for TMC278 with or without restraints, however the variations for TMC278 are consistently greater than those for TMC125.

In prior studies [28,29] we have identified as one major cause of slow convergence the long timescales to achieve conformational transitions between unbound conformations with large binding energies to bound conformations with low binding energies. In Figures 8 and 9 we show the time trajectories of the binding energies of some of the replicas of the HREM simulations of the complexes with TMC125 and TMC278 with and without torsional restraints. In general, a larger number of transitions from high binding energies to low binding energies and vice versa is indicative of better conformational

Figure 7. Binding free energy estimates from 1 ns sections of trajectories as a function of simulation time. Note that these are not cumulative estimates.

(a)(b)

Figure 8. Time trajectories of the binding energies of some of the replicas of the BEDAM simulation of TMC125 in complex with HIV-RT without (a) and with (b) torsional restraints.

(a)(b)

Figure 9. Time trajectories of the binding energies of some of the replicas of the BEDAM simulation of TMC278 in complex with HIV-RT without (a) and with (b) torsional restraints. Fewer long-lived binding/unbinding transitions occur compared to TMC125 (Figure 8). Note that without torsional restraints (a) the region of binding energies near –20 kcal/mol is rarely crossed. Also note that longer trajectories are shown for unrestrained TMC278 (a).

sampling and, in turn, faster convergence.

Based on this criterion, of the four calculations, the unrestrained TMC278 simulation appears to undergo the fewest conformational transitions (Figure 9(a)); several of the replicas remain in either bound conformations or unbound conformation for the entire duration of the simulation and never undergo transitions, while few replicas undergo one or two transitions. The corresponding simulation for TMC278 with restraints (Figure 9(b)) shows better behavior; all of the replicas display rapid transitions between low binding energies (around kcal/mol) and high binding energies (-kcal/mol). More rarely and similarly to the unrestrained simulation, long-lived transition to high binding energies occur (see for example bottom-center panels in Figures 9). The simulations for TMC125 binding (Figure 8) display analogous behavior except that the kinetics of the unrestrained and restrained simulations are more similar to each other and undergo more long-lived transitions than the TMC278 simulations.

Structural analysis shows that low-to-high and highto-low rapid transitions in binding energy space observed in the TMC125 simulations and the restrained TMC278 simulation (Figures 8 and 9(b)) correspond to excursions in which replicas rapidly traverse -space from to and concomitantly assume high binding energies, but with the ligand molecule remaining in the neighborhood of the bound orientation and getting quickly “recaptured” in the bound state at low binding energies. The infrequent long-lived transitions are cases in which instead during such excursions the ligand leaves the bound orientation and starts to freely diffuse in orientational space before eventually finding again the bound conformation. Compared to TMC125, the TMC278 ligand spends more time exploring orientations before finding an orientation suitable for binding. This is due in part of the presence of the extra orientational constraint due to the positioning of the cyanovinyl group in its receptor cavity and, for the unrestrained simulation, also because of the tendency of TMC278 to migrate to the E, L1, and L2 rotameric states, which are incompatible to the formation of strong interactions with the receptor.

It appears therefore that orientational degrees of freedom play an important role in the kinetics of binding and unbinding in these systems and that they strongly influence the rate of convergence. Our conclusion is that the convergence rate is for the most part dictated by slow orientational diffusion. An equivalent conclusion is that slow convergence is due to the difficulty of estimating precisely the size of the region of orientational space compatible with binding relative to the size of the overall orientational space. We predict therefore that imposing restraints on the ligand orientation would speed up convergence of binding free energy considerably [28]. On the other hand restraints must be designed based on prior structural knowledge so at to not interfere with the bound conformational state (or calculation of must be be considered, see above), hence they are not easily applicable to situations in which structural information about the complex is uncertain or when multiple binding modes are possible.

4.1.3. Can Restraints Slow Down Convergence?

One puzzling result of this work is that the statistical uncertainty for the binding free energy of TMC125 estimated from the restrained simulation is five times larger than that from the unrestrained simulation (see Table 1). Restraints are believed to generally aid convergence of binding free energy calculations [6] rather than worsen it, as in this case. We found indications that the poorer convergence of the restrained simulation of TMC125 binding is caused by the reduced number of binding pathways available to the ligand compared to the unrestrained simulation.

Indication of this phenomenon is apparent from the analysis of ligand conformational distributions at intermediate values of λ from the unrestrained simulation (Figures 10 and 11). As shown in Figure 5, TMC125 assumes mostly U and L1 conformations when dissociated from the receptor and exclusively U conformations when it is fully associated. However, the shift from U and L1 conformations to exclusively U conformations does not appear to occur as a direct population transfer from L1 to U. At intermediate values of the L2 and E conformational states become also significantly populated (Figure 10). This occurs mainly by population transfer from the L1 macrostate to the L2 and E macrostates (Figure 11); with increasing, the populations of the L2 and E states increase rapidly at the expense of the L1 population. At larger values of’s, population transfer from the L2 and E states to the U state then occurs (Figure 11). This observation is clearly indicative of the existence of binding routes whereby TMC125, initially in a L1 conformation, reach the bound U state by going through E and/or L2 intermediate states. These binding mechanisms are absent in the restrained simulation in which, by construction, only U conformations are allowed.

We therefore suggest that in this case conformational restraints oppose convergence rather than aid it because they block additional binding pathways and decreasing the chances of more rapid interconversion between unbound and bound conformations of TMC125. This conclusion goes contrary to the generally accepted wisdom that conformational restraints always speed up convergence by reducing the size of conformational space that needs to be sampled. This example shows that reducing the accessible conformational space may have the side effect of precluding some conformational transition pathways and slow down convergence.

4.1.4. Order-Disorder Pseudo-Phase Transition

It is apparent from Figures 8 and 9 that the unrestrained simulation of TMC278 (Figure 9(a)) behaves differently from the others in terms of binding/unbinding transitions. In this simulation a rarely crossed barrier exists in binding energy space separating unbound and bound conformations. The location of this barrier is clearly seen in the

Figure 10. The distribution of conformations of TMC125 in complex with HIV-RT at λ = 0.005 expressed in terms of the potential of mean force (in kcal/mol) with respect to the τ3 and τ4 torsional angles (see Figure 1). At this intermediate value of λ TMC125 visits all of the main rotamers (see Figure 3).

Figure 11. Conformational populations of TMC125 as a function of λ from the unrestrained binding simulation. At λ = 0 only the U and L1 conformations are significantly populated. As λ increases the L1 population decreases rapidly concomitantly with a rapid increase of the populations of the L2 and E populations. At even larger values of λ only the U (bound) conformation is significantly populated.

binding energy distributions of this system (Figure 11). For binding energies are predominantly above kcal/mol, whereas for much lower values of binding energies are sampled. Values of binding energies near kcal/mol are rarely visited at any and binding energy distributions straddling this region are bimodal (see distributions in Figure 12 at and 0.3125).

Figure 12. Probability densities of the binding energy, u, for the simulation of unrestrained TMC278 complexed with HIV-RT as a function of λ near the critical state. The distributions near the critical λ () are bimodal and the intermediate region of binding energies around –20 kcal/mol is infrequently visited at any λ.

This behavior is the tell-tale indication of an order-disorder pseudo-phase transition [45]; a situation in which a macrostate at high energy and high entropy is in equilibrium with another macrostate at low energy and low entropy. In this specific case the two macrostates correspond to the L1, L2, and E rotamers of TMC278, which are not capable of binding strongly to HIV-RT, and the U rotamer that becomes predominant when interactions with the receptor are turned on. The switch to U conformations occurs rapidly as a function of (Figure 13), a consequence of the large difference in average binding energies between bound and unbound conformations [28]. Pseudo-phase transitions of similar strength are not observed in the simulations of TMC125 and the simulation of TMC278 with restraints because of reduced conformational space and entropy of the unbound macrostate in these cases. Order-disorder transitions are also not observed for unrestrained TMC125 because of its relative rigidity compared to TMC278.

As previously observed [28], the rarely visited intermediate region of binding energies (Figure 12) constitutes an obstacle to binding/unbinding transitions, thereby hampering convergence of binding free energy calculations. Advanced conformational sampling tools are being developed in this and other research fields to address conformational sampling challenges of this kind [45,46].

4.2. Insights on Thermodynamic Driving Forces for Binding

Our original expectation was that, due to its higher flexi-

Figure 13. Conformational populations of the U, L1, L2 and E rotamers of unrestrained TMC278 in complex with HIV-RT as a function of the coupling parameter λ. The transition from rotamer states prevalent near the decoupled state (λ = 0) to the U rotamer occurs in a relatively small region of λ. As we have previously shown [28], the gradient of the population profile for the U state at 50% population is proportional to the difference in average binding energy between the bound and unbound macrostates of the complex.

bility in solution, binding of TMC278 would be characterized by a larger reorganization free energy penalty than TMC125. The calculations however indicate the opposite. As it can be seen by the values in Table 2 the binding reorganization free energy of TMC125 is approximately kcal/mol; kcal/mol greater than the reorganization of TMC278. This is despite the fact that reorganization of the and torsional angles alone within the defined ranges (see above) opposes binding of TMC278 by approximately kcal/mol (see in Table 1, 5th column). Clearly, collective reorganization of other degrees of freedom of TMC125 must be contributing more to the overall reorganization free energy than the and angles which account for only approximately kcal/mol of. Hints concerning the additional possible sources of reorganization come from a comparative analysis of the bound conformational ensembles of TMC125 and TMC278.

The motion of the DAPY compounds within the non-nucleoside binding site of HIV-RT has been described as “wiggling”, that is internal conformational rearrangements such as rotations of sidechains, and “jiggling”, or overall translation and rotational motion of the molecule. As illustrated in Figure 14, analysis of conformational ensembles of the complexes at reveals that, when bound to HIV-RT, TMC278 is significantly more mobile than TMC125 in both respects. Figure 14(a) represents the distribution of and torsional angles characterizing the “wiggling” motion of the left sidechain of the ligands in Figure 2, which contains the oxy-dimethylbenzonitrile moiety of TMC125 or the

(a)(b)

Figure 14. (a) Distribution of the τ1 and τ3 torsional angles of TMC125 (blue triangles) and TMC278 (red squares). The τ1 and τ3 torsional angles probe the motion of the ligand sidechain which contains the oxy-dimethylbenzonitrile moiety of TMC125 or the amino-dimethylbenzocyanovinyl moiety of TMC278 (“wiggling” motion); (b) Probability density (in Å–1) of displacement δd of the N1 atom of the central pyrimidine ring (the nitrogen atom in between the two substituents) relative to the its position in the corresponding crystal structure (“jiggling” motion).

amino-dimethylbenzocyanovinyl moiety of TMC278. This data clearly indicates that the sidechain of TMC278 undergoes a wider range of libration motion than TMC125, despite the fact that the cyanovinyl group specific to TMC278 is known to be anchored in what has been a described as a hydrophobic tunnel of the receptor [21]. In comparison (data not shown), we observed little difference in the distribution of and torsional angles of the amino-benzonitrile group, which is the sidechain common to both compounds (see Figure 2).

We have also observed a larger amount of “jiggling” motion of TMC278 relative to TMC125. As an illustration, Figure 14(b) shows the distribution of distances of the position of the N1 atom of the central pyrimidine ring (the nitrogen atom in between the two substituents) relative to the its position in the corresponding crystal structure. The variance of the distance distribution for TMC278 is 60% greater than that for TMC125, underscoring greater translational motion, or “jiggling”, of TMC278. Collectively these observations indicate that within the bound conformational macrostate the motion of TMC125 is more restricted than TMC278 thereby explaining in part the greater reorganization free energy penalty for binding of the former. A large contribution to reorganization may also originate from librational and vibrational degrees of freedom that we have not analyzed [47].

The larger unfavorable reorganization free energy for TMC125 is counterbalanced by a more favorable average binding energy; the latter is approximately kcal/mol more negative than for TMC278 (see Table 2, 4th column), compared to a deficit of roughly kcal/mol for the reorganization free energy component (Table 2, 6th column), resulting in a net ~1.5 kcal/mol advantage for TMC125 in terms of binding free energy. Among the energetic contributions to the binding energy (ligandreceptor electrostatic and van der Waals interactions, and desolvation free energy), stronger receptor-ligand electrostatic interactions accounts for the majority of the difference in average binding energy in favor of binding of TMC125.

Given its weaker interactions with the receptor (see Table 2), the affinity of TMC278 for HIV-RT would have been even weaker relative to TMC125 if not for its smaller (and unexpected) reorganization free energy. The present results indicate that the influence of the cyanovinyl group on the thermodynamics of binding of TMC278 can be ascribed less to increased strength of receptor-ligand interactions and more to increases redundancy and flexibility of these interactions, so that they are maintained in a wider range of conformations. This hypothesis is supported by the less favorable average binding energy of TMC278 relative to TMC125 (Table 2) and the lower reorganization free energy penalty for binding of TMC278, due in part to the greater freedom of motion of the sidechain containing the cyanovinyl group (Figure 14(a)). Based on structural studies of complexes of TMC278 with mutant forms of HIV-RT, Das et al. [21]. have observed that the ability of the cyanovinyl group to reposition and switch contacts between receptor residues is key to the greater potency of TMC278 against a wide range of drug-resistant mutants of HIV-RT.

Structure-based drug optimization strategies are often aimed at increasing the number and strength of ligandreceptor interactions. This example shows that both greater affinity and improved resistance profiles can be achieved by focusing instead on increasing the redundancy and flexibility of ligand-receptor interactions so that they are maintained over a wider range of conformations of the complex.

5. Conclusions

We have investigated some aspects of the complexation equilibrium between HIV-RT and TMC125 (etravirine) and TMC278 (rilpivirine), two antiviral drugs recently approved for the treatment of AIDS. We have confirmed that the BEDAM computational protocol we employed yields consistent results regardless of applied conformational restraining potentials on torsional degrees of freedom. The presence of the restraints is shown however to have significant influence on the convergence rate of the simulations and the conformational transition mechanisms. Unexpectedly, torsional restraints do not enhance convergence for these systems. Torsional restraints on the ligands slow down convergence by precluding some pathways for binding/unbinding transitions and moreover equilibration along ligand orientational degrees of freedom, which are not limited by the applied restraints, seem to dominate the rate of binding/unbinding transitions and ultimately the rate of convergence of binding free energies. We have also shown that restraints have little influence on the ligand-receptor interaction energy driving force for binding, while they lower the reorganization free energy penalty.

Physical insights concerning the thermodynamic driving forces for or against binding have been obtained. Contrary to expectations TMC125 forms stronger interactions with the receptor than TMC278 and, despite the higher flexibility of the latter when dissociated by the receptor, TMC125 suffers a higher reorganization free energy penalty for binding than TMC278. This is rationalized by the wider range of conformations of the TMC278/HIV-RT complex partly due to the ability of the cyanovinyl substituent of TMC278 to form interactions with the receptor in multiple orientations. Mindful of the approximations in the computational model, and in particular to the lack of modelling of receptor reorganization effects, this computational study proposes ligand design principles to increase affinity to HIV-RT and possibly to achieve more favorable inhibition profiles across mutant forms of the receptor. An ideal inhibitor, if chemically realizable, would possess the electrostatic charge distribution of TMC125 to form strong interactions with the receptor and the large and flexible substituents of TMC278 to minimize reorganization free energy penalties and the effects of resistance mutations, suitably modified as in TMC125 so as to disfavor the formation of non-binding competent extended conformations (L1, L2 and E) in solution.

6. Acknowledgements

The author expresses sincere gratitude to Dr. Ronald Levy for his contributions to the development of the BEDAM method and stimulating collaborative work on the statistical thermodynamics of protein-ligand binding. This work has been supported by grants by the National Institute of Health (GM30580) and the National Science Foundation (CDI type II no. 1125332). The calculations reported in this work have been performed at the BioMaPS High Performance Computing Center at Rutgers University funded in part by the NIH shared instrumentation grants no. 1 S10 RR022375 and 1 S10 RR027444, and on the Trestles cluster at the San Diego supercomputing center under XSEDE National Science Foundation allocation grant no. TG-MCB100145.

REFERENCES

  1. W. L. Jorgensen, “The Many Roles of Computation in Drug Discovery,” Science, Vol. 303, No. 5665, 2004, pp. 1813-1818. doi:10.1126/science.1096361
  2. O. Guvench and A. D. MacKerell, “Computational Evaluation of Protein-Small Molecule Binding,” Current Opinion in Structural Biology, Vol. 19, No. 1, 2009, pp. 56-61. doi:10.1016/j.sbi.2008.11.009
  3. M. K. Gilson and H.-X. Zhou, “Calculation of Protein-Ligand Binding Affinities,” Annual Review of Biophysics and Biomolecular Structure, Vol. 36, 2007, pp. 21-42. doi:10.1146/annurev.biophys.36.040306.132550
  4. M. R. Shirts, D. L. Mobley and J. D. Chodera, “Alchemical Free Energy Calculations: Ready for Prime Time?” Annual Reports in Computational Chemistry, Vol. 3, 2007, pp. 41-59. doi:10.1016/S1574-1400(07)03004-6
  5. D. L. Mobley and K. A. Dil, “Binding of Small-Molecule Ligands to Proteins: ‘What You See’ Is Not Always ‘What You Get’,” Structure, Vol. 17, No. 4, 2009, pp. 489-498.
  6. Y. Q. Deng and B. Roux, “Computations of Standard Binding Free Energies with Molecular Dynamics Simulations,” The Journal of Physical Chemistry B, Vol. 113, No. 8, 2009, pp. 2234-2246. doi:10.1021/jp807701h
  7. J. D. Chodera, D. L. Mobley, M. R. Shirts, R. W. Dixon, K. Branson and V. S. Pande, “Alchemical Free Energy Methods for Drug Discovery: Progress and Challenges,” Current Opinion in Structural Biology, Vol. 21, No. 2, 2011, pp. 150-160.
  8. E. Perola and P. S. Charifson, “Conformational Analysis of Drug-Like Molecules Bound to Proteins: An Extensive Study of Ligand Reorganization upon Binding,” Journal of Medicinal Chemistry, Vol. 47, No. 10, 2004, pp. 2499-2510. doi:10.1021/jm030563w
  9. J. Liebeschuetz, J. Hennemann, T. Olsson and C. Groom, “The Good, the Bad and the Twisted: A Survey of Ligand Geometry in Protein Crystal Structures,” Journal of Computer-Aided Molecular Design, Vol. 26, No. 2, 2012, pp. 169-183.
  10. M. Lapelosa, E. Gallicchio, G. F. Arnold, E. Arnold and R. M. Levy, “In Silico Vaccine Design Based on Molecular Simulations of Rhinovirus Chimeras Presenting HIV-1 GP41 Epitopes,” Journal of Molecular Biology, Vol. 385, No. 2, 2009, pp. 675-691. doi:10.1016/j.jmb.2008.10.089
  11. M. Lapelosa, G. F. Arnold, E. Gallicchio, E. Arnold and R. M. Levy, “Antigenic Characteristics of Rhinovirus Chimeras Designed in Silico for Enhanced Presentation of HIV-1 GP41 Epitopes,” Journal of Molecular Biology, Vol. 397, No. 3, 2010, pp. 752-766. doi:10.1016/j.jmb.2010.01.064
  12. J. E. DeLorbe, J. H. Clements, M. G. Teresk, A. P. Benfield, H. R. Plake, L. E. Millspaugh and S. F. Martin, “Thermodynamic and Structural Effects of Conformational Constraints in Protein-Ligand Interactions. Entropic Paradoxy Associated with Ligand Preorganization,” Journal of the American Chemical Society, Vol. 131, No. 46, 2009, pp. 16758-16770. doi:10.1021/ja904698q
  13. M. Bizzarri, S. Marsili and P. Procacci, “Intraligand Hydrophobic Interactions Rationalize Drug Affinities for Peptidyl-Prolyl Cis-Trans Isomerase Protein,” The Journal of Physical Chemistry B, Vol. 115, No. 19, 2011, pp. 6193-6201. doi:10.1021/jp110585p
  14. M. Bizzarri, E. Tenori, M. R. Martina, S. Marsili, G. Caminati, S. Menichetti and P. Procacci, “New Perspective on How and Why Immunophilin fk506-Related Ligands Work,” The Journal of Physical Chemistry Letters, Vol. 2, No. 22, 2011, pp. 2834-2839. doi:10.1021/jz201037u
  15. C.-Y. Yang, H. Y. Sun, J. Y. Chen, Z. Nikolovska-Coleska and S. M. Wang, “Importance of Ligand Reorganization Free Energy in Protein-Ligand Binding-Affinity Prediction,” Journal of the American Chemical Society, Vol. 131, No. 38, 2009, pp. 13709-13721. doi:10.1021/ja9039373
  16. C. Gao, M.-S. Park and H. A. Stern, “Accounting for Ligand Conformational Restriction in Calculations of Protein-Ligand Binding Affinities,” Biophysical Journal, Vol. 98, No. 5, 2010, pp. 901-910. doi:10.1016/j.bpj.2009.11.018
  17. K. Andries, H. Azijn, T. Thielemans, D. Ludovici, M. Kukla, J. Heeres, P. Janssen, B. De Corte, J. Vingerhoets, R. Pauwels and M.-P. de Bthune, “Tmc125, a Novel Next-Generation Nonnucleoside Reverse Transcriptase Inhibitor Active Against Nonnucleoside Reverse Transcriptase Inhibitor-Resistant Human Immunodeficiency Virus Type 1,” Antimicrobial Agents and Chemotherapy, Vol. 48, No. 12, 2004, pp. 4680-4686. doi:10.1128/AAC.48.12.4680-4686.2004
  18. P. A. J. Janssen, P. J. Lewi, E. Arnold, F. Daeyaert, M. de Jonge, J. Heeres, L. Koymans, M. Vinkers, J. Guillemont, E. Pasquier, M. Kukla, D. Ludovici, K. Andries, M.-P. de Bèthune, R. Pawels, K. Das, A. D. Clark Jr., Y. V. Frenkel, S. H. Hughes, B. Medaer, F. De Knaep, H. Bohets, F. De Clerck, A. Lampo, P. Williams and P. Stoffels, “In Search of a Novel Anti-HIV Drug: Multidisciplinary Coordination in the Discovery of 4-[[4-[[4-[(1e)-2-Cyanoethenyl]-2,6-dimethylphenyl]amino]-2-pyrimidinyl]amino]benzonitrile (r278474, Rilpivirine),” Journal of Medicinal Chemistry, Vol. 48, 2005, pp. 1901-1909. doi:10.1021/jm040840e
  19. K. Das, A. D. Clark, P. J. Lewi, J. Heeres, M. R. De Jonge, L. M. H. Koymans, H. M. Vinkers, F. Daeyaert, D. W. Ludovici, M. J. Kukla, B. De Corte, R. W. Kavash, C. Y. Ho, H. Ye, M. A. Lichtenstein, K. Andries, R. Pauwels, M.-P. De Bthune, P. L. Boyer, P. Clark, S. H. Hughes, P. A. J. Janssen and E. Arnold, “Roles of Conformational and Positional Adaptability in Structure-Based Design of Tmc125-r165335 (Etravirine) and Related Non-Nucleo0 side Reverse Transcriptase Inhibitors That Are Highly Potent and Effective against Wild-Type and Drug-Resistant HIV-1 Variants,” Journal of Medicinal Chemistry, Vol. 47, No. 10, 2004, pp. 2550-2560. doi:10.1021/jm030558s
  20. K. Das, P. J. Lewi, S. H. Hughes and E. Arnold, “Crystallography and the Design of Anti-AIDS Drugs: Conformational Flexibility and Positional Adaptability Are Important in the Design of Non-Nucleoside HIV-1 Reverse Transcriptase Inhibitors,” Progress in Biophysics and Molecular Biology, Vol. 88, No. 2, 2005, pp. 209-231. doi:10.1016/j.pbiomolbio.2004.07.001
  21. K. Das, J. D. Bauman, A. D. Clark, Y. V. Frenkel, P. J. Lewi, A. J. Shatkin, S. H. Hughes and E. Arnold, “HighResolution Structures of HIV-1 Reverse Transcriptase/tmc278 Complexes: Strategic Flexibility Explains Potency against Resistance Mutations,” Proceedings of the National Academy of Sciences, Vol. 105, No. 5, 2008, pp. 1466-1471. doi:10.1073/pnas.0711209105
  22. Y. V. Frenkel, A. D. Clark, K. Das, Y.-H. Wang, P. J. Lewi, P. A. J. Janssen and E. Arnold, “Concentration and pH Dependent Aggregation of Hydrophobic Drug Molecules and Relevance to Oral Bioavailability,” Journal of Medicinal Chemistry, Vol. 48, No. 6, 2005, pp. 1974- 1983. doi:10.1021/jm049439i
  23. Y. V. Frenkel, E. Gallicchio, K. Das, R. M. Levy and E. Arnold, “Molecular Dynamics Study of Non-Nucleoside Reverse Transcriptase Inhibitor 4-[[4-[[4-[(e)-2-Cyanoethenyl]-2,6-dimethylphenyl]amino]-2-pyrimidinyl]amino]benzonitrile (Tmc278/Rilpivirine) Aggregates: Correlation between Amphiphilic Properties of the Drug and Oral Bioavailability,” Journal of Medicinal Chemistry, Vol. 52, No. 19, 2009, pp. 5896-5905. doi:10.1021/jm900282z
  24. H. Okumura, E. Gallicchio and R. M. Levy, “Conformational Populations of Ligand-Sized Molecules by Replica Exchange Molecular Dynamics and Temperature Reweighting,” Journal of Computational Chemistry, Vol. 31, 2010, pp. 1357-1367.
  25. Y. Frenkel, “The Roles of Structural Variability and Amphiphilicity of TMC278/Rilpivirine in Mechanisms of HIV Drug Resistance Avoidance and Enhanced Oral Bioavailability,” PhD thesis, Rutgers University, Piscataway, 2009.
  26. E. Gallicchio and R. M. Levy, “Advances in All Atom Sampling Methods for Modeling Protein-Ligand Binding Affinities,” Current Opinion in Structural Biology, Vol. 21, No. 2, 2011, pp. 161-166. doi:10.1016/j.sbi.2011.01.010
  27. D. L. Mobley, “Let’s Get Honest about Sampling,” Journal of Computer-Aided Molecular Design, Vol. 26, No. 1, 2012, pp. 93-95. doi:10.1007/s10822-011-9497-y
  28. M. Lapelosa, E. Gallicchio and R. M. Levy, “Conformational Transitions and Convergence of Absolute Binding Free Energy Calculations,” Journal of Chemical Theory and Computation, Vol. 8, 2012, pp. 47-60. doi:10.1021/ct200684b
  29. E. Gallicchio and R. M. Levy, “Prediction of Sample Host-Guest Affinities with the Binding Energy Distribution Analysis Method (BEDAM),” Journal of Computer-Aided Molecular Design, 2012, in Press. doi:10.1007/s10822-012-9552-3
  30. E. Gallicchio, M. Lapelosa and R. M. Levy, “Binding Energy Distribution Analysis Method (BEDAM) for Estimation of Protein-Ligand Binding Affinities,” Journal of Chemical Theory and Computation, Vol. 6, No. 9, 2010, pp. 2961-2977. doi:10.1021/ct1002913
  31. D. L. Mobley, J. D. Chodera and K. A. Dill, “On the Use of Orientational Restraints and Symmetry Corrections in Alchemical Free Energy Calculations,” Journal of Chemical Physics, Vol. 125, No. 8, 2006, Article ID: 084902. doi:10.1063/1.2221683
  32. D. L. Mobley, J. D. Chodera and K. A. Dill, “The Confine-and-Release Method: Obtaining Correct Binding Free Energies in the Presence of Protein Conformational Change,” Journal of Chemical Theory and Computation, Vol. 3, No. 4, 2007, pp. 1231-1235. doi:10.1021/ct700032n
  33. E. B. Lansdon, K. M. Brendza, M. Hung, R. Wang, S. Mukund, D. Jin, G. Birkus, N. Kutty and X. H. Liu, “Crystal Structures of HIV-1 Reverse Transcriptase with Etravirine (TMC125) and Rilpivirine (TMC278): Implications for Drug Design,” Journal of Medicinal Chemistry, Vol. 53, No. 10, 2010, pp. 4295-4299. doi:10.1021/jm1002233
  34. M. K. Gilson, J. A. Given, B. L. Bush and J. A. McCammon, “The Statistical-Thermodynamic Basis for Computation of Binding Affinities: A Critical Review,” Biophysical Journal, Vol. 72, 1997, pp. 1047-1069. doi:10.1016/S0006-3495(97)78756-3
  35. F. P. Buelens and H. Grubmüller, “Linear-Scaling Soft-Core Scheme for Alchemical Free Energy Calculations,” Journal of Computational Chemistry, Vol. 33, No. 1, 2012, pp. 25-33. doi:10.1002/jcc.21938
  36. M. R. Shirts and J. D. Chodera, “Statistically Optimal Analysis of Samples from Multiple Equilibrium States,” Journal of Chemical Physics, Vol. 129, No. 12, 2008, Article ID: 124105. doi:10.1063/1.2978177
  37. E. Gallicchio and R. M. Levy, “Advances in Protein Chemistry and Structural Biology,” Recent Theoretical and Computational Advances for Modeling ProteinLigand Binding Affinities, Academic Press, London, Vol. 85, 2011, pp. 27-80.
  38. E. Gallicchio and R. M. Levy, “AGBNP: An Analytic Implicit Solvent Model Suitable for Molecular Dynamics Simulations and High-Resolution Modeling,” Journal of Computational Chemistry, Vol. 25, 2004, pp. 479-499. doi:10.1002/jcc.10400
  39. E. Gallicchio, K. Paris and R. M. Levy, “The AGBNP2 Implicit Solvation Model,” Journal of Chemical Theory and Computation, Vol. 5, No. 9, 2009, pp. 2544-2564. doi:10.1021/ct900234u
  40. W. L. Jorgensen, D. S. Maxwell and J. Tirado-Rives, “Developement and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids,” Journal of the American Chemical Society, Vol. 118, 1996, pp. 11225-11236. doi:10.1021/ja9621760
  41. G. A. Kaminski, R. A. Friesner, J. Tirado-Rives and W. L. Jorgensen, “Evaluation and Reparameterization of the OPLS-AA Force Field for Proteins via Comparison with Accurate Quantum Chemical Calculations on Peptides.” The Journal of Physical Chemistry B, Vol. 105, 2001, pp. 6474-6487. doi:10.1021/jp003919d
  42. J. L. Banks, J. S. Beard, Y. Cao, A. E. Cho, W. Damm, R. Farid, A. K. Felts, T. A. Halgren, D. T. Mainz, J. R. Maple, R. Murphy, D. M. Philipp, M. P. Repasky, L. Y. Zhang, B. J. Berne, R. A. Friesner, E. Gallicchio and R. M. Levy, “Integrated Modeling Program, Applied Chemical Theory (IMPACT),” Journal of Computational Chemistry, Vol. 26, 2005, pp. 1752-1780. doi:10.1002/jcc.20292
  43. R. Pauwels, J. Balzarini, M. Baba, R. Snoeck, D. Schols, P. Herdewijn, J. Desmyter and E. De Clercq, “Rapid and Automated Tetrazolium-Based Colorimetric Assay for the Detection of Anti-HIV Compounds,” Journal of Virological Methods, Vol. 20, No. 4, 1988, pp. 309-321. doi:10.1016/0166-0934(88)90134-6
  44. M. Mihailescu and M. K. Gilson, “On the Theory of Noncovalent Binding,” Biophysical Journal, Vol. 87, No. 1, 2004, pp. 23-36. doi:10.1529/biophysj.103.031682
  45. J. Kim and J. E. Straub, “Generalized Simulated Tempering for Exploring Strong Phase Transitions,” Journal of Chemical Physics, Vol. 133, No. 15, 2010, Article ID: 154101. doi:10.1063/1.3503503
  46. D. H. Min, M. Chen, L. Q. Zheng, Y. H. Jin, M. A. Schwartz, Q.-X. A. Sang and W. Yang, “Enhancing qm/mm Molecular Dynamics Sampling in Explicit Environments via an Orthogonal-Space-Random-Walk-Based Strategy,” The Journal of Physical Chemistry B, Vol. 115, No. 14, 2011, pp. 3924-3935. doi:10.1021/jp109454q
  47. C. A. Chang, W. Chen and M. K. Gilson, “Ligand Configurational Entropy and Protein Binding,” Proceedings of the National Academy of Sciences, Vol. 104, No. 5, 2007, pp. 1534-1539. doi:10.1073/pnas.0610494104

Journal Menu >>