computational probe mapping is frequently 
 applied in structure - based 
 drug design ( sbdd ) to identify potential binding pockets along a protein 
 surface ( i.e. , hot spots ) . however , mapping is usually limited by the implementation 
 of gas - phase minimizations , wherein protein flexibility and solvent 
 competition are ignored .
this leads to a rugged potential surface 
 and many local energy minima , where druggable sites are indistinguishable 
 from irrelevant minima .
furthermore , predictions based on a static 
 receptor structure often neglect important binding - site features , 
 such as the presence of transient cavities or bridging water molecules . as a result , 
 appropriately modeling the binding potential continues to be a challenge 
 for sbdd .
several computational and experimental models have 
 emerged that 
 address the limitations of traditional solvent mapping .
the multiple 
 protein structure ( mps ) method addressed the issue of protein flexibility 
 by creating pharmacophore models from consensus mapping of conformational 
 ensembles with solvent probes .
while this technique has demonstrated success in 
 mapping the binding sites of pharmaceutically relevant targets , it can not account for desolvation 
 penalties or water - bridging contacts , which may be essential to accurate 
 prediction .
the multiple solvent crystal structure ( mscs ) technique uses x - ray crystallography to determine receptor structures 
 in the presence of competing water and organic solvent .
this paved 
 the way for fragment - based methods by providing the first experimental 
 identification of preferred sites for probe binding .
mscs results 
 have demonstrated high potential for use in drug development and have 
 inspired several new computational approaches .
simultaneous 
 protein flexibility and solvent mapping was first 
 implemented in 2009 when molecular dynamics ( md ) simulations were 
 used as a tool to map hot spots .
seco 
 et al . performed md simulations with seven proteins that were solvated 
 by a 20% volume / volume ( v / v ) solution of water and isopropanol ( ipa ) . 

although their method sought to detect binding sites and predict druggability , 
 it did not address aromatic hotspots and was unable to reproduce some 
 of the experimental binding sites .
yang and wang have conducted similar 
 studies that included phenol with the alcohol / water mix to address 
 aromatic interactions .
another md - based approach , site 
 identification by ligand competitive saturation ( silcs ) , was implemented 
 with a ternary - solvent box composed of 1 m benzene ( bnz ) , 1 m propane , 
 and water to locate hot spots and to reproduce experimental binding 
 sites .
the authors introduced a dummy atom with 
 a virtual repulsion term for both bnz and propane in order to solubilize 
 the hydrophobic solvents and prevent probe 
their fragment 
 maps for bnz along the trypsin surface found abundant local minima 
 that were mapped equally or better than the valid binding site .
in addition , guvench and mackerell noted that 
 the repulsive term may also prevent the mapping of secondary binding 
 sites due to the unnatural physical interactions that are enforced 
 between the hydrophobic probes .
the most 
 recent technique was developed by bakan et al . , which uses mixtures of water with multiple probes simultaneously 
 ( ipa , acetamide , acetic acid , isopropylamine ) .
although each of these 
 md - based techniques have shown some success in identifying binding 
 sites , they were limited in their ability to selectively map hot spots 
 without the use of either a weighting or artificial repulsion term . 

many spurious minima were also identified . over the past few 
 years , we have worked to develop mixed - solvent 
 molecular dynamics ( mixmd ) using amber , with an aim of solvent mapping 
 across a wide range of targets . instead of simply showing 
 that probes could occupy regions of known binding sites
, we sought 
 to establish a stronger foundation that would allow for application 
 of our method to target systems without requiring a priori knowledge of the binding sites .
our initial study with mixmd used 
 a 50% w / w solvent box composed of acetonitrile ( acn ) and water , and 
 the results showed that mixmd reproduced binding sites from mscs studies 
 with excellent convergence to the true hot spots and no spurious minima .
we were interested in extending the probe set for mixmd to 
 fully 
 enable identification of hydrogen - bonding sites , hydrophobic contacts , 
 protein  protein interactions , and aromatic pockets .
here , we 
 show what happens when poor parameters are used in mapping a protein . 

