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

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.


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][4][5][6][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][13][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 d tients 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][20][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 com ake 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][23][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 confo co

Methods r Systems
the two DAPY NNRTI in-xes of HIV-RT in co

Abso
f the complexes with th which follows, without approximations, fr rmational 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 nformational 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.

Molecula
The chemical structures of hibitors 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 comple mplex 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.

Computational Protocol
lute binding free energies o 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 om a wellestablished statistical mechanics theory of molecular association [34].In Equation ( 1) M, or equivalently 1668 Å -3 ), ite V s is the volf the binding site (see below), a b G  is the excess free energy of binding defined as the d ence 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 molecu nd ume o iffer the method is the concept of a binding ener (2) defined for each conformation the effec of the complex as the difference between tive potential energies of, respectively, the lig und an ed conformations of the complex without internal conformational rearrangements.BEDAM is based on a λ-dependent hybrid potential energy function of the form is the potential energy of the uncoupled state of the complex.It is easy to see from Equations ( 2)-( 4) that is by definition the difference in free energy ben the where   is the difference d and u  is the difference in binding energies between the icas 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 in λ's being exchange repl  at which they are most likely to occur and to be then pr pagated to other states [26].
To improve c o onvergence of the free energy near 0 =  , in this work we employ a modified "soft-core" where is some large positive value (set in this work kcal/mol).This modified binding energy n, w ratio estimator (MBAR) [36] to estimate the excess fr max u as 10 6 functio hich 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 ee energy of binding b G  , using as input binding energy samples obtained from the HREM simulations.Binding free energies a omputed directly from the MBAR dimensionless free energies  f ˆ using the relationship The MBAR dimensionless free energ are defined as the negative of the loga e ies, , rithm of th  -dependent biased partition functions  Z the dimensionless free energies are estimated by the lf-consistent solution of the set of equat s [36] .In this case se ion where is the h binding energy sample from replica jn u , n t j K is th ber of replicas and is the he B p e num ding j N total nu er of bin energy samples from replica j .For t M AR analysis we employed the code rovided 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] mb

