rna performs a vast 
 array of functions in biological systems , including 
 genetic encoding , regulation , and catalysis , and yet very 
 few small - molecule drugs that target rna exist . this may be the result of many factors , including the relatively 
 recent discovery of rna s many biological roles and the difficulty 
 in preventing rna degradation during experiments , particularly by 
 ribonucleases . likewise , computational investigations 
 of rna  ligand binding are comparatively rare ( a pubmed search 
 of  protein binding simulations  as of january 2014 
 yielded 7633 results , and a search of  rna binding simulations  
 yielded 488 results ) . in order to address this paucity , 
 the current study reports the results of molecular dynamics ( md ) simulations 
 on a specific rna ligand system and aims to provide a more 
 reliable foundation for future studies involving highly charged rna the target of this research 
 is the domain iia rna sequence from 
 the hepatitis c virus internal ribosome entry site ( hcv ires ) . experimental structures exist for the unbound 
 ( or free ) structure and also of the rna in complex 
 with 2-aminobenzimidazole inhibitors . inhibitor 
 complexes are attractive structures to study because they involve 
 a relatively short rna sequence bound to druglike molecules . this 
 contrasts with typical structures that are often larger and more complex , 
 such as rna or riboprotein molecules in complex with aminoglycosides . moreover , a distinct structural difference between the free and 
 bound hcv ires is observed , and this is most notably characterized 
 by the loss of a critical bend in the rna upon ligand binding that 
 explains the inhibition mechanism . biologically , 
 the structure is of interest because of both the high degree of sequence 
 conservation in ires elements and its importance in hcv genome translation 
 and viral replication . rather than using 
 the 5 cap - dependent mechanism to initiate translation at the 
 ribosome , as is typical in eukaryotes , thus , the development of 
 inhibitors of the ires machinery could be useful in treating hepatitis 
 c virus infections . the 2-aminobenzimidazole inhibitors used 
 in the experimental structures 
 were developed by isis pharmaceuticals , inc . using a structure these rna 
 binding inhibitors were confirmed to reduce hcv rna levels in a viral 
 rna replication assay . as part of the 
 exploration of sars , a number of different derivatives were synthesized 
 and binding constants estimated ( those studied in this work are described 
 in figure 1 and table 1 ) . this provides a series of related inhibitors studied by the same 
 laboratory with equivalent and comparable experiments that can be 
 investigated by simulations to assess biomolecular simulation protocols . there are some drawbacks to this experimental data set , including 
 the following : ( 1 ) the protonation state of the inhibitor upon binding 
 is unknown ; ( 2 ) several inhibitors were synthesized as mixtures of 
 enantiomers or diastereomers , and the experimental binding data published 
 do not distinguish the effects from individual stereoisomers ; and 
 ( 3 ) the errors in the binding measurements were not reported . for example , 
 the protonation states can be estimated with reasonable accuracy using 
 pka estimation software ( see methods ) , and such calculations suggest that the 
 inhibitors are all fully protonated in solution at physiological ph 
 as depicted in figure 1 . with regard to stereochemistry , 
 it is very easy to perform separate calculations on each of the enantiomers 
 and diastereomers and , as an approximation , to report the mean value 
 for comparison with experimental data for the mixture . the lack of 
 error analysis in the experimental results , however , does suggest 
 the use of caution when making certain conclusions on the basis of 
 a comparison of the experimental and computational results . structures 
 of previously identified inhibitors , which 
 are studied further in this work , that bind to the 
 hcv ires subdomain iia . the locations of the 
 stereocenters in j3 , j4 , and j5 are indicated by a  symbol 
 on the chiral carbon . we refer to the identity of the stereoisomer 
 by adding r or s to the name designator ( i.e. , j3r , j3s , j4r , j4s , 
 j5r , and j5s ) . the two carbon stereocenters of the j6 diastereomers 
 are indicated with the  and  symbols , and denoted 
 by appending r or s to j6 in that respective order ( j6 ) : 
 j6rr , j6rs , j6sr , and j6ss . note : the j4 inhibitor was used in the 
 nmr study , whereas the j5 inhibitor was 
 used in the crystallography study . the free energy of binding was calculated 
 according to the relation g = rt ln  kd at 298.15 k. the exact dissociation constant 
 is not known for j1 , and thus , 5.45 kcal / mol represents the 
 lower bound for the binding energy . the available experimental structures of the rna  inhibitor 
 complexes derived from nmr data and x - ray 
 crystallography exhibit distinctly different 
 binding modes . these differences can not easily be attributed to the 
 identity of the two different inhibitors used in the separate investigations , 
 as they differ by only a single ch2 group 
 ( j4 for nmr and j5 for x - ray ; see figure 1 ) . the differences also can not be attributed to the slight change in 
 rna sequence in the nmr structure versus the crystal structure , as 
 the sequence changes are found only at the ends of the molecule , as 
 shown and discussed in figure 2a . rather , the 
 differences must lie in either the experimental conditions or the 
 structure refinement procedures used in the two experimental methods . the nmr structure can be described as an open conformation with stacking 
 contacts formed below the inhibitor . contacts to the rna binding site 
 are formed by both dimethylamino arms of the inhibitor : one to the 
 phosphate group of g5 and one to the junction of the c11:g33 base 
 pair ( figures 2b and 3a ) . hydrogen - bonding contacts between the benzimidazole ring donor 
 sites and the rna are not observed . in contrast , the crystal structure 
 is characterized by a compact binding site in which stacking interactions 
 are formed both above and below the inhibitor . both dimethylamino 
 groups of the inhibitor interact with the rna phosphate backbone , 
 and two critical hydrogen bonds are formed between the benzimidazole 
 ring and residue g33 ( figures 2c , 3b , and 4 ) . suggest that these hydrogen bonds explain the 
 observation that the c11g : g33c base pair mutant does not bind the 
 benzimidazole inhibitors . despite similar sequences , the reported conformations 
 of the inhibitor - bound 
 hcv ires domain iia determined by nmr analysis and x - ray crystallography 
 differ . ( a ) secondary structure diagrams of the domain iia constructs 
 used in the nmr and crystallography studies . the hairpin sequence from the nmr study 
 was used for all of the simulations in this study . the residues colored 
 in red show the portions of the rna that are identical in the two 
 published structures . ( b , c ) representative models depicting the global 
 structures of the nmr ensemble and the crystal structure , respectively . 
 in each structure , the rna backbone is emphasized with heavier width , 
 and the inhibitor is highlighted in green . the structural orientations 
 were chosen to emphasize the global differences in the binding conformations . stereoviews ( wall - eyed ) of the inhibitor binding 
 conformations 
 observed in ( a ) the experimental nmr ensemble and ( b ) the crystal 
 structure . schematic of the crystal structure binding site in the plane of 
 the j5 inhibitor . residue labels are numbered according to the nmr 
 hairpin sequence ( figure 2a ) . critical contacts 
 are indicated by colored dotted lines , and the colors correspond to 
 the distances depicted in si figure s2 . in addition to the differences 
 in the conformations of the nmr 
 and crystal structures , there are potentially conflicting reports 
 regarding the cation requirements necessary for the formation of the 
 inhibitor - bound complex . magnesium is observed at core positions in 
 both the unbound and bound crystal structures . it has also been found that removal of magnesium from the frster 
 resonance energy transfer ( fret ) binding assay of the crystal structure 
 yields an approximately 30-fold decrease in binding affinity . although these findings are consistent with 
 the well - known relationships between rna structure and cation binding , they should not be interpreted to suggest that coordinated magnesium 
 exclusively performs the stabilization role . a variety of rna tertiary structures are known to form in moderate 
 levels of monovalent salt , and magnesium 
 is known to compete with monovalent cations in stabilizing rna . in the case of the hcv ires domain iia rna , 
 although the addition of magnesium stabilizes the unbound solution 
 structure , no changes were observed in 
 the nmr spectra of the bound complex upon addition of magnesium to 
 a solution with a relatively high monovalent salt concentration . also , the first reported dissociation constants 
 for the benzimidazole inhibitors were determined in the absence of 
 magnesium with modest levels of ammonium acetate . these values ( table 1 ) are comparable 
 to those determined by fret assays using the crystallography construct 
 with magnesium . additionally , a fluorescence binding assay of the 
 nmr sequence conducted in 0.15 m kcl and 0.15 m nacl using the j4 
 inhibitor ( figure 1 ) yielded a dissociation 
 constant of 2.4 m , which is in 
 the equivalent range as the value determined by a fret assay for the 
 crystal sequence bound to the j5 inhibitor in 2 mm mg ( ec50 = 3.4 m ) . the goal of this work is to validate the computational methods 
 we plan to employ for hcv ires inhibitor discovery and/or optimization . 
 as an initial step in this process , we report results that attempt 
 to answer the following questions : ( 1 ) which of the two published 
 ires ( 2 ) can md simulations further inform whether the magnesium ions identified 
 in the crystal structure are required for inhibitor binding ? ( 3 ) is 
 the mm - gbsa / mm - pbsa method sufficiently accurate to predict the relative 
 order of binding affinity of several 2-aminobenzimidazole inhibitors ? ( 4 ) is there a significant difference in docking results when using 
 the partial charges assigned by the rigorous red method versus the more approximate am1-bcc method ? ( 5 ) do the md simulations of various 2-aminobenzimidazole inhibitors 
 suggest alternative inhibitor scaffolds that can be validated using 
 further md simulations ? all of the inhibitors in this 
 study were protonated at the dimethylamino and benzimidazole positions 
 as indicated in figure 1 . the fully protonated 
 state at ph 7.0 is consistent with pka estimates by two different pka prediction 
 programs , sparc and marvin sketch ( www.chemaxon.com/products/marvin ) . charge derivation was performed in a very careful manner because 
 of the highly charged nature of the inhibitors : ( 1 ) a hand - built inhibitor 
 model was geometry - optimized at the quantum - mechanical ( qm ) hf/6 - 31 g * 
 level consistent with the amber ff10 and ff12sb force fields , after 
 which the initial atomic charges were determined by restrained electrostatic 
 potential ( resp ) charge fitting ; ( 2 ) 
 50 ns of implicit - solvent generalized born ( gb ) md using the hawkins , 
 cramer , and truhlar model was performed 
 at 400 k to sample relevant inhibitor conformations , and the resulting 
 trajectory was clustered into 20 clusters 
 by the  averagelinkage  algorithm using amber s 
 ptraj program ; ( 3 ) the representative structures from clusters whose 
 occupancies were greater than 2% were then geometry - optimized at the 
 qm hf/6 - 31 g * level ; and ( 4 ) optimized structures whose energies were 
 within 0.5 kcal / mol of the minimum - energy structure were used 
 in a multiconformation , multiorientation resp fit using the red program to generate the final charges used in this study . bonds , 
 angles , torsions , improper torsions , and lennard - jones parameters 
 were assigned from the general amber force field ( gaff ) using the 
 antechamber and parmchk programs . some torsion 
 and improper torsion parameters were modified because the default 
 parameters did not maintain planarity at c2 of the 2-aminobenzimidazole 
 ring . mol2 files with the modified torsion parameters are provided in the supporting information ( si ) . all of the qm calculations 
 were performed using gaussian 09 ( www.gaussian.com ) , and all of the md simulations were performed using amber12 . in addition to the more detailed approach 
 to generate high - quality atom charges discussed above , we were interested 
 in the performance of more approximate charge parametrization methods . 
 to test this , we performed molecular docking studies ( discussed below ) 
 using both the resp charges from the above procedure and am1-bcc charges produced by amber s antechamber program . in the latter 
 case , experimentally 
 determined atomic - resolution structures exist for hcv ires rna complexed 
 with the j4r and j4s inhibitors ( nmr structures ) and also for the 
 j5r inhibitor ( crystal structure ) . to facilitate comparisons between 
 the md simulations , the crystal structure duplex was converted into 
 a hairpin of identical sequence to the nmr structure . to accomplish 
 this , the 3 dangling bases were removed , the c : g base pair 
 at the base of the lower stem was converted to a g : c base pair , and 
 a uucg tetraloop was added to the opposite stem ( figure 2a ) . with the exception of j4r and j4s bound in the nmr conformation 
 and j5r in the crystal conformation , binding poses for the inhibitors 
 in our test set were not available . in preliminary md simulations , 
 we noticed that the flexible portions of the inhibitor scanned the 
 nearby rna contacts on a short ( ns ) time scale , suggesting a dynamic 
 binding mode . while the core aromatic ring atoms of each inhibitor 
 are presumed to bind identically , the best pose for the flexible regions thus , we chose the following procedure to eliminate 
 bias in building the initial poses . in order to generate a diverse 
 set of binding conformations for the inhibitors in figure 1 , the 20 representative conformations of each inhibitor 
 identified by clustering during the charge derivation procedure were 
 root - mean - square ( rms)-fit to the benzimidazole core atoms of the 
 existing experimental structure , either nmr or crystal . for the nmr 
 structure , the first model from the protein data bank 2ku0 ensemble was used as the reference structure . this procedure 
 resulted in the generation of 20  nmr - like  rna  inhibitor 
 conformations and 20  crystal - like  rna  inhibitor 
 conformations for each of the 12 stereochemically distinct inhibitors . the selection of these conformations as initial structures for the 
 various simulation sets in this work is described in the following 
 section . the simulation that included magnesium used the exact coordinates 
 from the crystal structure for all of the atoms of the j5r inhibitor , 
 magnesium , and rna , except for the necessary sequence modifications 
 to the rna helix termini required to match the nmr sequence . in the 
 case of the novel ligands , a single inhibitor conformation was chosen 
 and rms - fit to the benzimidazole core atoms in the experimental crystal 
 structure . as described in 
 table 2 , inhibitor 
 studies , two strategies were employed : single long simulations and 
 multiple short simulations . for the single long simulation sets ( nmr1 
 and cry1 ) , a single initial structure was selected from the 20 initial 
 conformations for each of the 12 inhibitors on the basis of the minimum 
 gb energy of the complex . for the multiple simulation set ( cry2 ) , 
 all 20 initial conformations were used . in a few cases , bad initial 
 conformations with severe atom overlap simulation length represents the 
 simulation time for each simulation ( the product of the total number 
 of simulations and the simulation length yields the aggregate time ) . +  indicates that the stated time is the minimum from 
 among the simulation set . the domain iia rna 
 was parametrized using amber s ff12sb force field , which includes 
 a recent update of the backbone and  torsion parameters . the initial rna inhibitor conformations were first minimized 
 for 2500 cycles using the steepest - descent algorithm in implicit gb 
 solvent , and the resulting geometries were solvated . all of the simulations 
 described in table 2 were performed in tip3p 
 water with net - neutralizing potassium 
 ions and an additional 200 mm kcl as parametrized by joung 
 and cheatham . the number of waters added 
 was chosen to yield a periodic truncated octahedron with an approximately 
 12  minimum water shell between the solute and the box edge . 
 in order to facilitate energetic comparisons of inhibitors , the numbers 
 of solvent atoms for the systems in each simulation set were made 
 to be identical using an in - house perl script coupled to amber s 
 leap program . in the case of the j1 inhibitor , which has a net charge 
 of + 2 rather than + 3 , direct energetic comparisons with the other 
 inhibitors were not performed . following solvation , the monovalent 
 ion positions were randomized with amber s ptraj program using 
 the  randomizeions  command to remove bias from the 
 initial ion placement . in the case of the mg simulation ( table 2 ) , the crystallographic magnesium ions and water 
 molecules were included using the magnesium parameters of allnr , 
 nilsson , and villa in addition to 200 
 mm kcl . all of the solvated 
 simulations used a similar minimization , heating , and equilibration 
 procedure : ( 1 ) the entire system was minimized for 1000 steps using 
 the steepest - descent algorithm followed by 1000 steps of conjugate - gradient 
 minimization while 25 kcal mol  positional restraints were enforced on the rna and inhibitor benzimidazole 
 core atoms ; ( 2 ) the system was heated from 10 to 150 k at constant 
 volume with the langevin thermostat over the course of 100 ps with 
 25 kcal mol  positional 
 restraints on the rna and benzimidazole core atoms ; ( 3 ) further heating 
 from 150 to 298 k and initial equilibration were performed using constant 
 pressure and the langevin thermostat over the course of 100 ps with 
 5 kcal mol  positional 
 restraints on all of the solute atoms ; and ( 4 ) final equilibration 
 at 298 k was performed for 2 ns using constant pressure and the langevin 
 thermostat with 0.5 kcal mol  positional restraints on the rna and benzimidazole core atoms . production 
 simulations were performed at 298 k at constant pressure using the 
 weak - coupling algorithm for the thermostat and barostat . the pressure relaxation times were 1 ps for 
 the initial equilibration step , 5 ps for the final equilibration step , 
 and 10 ps for production . for heating and both equilibration steps , 
 a collision frequency of 2 ps was used for the 
 langevin thermostat . for production , the time constant of heat bath 
 coupling was 10 ps using the weak - coupling algorithm . for all heating , 
 equilibration , and production steps , the time step was 2 fs , the direct 
 space sum used a cutoff of 8.0  , and shake was applied to all 
 bonds involving a hydrogen atom . the 
 default particle - mesh ewald settings 
 ( which correspond to a grid spacing of 1  and a direct 
 space tolerance of 10 ) were used to determine 
 long - range charge interactions . coordinates were recorded every picosecond 
 during production simulations . with the exception of the simulation 
 performed with residual dipolar coupling ( rdc ) restraints , all of 
 the production simulations were unrestrained and performed using either 
 the cpu or gpu version of the pmemd program in amber12 . a single simulation of the crystal conformation 
 with the nmr rdc restraints enforced was performed using amber s 
 sander program ( the use of rdc restraints is not yet implemented in 
 the faster pmemd program ) . a short minimization was performed on the 
 equilibrated , solvated structure to best fit the rdc alignment tensor . the relative weighting of the alignment restraint was chosen to be 
 0.08 kcal / hz , which represents an empirical determination 
 of the largest value that did not produce simulation instability ( e.g. , 
 integration errors ) . mm - gbsa and mm - pbsa 
 are well - known 
 postprocessing techniques for computing binding energies from simulation 
 trajectories . in this work , mm - gbsa and mm - pbsa analyses were performed 
 with the mmpbsa.py program in amber12 using the single - trajectory 
 approach . for mm - gbsa , the hawkins , cramer , and truhlar gb model was used for implicit solvation with a salt 
 concentration of 200 mm approximated using debye  huckel screening . 
 for mm - pbsa , the following options were used : a level - set - based dielectric 
 model , a two - term nonpolar solvation free energy term based on a cavity 
 and dispersion calculation , an ionic strength of 
 200 mm , a solvent probe of 1.4  , a grid spacing of 0.25  , 
 and a value of 6.0 for the ratio between the longest grid dimension 
 and the solute . radii for inhibitors were chosen from a set of optimized 
 radii to best match the atom types present . for mm - gbsa , all of the frames from the trajectory were used , but 
 for mm - pbsa only every 100th frame was used for computational efficiency 
 reasons . solute entropy estimates were calculated separately ( see 
 below ) and were not included in the mm - gbsa / mm - pbsa values . in addition to the mm - gbsa / mm - pbsa framework of energy analysis , 
 we also calculated relative binding enthalpies that included the full 
 explicit solvent effects . the instantaneous enthalpy for any given 
 simulation frame is defined aswhere upot is 
 the potential energy , ek is the kinetic 
 energy , p is the pressure , and v is the volume . when the binding enthalpy was computed , the kinetic 
 energy and pressure  volume terms were assumed to be negligible 
 because of the use of the thermostat and barostat . thus , the relative 
 binding enthalpy was calculated by subtracting the solvated - inhibitor 
 mean potential energy ( obtained using simulations of the free ligands 
 in explicit solvent , denoted as lig ) from the solvated rna  inhibitor 
 mean potential energy ( obtained from the cry1 and cry2 simulations):the inhibitor j1 was excluded from these calculations 
 because its net charge differs from those of the other inhibitors , 
 which complicated direct comparisons because of differences in the 
 numbers of counterions . only the inhibitor entropy was considered , and entropy changes were 
 computed by subtracting the entropy of the free inhibitor from that 
 of the complexed inhibitor . the rna  inhibitor and free - inhibitor 
 conformations were taken from the cry1 and lig simulation sets , respectively . the first method , quasi - harmonic analysis , was computed using amber s ptraj program . the second method , 
 first - order configurational entropy analysis based on bond / angle / torsion 
 probability distribution functions , was computed using the accent 
 program developed by gilson and co - workers . because of the accumulative nature of these values and 
 the difficulty we had in converging them , error bars are not given 
 for the entropy estimates , and the data are not used in combination 
 with the other energetic analyses presented here . finally , unbound - inhibitor 
 solvation enthalpies were estimated 
 by subtracting gas - phase average potential energies from either gb 
 implicit solvation energies or from potential energies of explicitly 
 solvated trajectories . although one could compute the gas - phase energies 
 from solvated trajectories by stripping the solvent , we performed 
 100 ns gas - phase simulations for each of the 12 inhibitors represented 
 in figure 1 in order to ensure the independence 
 of these values . two approaches were 
 used to estimate 
 the error , depending on the simulation set . for the single long simulation 
 sets , briefly , a data set plotting the standard error 
 of the mean ( sem ) versus increasing block size was computed . given 
 sufficient sampling , the plot plateaus at a value that corresponds 
 to the sem . the data can be fit with a trend line to predict this 
 value , but the fit is not always accurate because of specific assumptions 
 about the type of correlation in the data . to be conservative , we 
 used the maximum value observed in the plot . for 
 the multiple short simulation set ( cry2 ) , we considered the average 
 values from the 20 separate simulations to be independent , uncorrelated 
 data points and computed the sem in the traditional way by dividing 
 the sample standard deviation by the square root of the number of 
 data points ( i.e. , the number of simulations ) . the error combining 
 rules were as follows : when the difference of two mean values was 
 computed , the errors were added ; when the average of two or more mean 
 values was computed , the errors were averaged . the regions of highest magnesium ion 
 occupancy were determined using the  grid  command in 
 amber s ptraj program . the simulation 
 trajectory frames were centered , imaged , and rms - fit using the heavy 
 atoms of residues 5 , 11 , 33 , and 34 that form the binding region . occupancy was determined by a three - dimensional histogram approach 
 using a 75   75   75  box with 0.5  
  0.5   0.5  resolution . inhibitor structure from the simulation 
 using the ucsf chimera package . to choose 
 the density surface contour level to be displayed , the contour level 
 was increased until magnesium occupancy in the bulk solvent region 
 was no longer observed , thus suggesting stable binding locations . docking was performed on the crystal receptor 
 structure ( modified to match the full nmr sequence ; figure 2a ) using dock 6.5 . the 
 top and bottom helical portions of the receptor were excluded from 
 consideration as they are known not to contain the binding site . this 
 exclusion did not inappropriately limit docking poses to the known 
 binding cavity since the entire backside region of the receptor was 
 explored for docking . in order to include various ring pucker conformations 
 in some of the inhibitors , for which dock 6.5 is not able to search 
 automatically , all of the inhibitor conformational clusters whose 
 occupancy was greater than 2% ( identified during the charge derivation 
 procedure ) where used as initial seed structures for docking . two 
 schemes were used to assign charges to the inhibitor for use during 
 docking . the first scheme simply used the resp charges that were derived 
 for use in the md simulations . in the second scheme , which resembles 
 a more typical docking procedure , charges were derived for each inhibitor 
 conformation using the semiempirical am1-bcc charge model , which can 
 be accessed through amber s antechamber program . the default 
 grid - based method in dock 6.5 was assigned as the primary scoring 
 function . nmr data were acquired 
 on 500 and 
 600 mhz varian nmr spectrometers using rna - optimized pulse sequences 
 from the standard varian biopack pulse sequence library ( www.agilent.com ) . the f1-filtered / f2-filtered 
 nuclear overhauser spectroscopy ( noesy ) measurements were performed 
 as described by zwahlen et al . the nmr 
 sample conditions and data acquisition were as described by paulsen 
 et al . , and the data were processed using 
 standard vnmr processing software and then visualized in sparky ( www.cgl.ucsf.edu/home/sparky ) . the nmr1 and cry1 
 simulation sets ( table 2 ) were intended to 
 evaluate the simulation stability 
 of the two available experimental conformations in the context of 
 12 related inhibitors over a fairly long time scale ( 200 + ns ) . visual 
 inspection of the simulation trajectories reveals a stark contrast 
 in the stability of the binding region . a quantitative measure of 
 this difference is shown in figure 5 , where 
 the binding region rms deviation ( rmsd ) is plotted versus simulation 
 time using the initial conformation of each production simulation 
 as a reference . for all 12 simulations in the cry1 simulation set 
 ( shown in black ) the rmsd value is low and very steady . the small 
 fluctuations arise from inhibitor ring transitions and conformational 
 searching by the dimethylamino groups . throughout the cry1 simulations , in contrast to the cry1 results , 
 the nmr1 rmsd results ( shown in red ) reveal a high degree of fluctuation 
 and departure from the initial structure . a variety of rna  inhibitor 
 poses are adopted during the 12 nmr1 simulations and do not point 
 to a consensus alternative to the original nmr pose . much of the structural 
 instability is due to conformational transitions in the bulge residues 
 610 . to visualize the structural difference between the nmr1 
 and cry1 simulation sets , the final frames of each simulation were 
 overlaid to produce an ensemble ( figure 6 ) . 
 as a result of unfolding of the bulge residues , many of the nmr1 conformations 
 do not retain the linear rna orientation known to inhibit viral replication 
 but rather adopt an additionally , 
 the inhibitor poses are smeared into a variety of orientations ( figure 6c ) . in contrast , the cry1 simulations maintain the 
 linear orientation and produce a rather tight ensemble of structures 
 and inhibitor poses ( figure 6b , d ) . we conclude 
 from these data that the nmr conformation is highly unstable in the 
 context of the force field energy landscape , whereas the crystal conformation 
 appears to be in an energetically favorable minimum . plot of binding 
 region rmsd vs simulation time , which reveals that 
 the 12 simulations in the cry1 simulation set ( black ) are much more 
 stable than those of the nmr1 set ( red ) . the simulation sets contain 
 12 simulations , one each for the 12 inhibitor - rna complexes ( figure 1 ) bound in either the crystal conformation ( cry1 ) 
 or the nmr conformation ( nmr1 ) . the atoms considered in the binding 
 region are dened to be the heavy atoms in residues 5 , 6 , 32 , 
 33 , and 34 and the inhibitor . the rst frame of each production 
 simulation was used as the rmsd reference structure for that simulation . 
 for clarity , the rmsd values ( a , b ) structural ensembles 
 made from the final frames of the 12 
 simulations in the ( a ) nmr1 and ( b ) cry1 simulation sets as defined 
 in table 2 . the ensembles were generated by rms fitting of residues 
 5 , 6 , 33 , and 34 of each frame . ( c , d ) final inhibitor position in 
 each simulation for the ( c ) nmr1 and ( d ) cry1 simulation sets using 
 the same rms fits as for the full ensembles . one of the nmr1 simulations , the complex with the j5s inhibitor , 
 is notable because it partially converts to the conformation of the 
 crystal structure . all of the critical contacts in figure 4 are observed , including the hydrogen bonds between 
 the inhibitor and g33 and the contacts between the dimethylamino groups 
 and the phosphate backbone ( si figure s2a ) . these distances are similar to those observed in the cry1 simulation 
 set ( si figure s2b ) . the only interaction 
 that is not formed is the base - triple interaction between residue 
 a6 and the hoogsteen edge of a32 , which forms a  roof  
 above the inhibitor . after 232 ns of simulation , the conformation 
 of residues 79 continues to restrain the flexibility of a6 
 in such a way as to prevent the full formation of the base triple , 
 although a6 and a32 are near enough to form direct contacts . other 
 simulations in the nmr1 set also partially adopt the crystal - like 
 binding mode , usually by forming the hydrogen bonds between the inhibitor 
 and g33 , but none are as stable as the aforementioned simulation with 
 j5s . further research , likely involving enhanced sampling , is necessary 
 to determine whether a complete transition is accessible on a reasonable 
 time scale . in addition to using long simulations , we also investigated 
 whether 
 using many shorter simulations with diverse initial inhibitor conformations 
 could produce comparable data . because of the instability of the nmr 
 conformation in the simulations , we discuss this approach only for 
 the crystal conformation , and we term this the cry2 simulation set . this simulation set is primarily used for comparisons in the energetic 
 analysis portion of the results , but we also wanted to determine whether 
 the rmsd spaces explored by the inhibitor in the two approaches were 
 different . si figure s3 compares the mean 
 rmsd values of the binding region , as well as the maximum and minimum 
 values , for the cry1 and cry2 simulation sets using the same reference 
 structure for each simulation of a given inhibitor . for most of the 
 inhibitors , the mean , minimum , and maximum rmsd values are similar . 
 only in the case of the weak - binding j1 inhibitor these results do not guarantee that exactly 
 the same conformations are sampled by the two approaches , nor does 
 it indicate that the proportions of conformations sampled are similar . however , the results do indicate that long simulations do not explore 
 rmsd space that is significantly farther away ( i.e. , has higher rmsd 
 values ) than that sampled with an ensemble of diverse short simulations . it is important to note that the published crystal structure contains 
 six magnesium ions , five of which are near the binding region . the authors note three specific magnesium ions 
 that seem particularly important structurally and also note that the 
 binding affinity dramatically decreases in the absence of magnesium . we chose not to include magnesium ions in this study , with the exception 
 of the mg simulation ( table 2 ) , because previous 
 nmr data and binding assays have indicated that moderate levels of monovalent 
 salt are sufficient to stabilize the rna . to investigate whether simulations 
 with magnesium differed from those without , a single 232 ns simulation 
 was performed using the experimental coordinates for the j5r inhibitor , 
 magnesium ions , crystallographic waters , and rna ( with the necessary 
 sequence modifications to create a tetraloop consistent with the nmr 
 structure ) . in addition to the crystallographic magnesium ions , 200 
 mm kcl was included in the bulk solvent . no changes in the ligand 
 binding mode were observed upon the inclusion of coordinated magnesium 
 ions within the rna structure . si figure 
 s4 compares the regions of highest magnesium ion occupancy with the 
 crystallographic locations of the magnesium ions . although qualitatively 
 similar , none of the highest - occupancy locations observed in the simulation 
 reproduce the exact coordination contacts observed in the crystal 
 structure . the magnesium coordinated to atom n7 of a6 in the crystal 
 structure moves to coordinate the phosphate groups of g30 and c31 
 in the simulation . the magnesium coordinated to both o4 of u14 and 
 op1 of u9 in the crystal structure loses the interaction to the phosphate 
 and forms an interaction with o2 of c8 instead . finally , both magnesium 
 ions located near g5 in the crystal structure move elsewhere and are 
 replaced , to some extent , by potassium ions . these results , combined 
 with the observation of stable simulations without magnesium ions , 
 suggest that magnesium is not critical in simulations with moderately 
 high concentrations of monovalent salt ( i.e. , both had 200 mm kcl ) . this is consistent with in vitro binding assays at relatively high 
 concentrations of monovalent salt with no magnesium where high - affinity 
 binding was detected . our simulation results 
 support the crystal structure as the predominant 
 conformation in solution . however , specific interactions indicated 
 by the nmr data are inconsistent with the crystal structure , indicating 
 that under the nmr solution conditions alternative binding modes may 
 be significantly populated . in the following section , we attempt to 
 reconcile the nmr data with the key binding interactions identified 
 in the crystal structure . we were first curious whether the crystal structure 
 was actually compatible with the nmr data and therefore that the nmr 
 refinement procedure was to blame for the inconsistency in conformation . 
 one way to investigate this is to determine whether the nmr distance 
 restraints are satisfied by the crystal structure . to determine this , the r - averaged distances for all of the 
 atom pairs that were assigned noe distance restraints were computed 
 for the 200 + ns simulation of j4r from the cry1 simulation set ( this 
 inhibitor structure is identical to that of the published nmr structure ) . 
 a rotating structure with indications of noe violations greater than 
 1.0  is provided as a movie file in the si , and a list of all restraint violations is given in si table s3 . most of the 
 large violations occur between inhibitor and rna atom pairs , with 
 the largest being those between the dimethylamino arms of the inhibitor 
 and rna residues a6 , c11 , and u12 . next , we attempted to perform md 
 simulations where we gently enforced the nmr distance restraints on 
 the crystal conformation rna . inevitably , the crystal - like binding 
 pocket conformation was destroyed , including the hydrogen bonds between 
 the ligand and residue 33 . we also tested this same approach starting 
 with the conformation from the nmr1 j5s simulation ( which partially 
 converted to the crystal conformation ) . this procedure was further repeated , except 
 that all of the ligand restraints to residue a6 were removed , and 
 again the crystal - like conformation was lost . from our attempts 
 to reconcile the nmr and crystal structures , we concluded that the 
 two distance restraints between the ligand and residues c11 and u12 
 are the key driving forces in maintaining the nmr conformation during 
 the refinement procedure . to confirm this in one trial we included 
 all of the noe distance restraints , and in the other we removed the 
 two restraints between the ligand and residues c11 and u12 . these 
 simulated annealing calculations differed from the originally published 
 approach in two ways . first , in order to simplify the procedure , only 
 the distance restraints were used while torsion and rdc restraints 
 were ignored . second , in the original nmr refinement , the ligand cyclic 
 ring was protonated on the primary amine rather than the more favorable 
 secondary amine . the annealing 
 results clearly confirmed that the two restraints between the ligand 
 and residues c11 and u12 are a key determinant in the conformational 
 outcome . with all of the distance restraints included , the ligand 
 is unable to form hydrogen bonds with residue 33 , and the resulting 
 binding pocket resembles that in the published nmr structure . when 
 these two restraints are excluded , the hydrogen bonds between the 
 ligand and residue 33 are clearly formed , and the binding pocket resembles 
 that in the crystal structure . in addition to the fine details 
 of local structure , we were interested 
 in whether the crystal conformation was consistent with the rdc data 
 published with the nmr structure . to 
 test this , we performed a 5 ns simulation with rdc restraints enforced 
 on the crystal conformation . this suggests 
 that although significant noe distance violations are observed for 
 the crystal conformation , its global conformation is consistent with 
 the rdc data . finally , we were interested in the possibility 
 of a second binding 
 site . the conditions necessary to saturate the rna  ligand complex 
 and minimize exchange broadening of the nmr spectra required significant 
 excess ligand , creating the possibility of an rna with multiply bound 
 ligands . to test this , we extracted conformations of the rna  ligand 
 complex from the cry1 j4r and j4s simulations and performed docking 
 of a second j4r or j4s ligand . in all cases , the best - scoring docking 
 pose was in a pocket near residues c11 and u12 , which also happens 
 to be a magnesium ion binding site in the crystal structure . the exact 
 pose of the second ligand was dependent on the specific receptor conformation , 
 but the result raises the possibility that at high concentrations 
 a second highly charged ligand could bind in the region near residues 
 c11 and u12 and produce additional noe signals . the base conformation 
 of residue a6 and its relationship to the 
 restrained dimethylamino chain of either inhibitor j4 or j5 is another 
 key difference between the nmr and crystal structures . in the nmr 
 structure , the a6h2 proton is proximal to the dimethylamino group , 
 while in the crystal structure the a6h8 proton is proximal . a careful 
 analysis of the nmr data for selectively isotopically labeled rna j4 
 complexes appears to uphold the original nmr report and suggests that 
 there is a fundamental difference between the orientations of a6 in 
 solution and in the solid state . si figure 
 s5a shows that there are strong noes from a6h2 to the methyl protons 
 and to the methylene protons on the dihydrofuran side chain , while 
 the corresponding noes to h6h8 are weak . c hsqc data ( si figure s5b , c ) , which show no adenosine h8 
 protons resonating at 7.92 ppm . while the problematic noes to c11 
 and u12 could result from binding of a second ligand to the rna , the 
 extremely weak noe from a6h8 to the ligand is incompatible with the 
 crystal structure , suggesting that the conformation at this site under 
 solution conditions differs from that solved by crystallography . taken together , at present the data suggest the following . at least 
 two nmr distance restraints are incompatible with the crystal structure 
 and may have arisen from a second binding mode . however , the overall 
 conformation of the crystal structure is consistent with the rdc alignment 
 restraints derived from nmr analysis . the nmr noe data involving residue 
 a6 indicate that its conformation in solution differs from the form 
 captured in the crystal . this could result from the effect of magnesium , 
 a second ligand , or some combination of the two . a major goal of this study was to determine 
 whether energetic binding analyses of the simulation trajectories 
 could reproduce the experimental trends in binding energy . this is 
 a challenging problem because of the highly charged ligands and the 
 relatively tight range of kd or gbinding values estimated in the experiments . given the instability of the nmr1 simulation sets , we considered only 
 simulations performed with the crystal conformation . one of the quickest 
 methods for obtaining binding energy data from md simulations is to 
 perform mm - gbsa / mm - pbsa analyses using the single trajectory method . 
 from the mm - gbsa ( figure 7 ) and mm - pbsa ( si figure s6 ) results for both the cry1 and cry2 , 
 several observations can be made . first , the magnitudes of the binding 
 energies are significantly larger than those observed experimentally . 
 for mm - gbsa , the range of binding energies for the various inhibitors 
 is around 45 to 65 kcal / mol . for mm - pbsa , the magnitude 
 is slightly smaller , with binding energy ranges of around 30 
 to 55 kcal / mol . obviously there is a significant difference 
 between these values and the 5 to 8 kcal / mol range 
 identified experimentally ( table 1 ) . one 
 explanation is that the free energy change upon rna reorganization 
 that accompanies inhibitor binding is not taken into account in the 
 single - trajectory approach . we performed md simulations and mm - gbsa 
 analysis on the apo - rna using the published nmr structure as the initial conformation and can estimate 
 that this would add around + 10 kcal / mol to the mm - gbsa binding energies . an additional energetic component that we did not include in our mm - gbsa / mm - pbsa 
 results is the change in solute entropy upon binding . we attempted 
 to estimate the conformational entropy change for the inhibitors alone 
 ( the estimates are between 0 and 25 kcal / mol , as discussed below ) , 
 but the size of the rna presents a challenge to obtaining a full solute 
 estimate . additionally , the rotational and translational entropy loss 
 upon ligand binding would likely add in the range of 312 kcal / mol . together , these corrections are in the range of + 1347 kcal / mol 
 and could , in theory , bring at least the mm - pbsa results close to 
 the experimental ranges . however , the raw binding energies from mm - gbsa / mm - pgbsa 
 analyses have not only a large magnitude but also a large spread in 
 the values . for example , if we exclude the results for the j1 inhibitor 
 ( because of uncertainty about its exact experimental binding energy ) , 
 the range for mm - gbsa is 10 kcal / mol and for mm - pbsa around 
 20 kcal / mol . given the similarity of the ligands and their 
 binding mode with the rna , it is unlikely that the aforementioned 
 corrections for receptor reorganization and solute entropy could account 
 for the large range in values . finally , even if we ignore the incorrectly 
 large range of binding energies , the relative trend in binding energy 
 does not match what is observed experimentally ( table 1 ) . the weakest binder , j1 , can be distinguished from the other 
 inhibitors , although less so with mm - pbsa than mm - gbsa . as for the 
 other inhibitors , the computed values do not correlate well with experiment . 
 in particular , both mm - gbsa and mm - pbsa rank j2 as having a binding 
 energy equal to or lower than that of the j5 enantiomers , despite 
 the fact that experimentally j2 is among the weakest binders while 
 j5 is among the strongest binders . it is not immediately clear why 
 the mm - gbsa / mm - pbsa results do not match the experimental data , but 
 this outcome may be related to the failure of implicit solvent models 
 to model highly charged systems accurately . in fact , it is known that 
 the hawkins , cramer , and truhlar gb model does not produce accurate 
 values for salt bridges and incorrectly 
 models dna helices . mm - gbsa binding energy 
 results for the ( top ) cry1 and ( bottom ) 
 cry2 simulation sets . the mean values for the 20 individual cry2 simulations 
 are depicted by the    symbols ( bottom ) , whereas 
 the bar shows the overall mean value . all of the data are in units 
 of kcal / mol , and the specific values for both panels are given in si table s4 . stereoisomers / diastereomers are 
 grouped by color . given the poor performance 
 of these implicit models , it is reasonable 
 to seek an alternative calculation of the binding energy that explicitly 
 includes solvent contributions . inhibitor systems using 
 the total potential energy of the explicitly solvated trajectories 
 directly . this was possible because all of the systems in the cry1 
 and cry2 simulation sets , with the exception of j1 , had by design 
 identical numbers of atoms ( excluding the inhibitor atoms ) . likewise , 
 the lig simulation set , which contained just the free inhibitors and 
 solvent , also contained identical numbers of solvent atoms . since 
 the j1 inhibitor has a different charge and therefore a different 
 number of counterions , it was excluded from consideration . the relative 
 energetic contributions to binding for the trivalent inhibitors ( j2j6 ) 
 can be calculated by subtracting the average potential energy of the 
 free - inhibitor simulations ( lig set ) from the average potential energy 
 of the bound - inhibitor simulations ( cry1 or cry2 set ) . only the single 
 long simulation strategy was used for the free inhibitor simulations 
 ( lig ) because the potential energies were tightly converged . confirmation 
 of this is demonstrated by the fact that the mean potential energies 
 of all the enantiomers are within the expected error of one another 
 ( si table s5 , fourth column ) . it should be noted that the absolute value of the binding enthalpy 
 is meaningless because it represents the force field energy difference 
 of two systems with unequal atoms ( i.e. , ligand in solution and ligand 
 bound to rna ) . however , the relative enthalpy differences can be compared 
 because the relative differences between the systems are identical 
 with the exception of the inhibitor atoms of interest . several observations 
 suggest that the relative binding enthalpy values are a better predictor 
 of binding free energy than the mm - gbsa or mm - pbsa values . first , 
 the range of binding enthalpies is much smaller ( 67 
 kcal / mol ) . second , the binding enthalpies calculated from the cry1 
 and cry2 simulation sets are within the error of one another , except 
 for j6ss , which has a 0.9 kcal / mol difference between the error bounds . this was not the case for the implicit solvent approach , although 
 it should be noted that the error values are larger for calculations 
 involving explicit solvent . third , the binding enthalpies are internally 
 consistent : the four j6 diastereomers , which contain two constrained 
 rings , of which one each is present in the j3 and j5 enantiomers , 
 produce a binding trend that can be predicted from the j3 and j5 results . the binding enthalpy of the s enantiomer is preferred 
 for j5 , whereas the r enantiomer is preferred for 
 j3 . the lone exception to this internal consistency 
 is again j6ss of the cry1 simulation set ( the prediction is correct 
 in the cry2 results ) . finally , if one averages the values of the enantiomers / diastereomers 
 as an approximation of the experimental conditions ( where stereochemistry 
 was not considered ) , the trends match well with experiment ( figure 9 ) . the range of the binding enthalpies calculated 
 from simulation is somewhat larger than the range of binding energies 
 from experiment ( 4 vs 2 kcal / mol , respectively ) . however , 
 the correct trend as well as agreement between the single long trajectory 
 and multiple short trajectory approaches is observed . two completely 
 independent approaches produce very similar results , suggesting that 
 adequate sampling was obtained and that the results accurately reflect 
 the underlying force field terms used to model the system . relative binding 
 enthalpies calculated using the solvated potential 
 energy dierence ( see the text ) from the ( top ) cry1/lig and 
 ( bottom ) cry2/lig simulation sets . data for the j1 inhibitor are not 
 shown because it has a dierent charge than the rest of the 
 inhibitors , which prevents a direct comparison of the results . all 
 of the data are in units of kcal / mol , and the specic values 
 for both panels are given in si table s5 . ( top ) stereochemically averaged relative binding energies calculated 
 using the solvated potential energy difference ( see the text ) for 
 the cry1/lig ( striped bars ) and cry2/lig ( solid bars ) simulation sets . data for the j1 inhibitor are not shown because it has a different 
 charge than the rest of the inhibitors , which prevents a direct comparison 
 of the results . all of the data are in units 
 of kcal / mol , and the specific values for both the top and bottom charts 
 are given in si table s6 . color scheme 
 follows figures 7 and 8 the difference between the implicit 
 solvation and explicit solvation 
 results can be highlighted by looking at cases where the inhibitor 
 was improperly bound to the rna . as noted in methods , the initial inhibitor poses were not chosen via a docking algorithm 
 but rather by a simple rms fit of the benzimidazole core atoms . while 
 this ensured that critical hydrogen bonds were formed between the 
 inhibitor and residue g33 , it did nothing to prevent unfavorable clashes 
 of the flexible inhibitor arms with the rna . in nearly every initial 
 inhibitor conformation , potential clashes were eliminated during the 
 minimization step , and the resulting geometry resembled a reasonable 
 binding pose . however , in a few cases , one of the inhibitor arms was 
 improperly inserted through the back of the binding cavity . despite 
 the apparent strain in the geometry , some of these cases produced 
 stable simulations without integration errors or energy instability . these cases were not included in the previously discussed energy analysis , 
 but they provide a useful test case in that an accurate representation 
 of the system energy should distinguish highly strained systems from 
 those with an expected binding conformation . for the multiple short 
 simulation set cry2 , the range of mean gb total potential energies 
 of j2 in normal binding poses was 7891 to 7913 kcal / mol . the mean gb total potential energy of the incorrect strained conformation 
 was 7916 kcal / mol , suggesting that it was actually a lower - energy 
 conformation . the mm - gbsa binding energy computed using this strained 
 conformation was 77 kcal / mol , indicating stronger binding 
 by about 15 kcal / mol compared with the average mm - gbsa binding value 
 for j2 from the cry2 simulation set . apparently , some portion of the 
 implicit solvation model incorrectly modeled this interaction . in 
 contrast , using explicit - solvent potential energies accurately identifies 
 the strained j2 conformation as a high - energy outlier . the range of 
 mean explicit - solvent potential energies for j2 from the cry2 simulation 
 set was 118766 to 118790 kcal / mol , whereas the mean 
 explicit - solvent potential energy for the strained conformation was 
 118761 kcal / mol . this trend in which explicit - solvent potential 
 energies more reliably predicted strained conformations than implicit - solvent 
 potential energies was true for other cases of strained rna  inhibitor 
 conformations as well . in order to understand whether solvation 
 energies played a role 
 in the errors of the mm - gbsa / mm - pbsa results , we performed further 
 analysis of the lig simulation set ( in which the inhibitor was simulated 
 freely in explicit solvent without rna ) . after postprocessing of the 
 lig trajectories to extract the gb potential energy , a comparison 
 of the explicitly solvated and gb - solvated ligand potential energies 
 was made ( si figure s7 , top ) . as indicated 
 by the trend line fit , the relative potential energies are similar 
 whether explicit solvent energies or gb solvent energies are used . an rms fit of explicit - solvation energy values onto the implicit - solvation 
 values reveals that the differences between the two methods range 
 from 0.06 to 1.60 kcal / mol . first , this shows that the mm - gbsa values 
 are primarily enthalpic since they are similar to the explicit results , 
 which are purely enthalpic . second , it suggests that the likely origin 
 of the errors in the implicit solvation model is found not in computing 
 the solvation energy of the ligand alone but in calculating the solvation 
 energy of the ligand in complex with the rna . similarly , when the 
 solvation enthalpy is computed ( by subtracting the gas - phase inhibitor 
 energies from the solvated energies , either explicit or implicit gb ) , 
 very similar trends for the explicit and implicit solvation models 
 are observed ( si figure s7 , bottom ) . the 
 origin of the large difference ( 13 kcal / mol ) between the solvation 
 enthalpies for the j2 and j3 and those for the rest of the inhibitors 
 is unclear . the distinguishing feature of the two sets of inhibitors 
 is the absence or presence of a cyclic ring connection formed on the 
 oxygen side of the benzimidazole group . as a final piece of 
 energy analysis , we note that the experimental 
 binding energy trends suggest that the use of ring constraints to 
 reduce the flexibility improves the binding energy . thus , it is reasonable 
 to expect that entropic considerations play a role in the inhibitor 
 binding free energy . because of the difficulty in converging entropy 
 estimates for large molecules , only the inhibitor was considered in 
 the following entropy calculations , and the receptor entropy contribution 
 is assumed to be nearly identical for the inhibitors we evaluate . two methods were used : quasi - harmonic analysis and a configurational estimate based on bond , angle , and 
 torsion probability distributions . the 
 entropic energy penalties upon binding of the various inhibitors based 
 on calculations from the cry1 and lig simulation sets are given in si figure s8 . the convergence 
 plots demonstrate that the estimates are not reliably converged , even 
 at greater than 700 ns for the free ligands , and thus must be interpreted 
 very conservatively . one observation that is clear is that the quasi - harmonic 
 estimates have a much larger magnitude ( approximately 5 to 25 kcal / mol ) 
 than the conformational entropy estimates ( 0 to 5 kcal / mol ) . the large 
 values produced by the quasi - harmonic approximation are likely due 
 to an overestimation of the harmonic potential width for the free 
 inhibitors because no attempt was made to differentiate between conformational 
 energy minima , which necessarily results in a wider estimated width 
 in order to cover a broader space . in contrast , the configurational 
 estimate using probability distributions distinguishes both well - width 
 and conformational differences . in this case , the drawback is that 
 individual degrees of freedom ( bonds , angles , and torsions ) are assumed 
 to be uncorrelated when only the first - order approximation is used , 
 as was done here . it is important to note that neither of these estimation 
 techniques include the entropy loss due to changes in translational 
 and rotational entropy upon ligand binding , which can be estimated 
 to add an additional 310 kcal / mol . in view of the assumptions made in the estimates and the lack of 
 clear convergence , it is unwise to treat these values as anything 
 more than qualitative observations . given that caveat , we do find 
 that the average entropic penalty upon binding is smaller for the 
 fully constrained inhibitors ( j6 diastereomers ) than for the less 
 constrained ones ( j2 , j3 , j4 , and j5 ) . in order to test whether high - throughput 
 computational techniques could be used to predict accurate conformations 
 for binding between highly charged ligands and rna , we used dock 6.5 
 to perform docking analysis of the inhibitors shown in figure 1 on the crystal rna conformation . with the exception 
 of the weakest - binding inhibitor considered , j1 , the docking pose 
 with the best score for each inhibitor was consistent with the crystal 
 structure binding mode and formed all of the critical binding contacts 
 ( figure 10a , b ) . the j1 inhibitor binds incorrectly 
 to a major groove pocket on the opposite side of the correct rna ligand 
 pocket . it is unclear why j1 was so poorly docked , but this may reflect 
 its weak binding value . not surprisingly , the trend in docking scores 
 did not match the experimental binding energy trend ( si table s9 ) . however , the docking results suggest that fairly 
 accurate binding poses can be obtained for highly charged systems 
 at low computational cost . given the somewhat elaborate procedure 
 we used to obtain charges for the inhibitors , we were interested in 
 whether similar docking results could be obtained using a more traditional 
 approach in which the semiempirical am1-bcc charge model was used 
 to assign atomic charges to the inhibitors . the results suggest that 
 the more rigorous resp approach does yield a significant improvement 
 in the results ( figures 10c and s11 ) . this consideration is likely to be more 
 important for multivalent and highly charged ligands such as those 
 studied in this work . taken together , the results show that methods 
 such as dock 6.5 using reliable scoring functions can be applied to 
 generate reasonable binding modes . the resulting structures can then 
 be further explored using more detailed simulation and energetic analyses 
 such as those presented in this work . ( a ) overlay of the best - scoring docking poses using resp charges for 
 each of the 12 stereochemically distinct inhibitors on the crystal 
 structure receptor rna . the only inhibitor that did not bind in the 
 expected orientation was the weakest binder , j1 . ( b ) close - up view 
 of the inhibitor docking poses in the binding site ( using resp charges ) . ( c ) comparison of the am1-bcc ( black ) and resp ( red ) charge methods . the best - scoring pose for each of the 12 inhibitors is plotted with 
 its corresponding rmsd value ( using the experimental crystal structure 
 as the reference , benzimidazole core atoms only ) . the outlier data 
 point for the resp values ( red ) is for the weak - binding j1 inhibitor . given the robustness of the crystal conformation 
 across simulations of all the inhibitors in figure 1 , we were interested in how well ligands with a different 
 scaffold could be accommodated by the binding site . we designed four 
 novel ligands that both exploit the known inhibitor interactions and , 
 in certain respects , reduce the complexity of the ligand ( figure 11 ) . these novel ligands 
 have several advantages : less positive charge , fully aromatic rings 
 that are less flexible , and no chiral centers . one of these ligands , 
 n7 , was particularly stable in the receptor binding site during a 
 130 ns simulation ( si figure s12 ) . mm - gbsa 
 analysis on that trajectory yielded a binding energy of 61 
 kcal / mol , which is in the same approximate range as those of the known 
 inhibitors , although our caution regarding implicit solvent models 
 applies . additionally , the n2 and n3 ligands were moderately stable 
 and maintained most of the critical binding contacts . visual inspection 
 of the simulation trajectories suggests that further optimization 
 of the dimethylamino arm length and orientation would likely further 
 improve the binding . the 
 presented results are relevant for computational drug development 
 targeting the hcv domain ii site as well as for research on rna inhibitor complex is strongly supported over the 
 nmr conformation by simulation using the amber ff12sb force field 
 in explicit solvent . this result highlights the need for accurate 
 initial structures when performing md simulations . because of the 
 rugged energy landscape of rna , highly incorrect initial structures 
 will likely not reach the global minimum on tractable time scales . 
 since this concern is coupled with known force field flaws for rna 
 ( a subject of ongoing research in our lab ) , the force field does 
 not predict that the nmr conformations are near a local energy minimum . the result is a quick degradation of the initial structural contacts 
 during the simulation . in contrast , the crystal structure is clearly 
 in a minimum of the force field energy landscape and yields extremely 
 stable simulations that maintain fidelity in the binding region . it still remains unclear how best to resolve the apparent contradictions 
 between the nmr data and the crystal structure . the low proton density 
 of rna compared with proteins presents a challenge in the use of nmr 
 spectroscopy for high - resolution rna structure determination . the 
 addition or exclusion of just a few atom - pair restraints can make 
 a large difference in the resulting structural refinement . in this 
 case , close contacts that would result in detectable noes between 
 rna residue g33 and the inhibitor do not exist , so there was no justification 
 during the original nmr refinement to enforce the hydrogen - bonding 
 interactions observed in the crystal structure . in regard to noe distance 
 restraints that are violated by the crystal structure , possible explanations 
 include the presence of an additional inhibitor binding site near 
 the primary site or , alternatively , a third unidentified conformation 
 that is present in solution . at present , even if a simulation 
 model maintains structural fidelity , drug 
 development efforts require accurate estimates of relative binding 
 energy . traditionally , mm - gbsa / mm - pbsa trajectory postprocessing techniques 
 have been moderately successful at predicting the binding free energies 
 of protein  ligand systems , but studies with highly charged 
 ligands and highly charged receptors ( e.g. , rna ) are rare . as has 
 been noted elsewhere , the binding free 
 energy is largely determined by the difference between the desolvation 
 energy and the energy of the bound complex . for a highly charged ligand  receptor 
 interaction , both of these values will be very large , and thus , errors 
 in the method will dwarf the binding energy value . in this case , the 
 error is likely not related to insufficient sampling ( our estimated 
 errors are reasonably small ) but rather is due to an error in the 
 model used to describe the desolvation energy and the energy of the 
 bound complex . as additional evidence , we have obtained very poor 
 results from md simulations of various rna structures when using the 
 hawkins , cramer , and truhlar gb implicit solvent model ( data unpublished ) . this could explain why the binding energies calculated from the explicit - solvent 
 systems compare more favorably to experiment than the binding energies 
 obtained from implicit - solvent approximations . in view of the relative 
 success of explicit solvent in comparison to implicit solvent for 
 md simulations of rna , it is not surprising that energetic results 
 utilizing the former solvation terms would produce better results . similarly , the charge parametrization method appears to be crucial 
 for at least docking studies but likely simulations as well . this 
 is not surprising given the large net charge of the inhibitors used 
 in this study ; the careful procedure reported here provides a useful 
 framework for future charge parametrization on highly charged ligands . the docking results suggest that methods such as dock 6.5 with good 
 scoring functions can accurately predict ligand binding modes when 
 care is taken with ligand charge derivation . the data suggest 
 that the multiple short simulation approach offers 
 efficiency benefits over the single long simulation approach for energetic 
 analysis . the mean explicit - solvent binding enthalpy values computed 
 using the two approaches were within the error of each other for 10 
 out of the 11 inhibitors considered . in the case that differed , it 
 seemed likely that the multiple simulation approach was correct on 
 the basis of the argument for internal consistency of stereochemical 
 binding . the aggregate simulation time for the multiple simulation 
 approach was only 480 ns ( cry2 set ) , which is significantly less than 
 the 2616 ns of aggregate time used for the single simulation approach 
 ( cry1 set ) . the drawback to the multiple simulation approach is that 
 it does not ensure binding stability over longer time scales . the 
 2 ns trajectories used are not long enough to allow large structural 
 transitions . thus , if we had used only 2 ns trajectories , even 240 
 of them , the difference in stability between the nmr conformation 
 and the crystal conformation would not have been as clear . in view 
 of the ongoing questions about possible differences between the binding 
 modes for the crystal and solution conditions , the use of longer simulations 
 may be justified when analyzing highly charged rna  ligand complexes . finally , we conclude that computational simulations are not intended 
 to replace or  compete  with experimental observations . rather , they should aid in the evaluation of existing data and improve 
 the accuracy of hypotheses and predictions for future efforts . the 
 data presented in this work should suggest several avenues for additional 
 research , both experimental and computational .
despite the many 
 biological functions of rna , very few drugs have 
 been designed or found to target rna . here 
 we report the results of 
 molecular dynamics ( md ) simulations and binding energy analyses on 
 hepatitis c virus internal ribosome entry site ( ires ) rna in complex 
 with highly charged 2-aminobenzimidazole inhibitors . 
 initial coordinates 
 were taken from nmr and crystallography studies that had yielded different 
 binding modes . during md simulations , the rna 
 inhibitor complex 
 is stable in the crystal conformation but not in the nmr conformation . 
 additionally , we found that existing and standard md trajectory postprocessing 
 free energy methods , such as the mm - gbsa and mm - pbsa approaches available 
 in amber , seem unsuitable to properly rank the binding energies of 
 complexes between highly charged molecules . 
 a better correlation with 
 the experimental data was found using a rather simple binding enthalpy 
 calculation based on the explicitly solvated potential energies . in 
 anticipation of further growth in the use of small molecules to target 
 rna , 
 we include results addressing the impact of charge assignment 
 on docking , the structural role of magnesium in the ires 
 inhibitor 
 complex , the entropic contribution to binding energy , and simulations 
 of a plausible scaffold design for new inhibitors .