we then use md simulations of water - probe boxes to examine which remain 
 evenly mixed .
it is important to identify additional functional sites , 
 particularly aromatic sites , for mixmd to be most useful for sbdd . 

we also 
 identify additional probe solvents that expand the ability of mixmd 
 to locate important interaction profiles in protein 
ligand 
 binding : acetone ( ace ) , n - methylacetamide ( nma ) , 
 imidazole ( imi ) , and pyrimidine ( 1p3 ) .
we highlight the importance 
 of developing a solid foundation for mapping with md , particularly 
 as it relates to the validation of parameters through detailed evaluation 
 of metrics for probe  probe dispersion .
solvent probes were selected to represent 
 specific interaction types for use in mixmd ( table 1 ) .
parameters for ace , acn , and nma were obtained 
 from studies of liquid solvent with amber .
therefore , 
 we explored the suitability of opls parameters , which were developed 
 to work with tip3p water and to reproduce 
 the empirical data for pure liquids .
the 
 adaptation of nonbonded parameters from opls for use in amber simply 
 required the conversion of  to rmin , where 2  /2 = rmin .
this is necessary because boss and amber use different combining rules 
 for van der waals ( vdw ) parameters .
a comprehensive list of the parameters 
 used in our mixmd simulations is given in table s1 and figure s1 of 
 the supporting information .
pure solvent boxes of no fewer than 
 200 probe molecules were built using the tleap module 
 in ambertools .
md simulations were performed in amber10 through the sander module with 
 shake and a 1 fs time step .
each solvent box was subjected to 1000 steps 
 of steepest descent followed by 49,000 cycles of conjugate gradient 
 minimization and then heated over 20 ps from 10 to 300 k at constant 
 volume .
two nanoseconds of equilibration were followed by a 5 ns simulation 
 at constant pressure with the temperature held at 300 k by a weak - coupling 
 algorithm .
the berendsen coupling method 
 was chosen because it was used in the original publication of amber 
 parameters for the solvent molecules used in this study .
the behavior of the pure - probe boxes was used 
 to confirm proper behavior in our setup .
the final boxes of organic solvent were also used in the setup 
 of the mixed boxes .
mixed - solvent boxes were created in tleap by solvating a single probe molecule with a layer 
 of tip3p water and then a layer of the 
 equilibrated box of pure , organic solvent .
the size of the outer box 
 of probes is adjusted to achieve the desired solvent ratio .
at least 
 2500 water molecules and the equivalent mass of probe molecules necessary 
 to achieve a 50% w / w solution were used in order to ensure that realistic 
 solvent behavior would be observed ( table s3 , supporting information ) .
each layered mixed - solvent box underwent 
 the same procedure for minimization , heating , equilibration , and production 
 simulation as described for the pure - solvent boxes above .
it 
 should be noted that berendsen s coupling method can sometimes be problematic for simulations 
 of mixtures because differential heating can occur for different components 
 of the system .
our use of the method 
 with the small boxes could be considered a worst - case scenario , where 
 probes may ( or may not ) have a small bias to aggregate .
any mixture 
 of water and probe solvent that displays homogeneous behavior under 
 these circumstances is likely to be robust to a wide range of applications . 
 in our examination of thermolysin with ipa+water
, we chose the andersen 
 coupling method to be consistent with our other published applications 
 of mixmd .
there is no bias from the temperature coupling 
 that can cause the differences seen in the two parameter choices for 
 ipa when used in mixmd of thermolysin .
we characterized the 
 distribution of probe density by computing 
 atom  atom radial distribution functions ( rdfs ) with the radial 
 command in ptraj .
the final 1 ns of production time 
 was used to calculate the rdf with a bin size of 0.1  .
when 
 the solvent probes are fully solubilized into water , they are well 
 dispersed throughout the solvent box , and the rdf converges toward 
 unity at the vdw cutoff .
when phase separation occurs in simulations 
 of mixed solvent , too many probes and too little water occupy the 
 local microenvironment around each probe .
this makes the local environment 
 too dense with organic probes , and the rdf will not converge to unity 
 within the vdw cutoff .