B g Free Energy Calculations with
The concept of the reorganization free energy is usef formalize situations in which, by imposing suitable ational restraints, auxiliary sta which do not match exactly with the end point states of binding energy function of the form 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 restrain G  of restraining the ensembles of conformations of the receptor and the ligand in solution to chosen macrostates * R and * L .For instance, in the application belo * L 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 th population * w the e R L P  , defined as the probability of finding a conformation belonging to the specified macrostate, in absence of restraints and ligand-receptor interactions Following this step, we consider the binding free ener G  gy, , between the that is th nding free ene an , is the free energy d rence between these two states.In the present application we assume release = 0 G  because the remacrostate encompasses all of the conformational ensemble of the complex; that is the * ) (RL and RL species are virtually e Based on the above we have in otion that a binding free energy calculation can be conducted with or conformational restraints, and, if the binding free Equation (10) , obtained with restraints is properly combined with the free energy, restrain G   , of imposing the restraints, either approach should return the same binding Hence the use of restraints is purely a matter of efficiency; restraints should nsidered whenever free energy.be co ed the OPLS-AA/AGBNP2 effective potential in which the effect of the solvent is represented im-2 implicit solvent model t 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.

Computational Details
We employ plicitly by means of the AGBNP [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  , 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 si e volum was defined as any coformation whic 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, site V , is calculated as approximately 33.3 Å 3 and site  -hopping calculations were performed for 8 ns of molecular dynamics per replica (1 ns of mol la namics for each of four nding 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 n from replica j is ass 92 ecu r dy bi igned a  -dependent thermodynamic weight given by where the are the solution of the MBAR E (7) and th er symbols are as for Equation (7).T modynamic averages are formulated as weighted aver-  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 L 1 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 restrain G  defined by Equation ( 9) and reported in Table 1.The restraining free energy for TMC125 ( 0.30 kcal/mol) is much smaller than the one for TMC278 ( 3.21 kcal/mol) reflecting the relative abundance of U conformations for   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 ormational distributions in the bound state.
of restraints this preferentially adopts other conformations.
Despite the large differences in conformational propensities when free in solution, the two compounds have milar conf si 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 3  and 4  degrees of freedom in order to adapt to recept his is confirmed by the or.T the significantly smaller reorganization free energy ( 31.01 kcal/mol;  angles, forcing it to remain withing the U state, compared to unrestrained conditions ( 34.55 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 t distributions of 3 he  and 4  angles.

Binding Free Energy Estimates
The computed binding free energies for the complexes of wild-type HIV-RT with TMC125 and TMC278 inhibitors with and without torsional restraints are shown in Table 1.Based on Equations ( 1) and ( 9  e free e e en to d that di nts are, in addition, often used as computational devices to speed up the ee en ducing the size of conformational space of the complex that needs to be sampled to achieve convergence [6,28,31]  ergy calculations without restraints.The binding free energy estimates (6th column in Ta- ble 1) for TMC278 are in good agreement with the experimental measurement (2nd column).The computed binding free energy values for TMC125 are instead 1 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 note rect 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 50 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.

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 dis-  [30] combined with multi-state ation of conformational free 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 ( 14.96   kcal/mol) is significantly larger in magnitude than the binding free energy of unrestrained TMC278 ( 11.37   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, t two approaches yield equivalent binding free energy es-

he
timates 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 31 kcal/mol with restraints compared to 34.5 kcal/mol without them.The larger reorganization by 3.5 kcal/mol without restraints is clearly due to the reorganization free of the ligand, estimated at approximately 3.2 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).

Effects of Restraints on the Rate of Convergence
As discussed above increasing computational eff is the main motivation for using conformational restraints in binding free energy calculations.One measure of putational 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) iciency com .The estimated unce rg e (the results presented in Table 1 are from the sec nd half of the trajectory); lo nd the M ativ better conformational rtainty for the binding free energy of TMC125 with restraints is significantly larger than for the unrestrained calculation ( 0.5 vs 0.1 kcal/mol, respectively).The statistical uncertainties for the two calculations for TMC278 binding are very similar (both are approximately 0.65 kcal/mol) with the uncertainty ( 0.2 kcal/mol, see Table 1) of the restraining free energy, restrain G  , 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 ene estimates vary as the simulation progresses (Figure 7).In this figure we plot binding free energy estimates fro ata 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 equilibratio 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 w bi ing 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 complexes with TMC125 and T C278 with and without torsional restraints.In general, a larger number of transitions from high binding energies to low binding energies and vice versa is indic kcal/mol) and high binding energies ( kcal/mol).More rarely and similarly to the unr d simulation, long-lived transition to high binding occur (see for example bottom-center panels i 9).The simulations for TMC125 binding 8) display analogous behavior except that the k the unrestrained and restrained simulations are m re similar to each other and undergo more long-live sitions than the TMC278 simulations.
Structural analysis shows highto-low rapid transitions in binding ene observed in the TMC125 simulations and th ed TMC278 simulation (Figures 8 and 9 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 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 n orientation suitable for binding.This is due i a the presence of the extra orientational con 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, L 1 , and L 2 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 release ΔG 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.

Can Restraints Slow Down Convergence?
One puzzling result of this work is that the statistical uncertainty for the binding free energy of TMC125 esti-mated 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 L conformations w 1 hen dissociated from the receptor and exclusively U conformations when it is fully associated.However, the shift from U and L 1 conformations to exclusively U conformations does not appear to occur as a direct population transfer from L 1 to U. At intermediate values of  the L 2 and E conformational states become also significantly populated (Figure 10).This occurs mainly by population transfer from the L 1 macrostate to the L 2 and E macrostates (Figure 11); with increasing  , the populations of the L 2 and E states increase rapidly at the expense of the L 1 population.At larger values of  's, population transfer from the L 2 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 L 1 conformation, reach the bound U state by going through E and/or L 2 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.

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   This behavior is the tell- (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 o thereby inding ed field Our orig was that, due to its higher flexi-f 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) constiutes an obstacle to binding/unbinding transitions, t hampering convergence of b free energy calculations.Advanced conformational sampling tools are being develop in this and other research s to address conformational sampling challenges of this kind [45,46].

Insights on Thermodynamic Driving Forces
for Binding inal expectation  2 the binding reorganization free energy of TMC125 is approximately 36 kcal/mol; 2 kcal/mol greater than the reorganization of TMC278.This is despite the fact that reorganization of the 3

 and 4
 torsional angles alone within the defined ranges (see above) opposes binding of TMC278 by approximately 3 kcal/mol (see restrain G  in Table 1, 5th column).Clearly, collective reorganiz tion of other degrees of freedom of TMC125 must be contributing more to the overall reorganization free energy than the 3 a-

 and 4
 angles which account for only approximately 0.3 kcal/mol of reorg G   .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 1 =  reveals that, when bound to HIV-RT, TMC278 is significantly more mobile than TMC125 in both respects.amino-dimethylbenzocyanovinyl moiety of T 278.This data clearly indicates that the sidechain of TMC278 undergoes a wider range of libration motion than TMC125, despite the fact that the cyanovin l gro 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 2  torsional angles of the amino-benzonitrile gr , wh 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 midine 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 o 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 Ta , 4th column), compared to a deficit of roughly ol 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 rel e to TMC125 if not for its alty for is key to the greater potency of TMC278 against a wide riginate 3.5 ble 2 kcal/m 2 ativ 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) nd the lower reorganization free energy pen a 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 range of drug-resistant mutants of HIV-RT.
Structure-based drug optimization trategie are often aimed at increasing the number and strength of ligandreceptor interactions.This example shows that s s both gr tha for or against binding have been obtained.C he formation of ynamics of protein-ligand binding.supported by grants by the National (GM30580) and the National Science eater affinity and improved resistance profiles can be achieved by focusing instead on increasing the redundancy and flexibility of ligand-receptor interactions so t they are maintained over a wider range of conformations of the complex.

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 ontrary 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 penal-ties and the effects of resistance mutations, suitably modified as in TMC125 so as to disfavor t non-binding competent extended conformations (L 1 , L 2 and E) in solution.

Figure 1 .Figure 2 .
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 effective potential en funct of the coupled and decoupled states of the complex, respectively.The excess binding free energy b G  are limited to the chosen macrostates.In genera g state of the c plex, d noted by ( RL ) * in Figure3, does not match the full complex state RL because in the former the receptor and the ligand are limited to their respective macrostates.

Figure 3 .
Figure 3. Thermodynamic cycle illustrating the restrainand-release decomposition of the excess free energy of binding [Equation (10)].Although not indicated, he ligand

2 . 5 ,
bservations using the weights defined by Equation(11).For example the average binding energy at angle for instance) at  is expressed as jn and the value of the structural quantity of sam-) ( jn    nt falls within the is a function defined as 1 if the argum bin centered at e  with width   and zero AM d iga e co otherwise.BED binding free energy calculations were conducte for TMC125 and TMC278 with and without torsional angle restraining potentials designed to restrict the l nd within the bound ligand macrostate, d fined as any nformation with 3  and 4  torsional angles within a  120 interval around  0 =  (see Figures 1 and 3).Restraining potentials were implemented as flat-bottom harmonic restraints on the 3  and 4  torsional angles centered on 0 =  with a to rance of  60 on either side.opulation  L R P  of the ound macrostate at 0 = le b The p  , required to compute restrain G from Equation (9), was evaluated by numerical integ tion of he joint probability densityyed in terms of the corresponding potential of mean force (PMF) the smallest value of the PMF as the zero of energies.in the largest conformational variations3.1.Conformational Analysis of Free and Bound TMC125 anSince they re (see Figure4), we will focus first on the 3  and 4  torsional degrees of freedom of the ligands as defined in re The computed conform Figu ational distributions for TMC125 and TMC278 free in solution with respect to these torsional angles (see Methods) are shown in Fig- ures 5 and 6, respectively.As it can be seen from Figure TMC125 assumes predominantly U conformations in solution with a minor L 1 component.Remarkably, TMC278 displays the opposite behavior by occupying mostly extended (E) and L 2 conformations (Figure

Figure 5 .
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 Figure1).The regions corresponding to the main conformational states (see Figure3) are indicated.

Figure 6 .
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 Figure1).The regions corresponding to the main conformational states (see Figure3) are indicated.
complexes of HIV-RT ith the TMC125 a 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.