accordingly , we used these rdf descriptions 
 of mixed vs aggregated solvent systems to evaluate the structure of 
 each mixed - solvent box . the carbon
 carbon ( rcc ) , carbon  oxygen 
 ( rco ) , carbon  nitrogen ( rcn ) , nitrogen  nitrogen ( rnn ) , nitrogen  oxygen ( rno ) , and oxygen  oxygen ( roo ) distributions 
 were assessed between probe  probe , water  water , and 
 probe  water as was relevant for each probe molecule . for brevity , 

cumulative 
 density distributions of roo for water  water 
 distances in mixed solvent were also examined to determine the number 
 density of water within the volume sphere . for the case of a 50% w / w 
 mixed - solvent box , equitable distribution of solvent probes and water 
 within the system will give a cumulative plot for roo of 
 water
the distribution of water molecules 
 relative to one another is less than in pure solvent because obviously 
 the presence of a miscible probe reduces the number of water molecules 
 that can occupy the local microenvironment around each water .
earlier papers  ours included  have 
 used the behavior 
 at the edges of the simulation box to prove even mixing of probes 
 and water .
the reasoning was that the box edges are far from the protein , 
 outside the vdw cutoff , and should have minimal bias in the ratio 
 of probe to water .
if the edges of the box had the same ratio as the 
 setup solvent box , then there was even mixing .
however , taking a step 
 back to these basic simulations of solvent boxes provided a better 
 means of measuring mixing : the use of a rdf .
rdfs of the solvents 
 relative to one another can be calculated for simulations of proteins , 
 and we recommend they be used in all mixed - solvent simulations to 
 judge whether the behavior of the probes and waters are reasonable .
the rdf , or pair correlation function g(r ) , can be used to describe the structure of the solvent 
 based on atom / molecule pairs in a system .
the function g(r ) enables a precise description of how likely 
 a probe molecule will have another probe present at a separation distance 
 of r. this function may be simplified for a liquid 
 system in a solvent box , where g(r ) is proportional to the observed number density ( o ) divided by the expected number density ( e ) .
the 
 value o is the total number of atoms at a given 
 distance , and e is the number expected at that same 
 distance if the solvent were uniformly distributed .
if the probe solvent 
 and water are evenly mixed , the rdf will converge to 1.0 at long ranges . 

if it converges to a larger number , this indicates phase separation 
 ( further discussion below ) .
previously , mixmd was validated 
 with acn as a probe , which specifies 
 hot spots for amphipathic nitrogen - containing ligands .
we applied amber parameters for acn from grabuleda 
 et al . , and they accurately reproduce 
 experimental data from mscs . for completeness , we ran simulations 
 of acn+h2o boxes to confirm that the rdf demonstrated appropriate 
 convergence to unity ( figure s2 , supporting information ) . in trying to reproduce simulations of thermolysin in ipa+h2o published by seco et al . , we 
 found that the ipa phase - separated from water ( figure 1a ) .
it is unclear whether the authors observed this same behavior , 
 and they do not note which temperature - coupling method was used .
they 
 noted  partial phase separation  in the paper , and they 
 rescaled the reference values for probe density to account for the 
 solvent behavior .
however , the bulk phase separation we observed signifies 
 the implementation of inadequate parameters for liquid ipa , which 
 were derived from the amber parameter files for threonine with charges 
 assigned after resp calculations ( ipaseco ) .
snapshots of mixmd simulations 
 of thermolysin in ipa+h2o executed using the protocol outlined 
 in seco et al .
( a ) fully flexible simulations 
 of thermolysin 
 were performed with a mixed box of 50% w / w ipaseco+h2o .
five independent runs were calculated , and every simulation 
 resulted in the solvent layers separating between 4 and 5 ns and remaining 
 separated when simulations were extended to 10 ns .
( b ) 
 the same simulation using ipaopls resulted in solvent remaining well distributed for the entire equilibration 
 and 10 ns of production in all five independent simulations .
we turned to the work of jorgensen 
 et al . for alcohol parameters 
 because they were developed to specifically reproduce a variety of 
 pure - solvent behaviors and interactions with tip3p water . for small boxes without protein , we compared 
 our simulations of ipaopls+h2o to experiments 
 by langdon and keyes that determined the density for ipa+h2o systems at different temperatures .
we found that our mixed box of ipaopls+h2o 
 reproduced the experimental density at 308 k ( 0.830g / ml ) to within 
 0.06% .
our ability to replicate this fundamental data confirmed the 
 use of good parameters for our ipa+h2o systems .
indeed , 
 mixmd simulations of thermolysin resulted in proper mixing using the 
 ipa parameters from jorgensen et al .
the rdfs for the ipa probes ( roo and rcc ) and water ( roo ) in mixed - solvent boxes further justified 
 the use of opls parameters . 

figure 2 compares the o  o rdf for both 
 mixed - solvent and pure - solvent boxes of ipa .
the roo for ( ipa - ipa)opls converged to unity with 
 an appropriately shorter peak than that of pure liquid alcohol .
however , 
 the roo for ( ipa - ipa)seco remained 
 well above 1.0 at 8  and featured a large peak for the first 
 solvent shell , the same as seen in the simulation of a pure ipa box . 

all simulations with ipaseco feature rdfs that indicate 
 aggregation of the solvent molecules . for this first example , the 
 count data behind the rdfs will be discussed in detail to further 
 explain our interpretation , which is supported by viewing the boxes .
water 
 and ( b , c ) probe  probe in a mixed - solvent environment with ipa . 

the impact of different parameters for ipa is demonstrated by the 
 poor convergence to unity for o  o in simulations of 50% w / w 
 ipaseco ( b ) when compared to the appropriate convergence 
 obtained for the same mixed solvent with ipaopls parameters 
 ( c ) .
ipa and the great majority of 
 organic solvents chosen for this 
 study are miscible with water ( table 3 ) .
the 
 definition of miscibility is that the two solutions are homogeneously 
 distributed at all ratios of the two solvents .
it might be reasonable 
 to find some bias or local structure in the first and second solvent 
 shells , but bulk separation should not be seen at long ranges .
this 
 has critical correspondence with one of the basic underlying assumptions 
 of rdfs .
rdfs are founded on a uniform distribution of solvent to 
 describe e ( solvent count / box volume ) ; any deviations 
 from uniformity ( 1.0 ) reflects the local structure of the solvent .
in all the rdfs in 
 figure 2 , water has its 
 maximum g(r ) at 2.75  , and 
 the ipa maxima are at 2.85  . for pure water at 2.75  , 
 the maximum g(r )
is 2.6 , and the 
 cumulative count of water is 1.4 molecules ( table 4 ) .
if we define the first solvent shell between 2.45 and 3.15 
  , it contains 3.9 neighboring water molecules .
when simulated 
 in a mixture with ipaopls , there are 1.1 waters within 
 2.75  and 2.9 waters within 3.15  ( table 4 ) . of course , both counts are reduced because some of the 
 interactions are occasionally fulfilled by ipa .
this is also incorporated 
 into e , where the larger box volume reduces the 
 expected density .
it is that reduced expectation ( actually 0.244 at 
 2.75  ) that results in a g(r)of 4.6 .
this is higher than the g(r ) maximum from the pure water simulation , but it does not reflect 
 more interactions with water molecules . within a 7.95  radius 
 sphere ,
the observed and expected counts are nearly equal and g(r ) approaches 1.0 , showing that the molecules 
 have a uniform distribution over larger scales .
both o and e depend 
 on the volume of the simulation box ; this drops out 
 and makes each g(r ) = no(r)/ne(r ) .
the values in this table have been rounded ; the g(r ) in the figures are the exact values . 

values for ne(r ) are based on the volume 
 of the count 
 sphere , the number of solvent in the whole box , and size of the whole 
 box during the simulation . for the mixture of water and ipaseco ,
the 
 expected count of water was 34.0 , but the observed count was 51.6 
 water molecules ( table 4 ) . clearly , there are 
 too many waters close to one another ! furthermore ,