Figure 7 .Figure 8 .
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.

Figure 9 .
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).
(b))shows better behavior; all the replicas display ra d transitions between low bi ng energies (in the bound state at low binding energies.

Figure 10 .
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 Figure1).At this intermediate value of λ TMC125 visits all of the main rotamers (see Figure3).

Figure 11 .
Figure 11.Conformational populations of TMC125 as a function of λ from the unrestrained binding simulation.At λ = 0 only the U and L 1 conformations are significantly populated.As λ increases the L 1 population decreases rapidly concomitantly with a rapid increase of the populations of the L 2 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 ) (u p  0.2 <  binding energies are predominantly above 20  kcal/mol, whereas for 0.35 >  much lower val r Values of bindin ues of g ener binding ene gies near gies are sampled.20  kcal/mol are rarely visited at

Figure 12 .
Figure 12.Probability densities of the binding energy, u, for

Figure 13 .
Figure 13.Conformational populations of the U, L 1 , L 2 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

Figure 14 . 4 
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 m (b) shows the distribution of distances oup ich is d  of the position of the N 1 atom of the central pyri-

Table 2 ,
5th column) obtained for TMC278 binding with torsional restraints on the 3

Decomposition of computed binding free energies into average binding energy
a In kcal/mol.Binding free energy with restraints from Tab nization free energy with restraints (sam ganization re-b le 2. c Reorga e as total reor free energy for un strained calculations).d Total reorganization free energy including restrain ΔG from Table 1.