phase separation ( as in figure 1a ) is best 
 quantified by high g(r ) at long r , but even the local changes reflected the phenomenon . 
 within 2.75 and 3.15
 , there are 1.3 and 3.4 waters ( table 4 ) , respectively , which is much closer to the values 
 shown in the pure - water boxes .
the mixed box equilibrated to a slightly 
 larger volume ( 1% change ) , which is reflected in the slight differences 
 in the expected counts between ipaseco+water and ipaopls+water .
however , the 1% change does not explain the maximum g(r ) of 5.1 , which is based on the increase 
 in number of water molecules .
the same patterns are seen in 
 the rdfs of the ipa oxygens ( table 4 ) .
pure 
 ipa simulations show appropriate behavior 
 for both sets of parameters ( black lines in figure 2bc , table s2 and text in supporting information ) .
though the rdfs may imply ipaopls has slightly tighter 
 solvent shells , they both sum to the same number of ipa molecules 
 in the first solvent shell ( 2.1 for ipaopls and 2.0 for 
 ipaseco ) . for mixed simulations ,
both ipaseco and ipaopls show lower counts than pure - solvent boxes 
 at their g(r ) maxima ( 2.85  ) 
 and within their solvent shells ( 3.95  ) .
influenced by the discrepancies in solvent behavior 
 observed for ipaseco and ipaopls , we carefully 
 examined the behavior of existing parameters for additional solvent 
 probes to determine their applicability in mixmd .
the inclusion of 
 various probe types , including complicated molecular fragments , would 
 facilitate the development of accurate pharmacophores for druggable 
 hot spots , ligand binding sites , and protein  protein interactions . 

the selection of additional probes for expanding the functional groups 
 represented in mixmd was based on existing mscs data and common interaction 
 types found in protein  ligand systems .
ipa maps both hydrogen 
 bond - accepting and -donating locations , but additional hydrogen - bonding 
 probes are needed to develop robust pharmacophore maps of binding 
 sites .
ace and nma were selected because each molecule represents 
 a different hydrogen - bonding profile .
ace locates where the protein 
 donates hydrogen bonds , and nma represents protein  protein 
 or peptide binding .
specific 
 parameters for ace and nma had been developed for use with amber and opls .
simulations of ace+h2o and nma+h2o were performed using each parameter set 
 to assess their suitability .
the g(r ) was calculated for probe  probe , probe  water , and 
 water 
analysis 
 of the resulting rdf plots indicated that in each case all g(r ) converged to unity near 8  . 

we found that both the amber and opls parameters for liquid ace and 
 nma demonstrated correct behavior according to rdf plots and visual 
 inspection ( figure 3 ) .
the o  o rdf for 
 water from the mixed - solvent simulations also illustrate that solvent 
 mixing has occurred .
analysis of the cumulative plots for the roo of water in solution with ace or nma showed 
 that the number density of water within the local volume sphere was 
 half the number density observed in pure water simulations , which 
 is the expected result for a well - distributed system of 50% w / w probe 
 and water .
visual analysis of the trajectory snapshots served to corroborate 
 the occurrence of proper solvent mixing .
however , the rdf s 
 deviation from unity at 8  was slightly higher for the opls 
 descriptions of ace and nma ( table 2 ) . for 
 this reason
, we would recommend that the amber parameters be used 
 for ace and nma in our mixmd method , especially when the amber parameters 
 are also used to describe the protein .
( a , d ) o  o and 
 ( b , e ) n  n radial distribution functions 
 for probe  probe and o  o rdfs for water  water 
 ( c , f ) in a mixed - solvent environment with ace and nma .
the impact 
 of optimized parameters from amber versus general solvent parameters 
 from opls is shown in the slightly better convergence to unity for 
 simulations using amber parameters ( d  f ) compared to opls ( a  c ) 
 parameters .
in particular , we wanted 
 to find parameters for soluble heterocycles that could map aromatic 
 hot spots without requiring an artificial repulsive term .
no other 
 md approach for probe mapping has successfully integrated aromatic 
 probes without an artificial interaction term .
futhermore , we were 
 particularly interested in developing probes that matched common heterocycles 
 used in modern pharmaceuticals . as a result ,
a highly polar 
 diazole , imi , is miscible with water and is present as a functional 
 group in a wide range of biologically active molecules , including 
 mercaptopurine ( anticancer ) , ketoconazole ( antifungal ) , and moxonidine 
 ( antihypertensive ) .
imi parameters were derived from jorgensen and 
 mcdonald , and then simulations of pure 
 imi and imi+h2o were conducted to assess the potential 
 of imi as a probe for mixmd .
the proper convergence of the rdf plots 
 to 1.0 showed that imi was properly solubilized and well distributed 
 within the binary - solvent system ( figure 4 ) . 

our pure - solvent rdfs were consistent with the results of jorgensen 
 and co - workers and further validated our parameter 
 choice and the accuracy of our protocol .
also , visualization of simulation 
 snapshots confirmed that the imi probe was well dispersed in the aqueous 
 solution .
( c ) water  water in a mixed - solvent environment with ( a ) 
 imi and ( b ) 1p3 .
opls 
 parameters for pyr , 1p2 , 1p3 , and 1p4 were each simulated in a binary 
 water  probe solution .
visualization of the simulations of pyr+h2o , 1p2+h2o , and 1p4+h2o clearly showed 
 that they underwent phase separation .
this was proved by rdfs that 
 converged to values well above 1.0 at 8  ( figure 5 ) .
however , 1p3+h2o yielded a well - dispersed system 
 with well - behaved rdf plots ( figure 4 ) .
probe  probe 
 radial distribution functions for mixed solvent 
 environments composed of ( a ) bnz+h2o , ( b ) pyr+h2o , ( c ) iph+h2o , ( d ) 1p2+h2o , and ( e ) 1p4+h2o clearly showed phase separation , with g(r ) values well above 1.0 at 8  .
we also examined two nonpolar aromatic probes , 
 bnz and iph , simply 
 to confirm that they are not appropriate probes when an artificial 
 repulsive term is not used .
comparison of the resultant rdf plots 
 to the reference rdf plot of pure water showed that the g(r ) values for bnz+h2o and iph+h2o do not converge to unity within 8  ( figure 5 ) .
in fact , the g(r ) for bnz+h2o indicated convergence to a value of 2 . the 
 simulations with iph as a probe molecule showed a converged g(r ) of approximately 1.45 at 8  . 

visualization of snapshots from these simulations showed significant 
 aggregation of the probe molecules , verifying that these probes are 
 not be soluble enough for use with mixmd without inclusion of a repulsive 
 term . because bnz is one of the solvent probes used in silcs
, 
 we were 
 able to compare the rdf results for our mixmd simulations of bnz+h2o to the published results from silcs mapping studies .
the bnz probe used in silcs was parametrized 
 based on the charmm force field and contained a repulsive term to 
 correct for probe  probe aggregation .
although a nonbonded cutoff 
 of 8  with a 58  switching term was applied in 
 silcs , the probe  probe distance for bnz converged to 1.0 at 
 13  ( solvent box size was 72   58  
  43  ) . in comparison , the soluble aromatic probes applied 
 in mixmd converged at approximately 67  .
the first solvent 
 shell of bnz in silcs was observed at 9  , while we observed 
 the first peak at 45  .
the presence of the second hydrophobic 
 probe ( propane ) in silcs simulations may have influenced these differences 
 in rdf .
however , the use of a repulsive term to prevent probe  probe 
 aggregation does not correspond with the experimental properties of 
 bnz , suggesting that silcs simulations with high concentrations of 
 bnz may yield incorrect results . to determine whether bnz and 
 iph would disperse more readily in 
 the presence of an  intermediary  probe ,
ternary mixed - solvent 
 boxes were constructed to contain 1 m aromatic probe and 1 m ipa in 
 water .
we hypothesized that the introduction of the ipa molecule would 
 enhance available interactions and aid in the dispersion of the aromatic 
 probes into water .
however , this was not observed in simulations of 
 the ternary - solvent boxes . at 5
ns of production time , the bnz molecules 
 were predominantly clustered along one side of the solvent box , while 
 the ipa molecules were dispersed throughout the water .
an additional 
 20 ns of production time did not change the observed aggregation ( figure 6 ) .
ternary - solvent simulations with iph indicated 
 a similar result ; after 5 ns of production time , iph was separated 
 with some dispersion toward the interior . elongating the run time 
 out to 25 ns
representative snapshots of md trajectory for solvent boxes of 
 1 m aromatic probe ( yellow ) and 1 m ipa ( purple ) in water ( cyan ) for 
 ( a ) bnz+ipa+h2o at 25 ns and ( b ) iph+ipa+h2o 
 at 25 ns .
although 
 all of the probes selected for 
 use in our validation study had experimental data establishing their 
 solubility in water , not all of these probes were soluble in simulation 
 ( table 3 ) .
the solvation conditions applied 
 in a mixmd simulation included high probe concentrations , which may 
 have affected probe solubility . to determine whether the outcomes 
 we observed were justified based on experimental data , we considered 
 their miscibility .
the terms soluble and miscible are frequently used 
 interchangeably ; however , they are two distinct properties .
solubility 
 is dependent upon temperature and pressure and refers to the ability 
 of a solute ( solid , liquid , or gas ) to dissolve into a solvent ( solid , 
 liquid , or gas ) , thereby forming a homogeneous mixture .
miscibility 
 refers to the ability of two liquids to form a homogeneous solution , 
 independent of proportion .
all probes that have been experimentally 
 categorized as miscible were well dispersed in our computational simulations , 
 with two exceptions ( table 3 ) .
pyridine is 
 a weak base with a nontrivial concentration of protonated pyridinium 
 ions ( pka of 5.3 ) following solubilization 
 into water .
this effect was not simulated in our studies , and we suggest 
 that the charged species may be responsible for solubilizing the uncharged 
 pyridine molecules into water . a similar explanation may also be extended 
 to pyridazine .
one advantage of using silcs as a tool for solvent 
 mapping has been that the ternary - solvent box can enable the development 
 of a complete pharmacophore model from a single md run because the 
 three probe types identify aromatic , hydrophobic , and hydrogen - bonding 
 sites .
using three of the probes with parameter 
 sets we have validated for mixmd , we constructed a ternary system 
 with imi ( aromatic probe ) , ipa ( hydrophobic and hydrogen - bonding probe ) , 
 and water at a concentration of 33% w / w by probe .
visualization of 
 the trajectory data indicated appropriate solvent mixing , and our 
 rdf results showed that both organic probes were fully dispersed into 
 solution ( figure 7 ) .
total run time required 
 a comparable amount of simulation time to a binary mixmd simulation ; 
 the total cost was an additional 4 or 8 h of production time on an 
 8-core cpu as compared to a binary ipa+h2o or imi+h2o simulation , respectively .
( a ) final snapshot from the production 
 simulation of a mixed - solvent 
 box containing imi+ipa+h2o illustrates that the probes 
 are well mixed within the system .
( b  d ) probe  probe 
 and water  water radial distribution functions for imi and ipa 
 establish convergence to unity within the trisolvent environment .
we suggest that the use of this 
 ternary mixmd model may offer several 
 advantages over the silcs approach .
we have shown that mixmd can be 
 used to identify maximally occupied sites without recovering additional 
 spurious minima , which is a feature that was not consistently seen 
 in the silcs studies .
in addition , instead 
 of requiring water to act as a hydrogen - bonding probe , in mixmd , the 
 druggable hydrogen - bonding sites can be identified in competition 
 with water using ipa as the probe .
this is important for ascertaining 
 whether a drug - like molecule that contains a similar function group 
 could displace water molecules at the proposed hydrogen - bonding site . 

a further advantage of using mixmd is the array of available probes , 
 which can be tailored to an investigator s mapping needs . for 
 example
, either imi or 1p3 could be used in the ternary box to locate 
 aromatic hotspots that are specific to the size of the ring .
the ternary 
 mixture of 33% w / w ace+acn+h2o was also simulated , and 
 the results depicted a well - distributed solvent system for use in 
 hot spot mapping with mixmd ( figure s3 , supporting 
 information ) .
we have examined solvent 
 parameters for neat liquid from the literature 
 for use in hot spot mapping through mixmd .
our work highlights the 
 importance of parameter validation and analysis of simulation data 
 when performing computational solvent mapping in order to obtain appropriate 
 behavior .
influenced by the poor results obtained using a standard 
 set of parameters used in the field for liquid isopropanol in water , 
 we identified a set of functional group probes with specific parameters 
 available for liquid simulation .
we pursued validation of these parameters 
 for neat and mixed - solvent simulations and identified six probes for 
 use with tip3p in mixmd : acetonitrile , acetone , imidazole , isopropanol , n - methylacetamide , and pyrimidine .
of course , these probes 
 should be used with more proteins beyond thermolysin ( figure 1 ) to prove their applicability .
we have previously 
 applied opls ipa parameters to also examine elastase , hewl , hewl , 
 p53 core , and rnase a. furthermore , we 
 are currently using all of these probes to map 10 additional protein 
 systems , but these are separate studies in their own right .
the extent of ideal dispersion is most easily quantified by using 
 the cumulative rdf of roo for water .
although 
 published results for pyridine and pyridazine indicated that they 
 are miscible with water , these two probes were observed to separate 
 from the aqueous phase during our simulations , rendering them unsuitable 
 for solvent mapping in our context .
we hypothesized that this disparity 
 between experiment and computation was caused by the presence of the 
 corresponding ionized species , which were not simulated . in order 
 to avoid unphysical mapping or reduced accuracy ,
our simulations did 
 not employ the use of an artificial repulsion or weighted density 
 term .
we have successfully expanded the range of probes that can be 
 incorporated into mixmd studies to allow for the mapping of druggable 
 sites , including hydrogen - bonding regions , aromatic pockets , and protein 
, a full pharmacophore 
 model could be created from the simulation of a single protein+probes+water 
 system , which would greatly reduce the computational costs of mixmd 
 studies .
this is a clear advantage of the silcs method that we would 
 like to incorporate .
ultimately , investigators could  mix and 
 match  the probes used in mixmd to identify hotspots with a 
 greater degree of specificity .
 probe mapping is a common approach 
 for identifying potential binding 
 sites in structure - based drug design ; however , it typically relies 
 on energy minimizations of probes in the gas phase and a static protein 
 structure . 
 the mixed - solvent molecular dynamics ( mixmd ) approach was 
 recently developed to account for full protein flexibility and solvation 
 effects in hot - spot mapping . 
 our first study used only acetonitrile 
 as a probe , and here , we have augmented the set of functional group 
 probes through careful testing and parameter validation . 
 a diverse 
 range of probes are needed in order to map complex binding interactions . 
 
 a small variation in probe parameters can adversely effect mixed - solvent 
 behavior , which we highlight with isopropanol . 
 we tested 11 solvents 
 to identify six with appropriate behavior in tip3p water to use as 
 organic probes in the mixmd method . 
 in addition to acetonitrile and 
 isopropanol , we have identified acetone , n - methylacetamide , 
 imidazole , and pyrimidine . 
 these probe solvents will enable mixmd 
 studies to recover hydrogen - bonding sites , hydrophobic pockets , protein  
 protein 
 interactions , and aromatic hotspots . also , we show that ternary - solvent 
 systems can be incorporated within a single simulation . 
 importantly , 
 these binary and ternary solvents do not require artificial repulsion 
 terms like other methods . within merely 5 ns , layered solvent boxes 
 become evenly mixed for soluble probes . 
 we used radial distribution 
 functions to evaluate solvent behavior , determine adequate mixing , 
 and confirm the absence of phase separation . 
 we recommend that radial 
 distribution functions should be used to assess adequate sampling 
 in all mixed - solvent techniques rather than the current practice of 
 examining the solvent ratios at the edges of the solvent box . 