chemical potential 
 (  ) describes the equilibrium movement 
 of particles between two phases or states .
the driving force behind 
 this movement comes from ( 1 ) a concentration gradient 
particles 
 tend to move from high to low concentration regions to gain mixing 
 entropy ; ( 2 ) chemical affinity 
excess chemical 
 potential ( ex ) is the quasistatic work to bring 
 a particle ( e.g. , solute molecule ) from the gas phase to the solvent ; 
 ex=   id , where 
  and id are the chemical potential and the 
 ideal gas chemical potential of the solute , respectively . in the context of statistical mechanics ,
chemical 
 potential allows for the thermodynamic state of a system to be defined 
 in terms of a grand canonical ( gc ) ensemble ( vt ) that allows 
 for variation in the species concentrations across phases / states . 

simulation procedures have long evolved toward efficiently determining 
 gibbs free energy of hydration ( hfe ) , chemical affinity , and 
 other thermodynamically relevant properties of waters and other small 
 solute molecules from gc ensembles instead of the more conventional isothermal , isobaric ( npt ) , canonical ( nvt ) , or microcanonical ( nve ) ensembles where the concentration of the species is fixed . to date , 
 many of the gc ensemble strategies have employed monte carlo ( mc ) 
 simulations toward either driving the sampling of waters or individual 
 small molecules around proteins or crystal environments , or alternately 
 improving the accuracy of the relative hfe calculations in free energy 
 perturbation ( fep ) calculations .
however , 
 since gcmc simulations of systems containing explicit solvent to represent 
 the bulk - phase suffer from convergence problems due to low acceptance 
 rates encountered for the solute insertions simulations in the past were restricted toward investigation of 
 the chemical affinities of only the solvent , drive individual solute sampling in the 
 absence of explicit solvent , determine thermodynamic 
 properties in crystal conditions , or use 
 expanded ensemble strategies . in 
 the context of protein and macromolecular environments , chemical 
 fragment sampling simulation techniques
have been employed in the 
 past toward discovery or rational design of molecules that can bind 
 to macromolecular targets with high affinity to achieve a desired 
 biological outcome .
the site identification by ligand competitive saturation ( silcs ) 
 method is one such technique that identifies the location and approximate 
 affinities of different functional groups on a target macromolecular 
 surface by performing molecular dynamics ( md ) simulations of the target 
 in an aqueous solution of solute molecules representative of different 
 chemical fragments .
however , these md isothermal , isobaric ( npt ) ensembles 
 suffer from the long diffusion time scales of the solutes through 
 explicit solvent and macromolecule environments , especially when the 
 macromolecular binding sites are deeply buried and inaccessible to 
 the solvent ( i.e. , occluded ) .
sampling of the distribution of chemically 
 diverse solutes in proteins in gc ensembles with mc simulations offers 
 the potential to overcome these limitations , yielding more accurate 
 solute spatial distributions in aqueous macromolecular environments . in this paper , we present a methodology for solute sampling in 
 explicit solvent aqueous systems and solvated protein environments 
 using a variation of gcmc in which the ex values 
 are adjusted to attain a target concentration along with md simulations 
 in an iterative fashion ( oscillating-ex gcmc - md ) . 
 in this approach gcmc sampling is performed on both the solutes and 
 the waters .
a short time - scale md after the gcmc allows for both conformational 
 sampling of the solutes and configurational sampling of the aqueous 
 system and the macromolecule . to achieve satisfactory convergence
, 
 the process of gcmc - md is repeated through multiple iterations , with 
 the ex of the species being studied systematically 
 oscillated over the iterations to drive the solute and water exchanges . 

two types of aqueous systems were considered : ( 1 ) a system containing 
 1 m of only one type of solute in water , thereby replicating the standard 
 state of the solute , and ( 2 ) a dilute aqueous mixture containing 0.25 
 m of many types of solutes .
both systems are implemented in aqueous 
 solution alone and in the presence of the l99a mutant of t4-lysozyme 
 ( t4-l99a ) .
gcmc was run on the water and all the different solutes 
 in the system with the variation of their respective ex values through the gcmc - md iterations used to improve the solute 
 exchange probabilities while allowing for determination of the respective 
 ex values required to maintain a defined concentration 
 of the solutes and waters in the aqueous systems .
as the oscillating-ex gcmc - md approach involves variation of ex such that it is not formally the gc ensemble , we show that it produces 
 a representative ensemble of conformations as the average ex values approximate the hfe of solutes at a specified target 
 concentration in aqueous systems . in the presence of the protein , 

this oscillating-ex gcmc - md strategy could be used 
 for efficiently sampling the spatial distribution of the solutes in 
 and around the protein .
multiple solute sampling in the context of 
 silcs allows for determination of their affinity patterns to the protein 
 that can be used toward rational drug design .
there are four possible gcmc 
 moves on a molecule , m ( i.e. , solute or water ) : insertion  brings 
 m into the system a from the reservoir ; deletion  removes m 
 from the system a and moves it back into the reservoir ; translation 
 and rotation  m is translated / rotated within a subvolume surrounding 
 the original location of m in system a. the probabilities of these 
 moves as governed by the metropolis criteria are1where b = ex + ln n , n = , ex is the 
 excess chemical 
 potential , n is the expected number of molecules , 
  is the density , v is the volume 
 of system a , fn is the 
 fractional volume of the subspace where the insertion attempts are 
 made , e is the change in energy due to a move , 
  is 1/kbt , kb is the boltzmann constant , and t is temperature ( 300 k in the present study ) . through the gcmc simulation , 
 the volume of the simulation system a , the total energy , and the total 
 number of particles between the system a and its reservoir
as seen in eq 1 , 
 the target concentration 
 of each solute ( n ) , the interactions of each 
 solute with all the other molecules in the simulation system a ( e ) , and the supplied ex determines the 
 solute populations through gcmc simulations . since , the move probabilities 
 of the individual solutes or waters are driven by their n and ex values , when the ex of 
 a solute is less than the work needed to move a molecule from the 
 gas - phase reservoir to the system a , a decrease in concentration of 
 the solute from the system a would occur .
likewise , when the ex supplied is more than the needed work , an increase in solute 
 concentration will occur in system a. consequently , the value of ex supplied to the gcmc simulation may be oscillated based 
 on the concentration in the simulation system and the target n. with n used as a 
 target for the oscillating-ex gcmc calculation , 
 the simulation system a is defined to 
 contain water at a bulk - phase concentration of 55 m. through the iterative 
 oscillating-ex gcmc - md simulations , the concentration 
 of the solutes and the water in system a vs their target concentration 
 ( n ) is used as a guide to vary their respective 
 ex through each subsequent gcmc iteration .
the simulations 
 start with the ex of the solutes and the water set 
 to 0 in the first iteration .
over every subsequent iteration or subset 
 of iterations , ex of the respective solutes and 
 water is varied by a magnitude that is governed by the deviation of 
 the solutes and water in the system a from their respective target n. as the concentrations of the solutes and the waters 
 reach their target , the width of the variation of each ex is decreased and this defines the onset of convergence . 

thus , in standard state simulations with only one type of solute at 
 1 m in water , the approximate hfe associated with those conditions 
 can be calculated via the average ex . upon convergence , 

the oscillation of ex values for waters and the 
 solutes is continued to obtain the average ex and 
 to facilitate the sampling under equilibrium conditions .
the 
 standard state simulations with only one type of solute at 
 1 m in system a will be referred to as scheme i and the multiple solute 
 aqueous mixture simulations as scheme ii .
system a is a spherical 
 region of radius , ra , into which gcmc 
 moves are performed . as shown in figure 1 , 
 separate reservoirs for the solutes and the waters are coupled to 
 system a. system a is immersed in a larger system , b , which includes 
 additional water . for the bulk aqueous systems ,
system b is a larger 
 sphere of radius rb = ra + dr ( in the present study , dr = 5  ) .
the larger system b limits edge effects , 
 such as hydrophobic solutes from occupying the edge of the system 
 a ( supporting information ( si ) , figure 
 s1 ) .
alternatively , system b may be treated using periodic boundary 
 conditions and/or include other chemical entities in addition to water . 

an example of a periodic system b in which a spherical system a is 
 immersed , with the combined system containing the t4-l99a is shown 
 in figure s6 of the si .
setup for the standard 
 state ( scheme i ) and the aqueous solute 
 mixture ( scheme ii ) oscillating-ex gcmc - md simulations . 

water and the solute molecules are exchanged between their respective 
 reservoirs and the spherical simulation systems , a , indicated by the 
 dashed boundary , defined by radius ra . 

system a is immersed in a larger system , b defined by the solid boundary , 
 that included additional waters of radius rb set to ra+ 5  in the present study 
 to prevent hydrophobic solutes from occupying the edge of the system 
 a ( see figure s1 , si ) . alternatively
, the 
 larger system in which system a is immersed may be treated using periodic 
 boundary conditions and include other chemical entities in addition 
 to water ( see figure s5 , si ) .
the iterative oscillating-ex gcmc - md 
 procedure 
 is described as follows:(a)run i steps of gcmc 
 to exchange solutes and waters between their respective reservoirs 
 and the simulation system a. the i steps are divided 
 between each of the m solutes and water . in the present study i/2 and ( i)/(m + 1 )
gcmc moves ( insertion / deletion / translation / rotation ) 
 are attempted for the solutes ( f1 , f1 , f2 , ... , fm ) and the waters , in schemes i and ii , 
 respectively , though the number of moves for the different solutes 
 and water are not required to be equal .
the order in which the four 
 possible gcmc moves are attempted , and the molecule ( solutes or water ) 
 on which the move is performed is randomized . in the first iteration , 
 ex of the solutes and the water
the 
 radii of the water and the solute sphere(s ) subjected to gcmc moves 
 are set to ra through the gcmc process , 
 though the energetic interactions associated with the moves also includes 
 contributions from any waters , or other chemical entities , outside 
 of the gcmc sphere that is defined as system b.(b)after the gcmc , j steps of md are run on combined system a and b. for the finite spherical 
 systems , the solutes are retained within the spherical dimensions 
 of system a , ra through harmonic flat - bottom 
 spherical restraints , while the water 
 and any other molecules in system b encompassing system a are subjected 
 to ( i ) harmonic flat - bottom spherical restraints with rb = ra + dr when system b is spherical or ( ii ) periodic boundary conditions.(c)vary the value of ex of solute and water , p ( p = f1 , f2 , ... , fm solutes or w , water ) by dp .
the magnitude of dp is determined by the 
 deviation of the current concentration , np ( p = f1 , f2 , ... , fm , or w ) in 
 system a from the target concentration ( nptarget ) .whereelse2with the variation of dp system - dependent , various schemes ( standard state , aqueous 
 mixture , 
 heterogeneous systems ) can be generated using eq 2 .
we note that different schemes may be applied to vary dp.(d)perform k iterations 
 of the oscillating-ex gcmc - md ( steps a 
c ) 
 with new values of ex from c ) . run i steps of gcmc 
 to exchange solutes and waters between their respective reservoirs 
 and the simulation system a. the i steps are divided 
 between each of the m solutes and water . in the present study i/2 and ( i)/(m + 1 ) gcmc moves ( insertion / deletion / translation / rotation ) 
 are attempted for the solutes ( f1 , f1 , f2 , ... , fm ) and the waters , in schemes i and ii , 
 respectively , though the number of moves for the different solutes 
 and water are not required to be equal .
the order in which the four 
 possible gcmc moves are attempted , and the molecule ( solutes or water ) 
 on which the move is performed is randomized . in the first iteration , 
 ex of the solutes and the water
the 
 radii of the water and the solute sphere(s ) subjected to gcmc moves 
 are set to ra through the gcmc process , 
 though the energetic interactions associated with the moves also includes 
 contributions from any waters , or other chemical entities , outside 
 of the gcmc sphere that is defined as system b. after the gcmc , j steps of md are run on combined system a and b. for the finite spherical 
 systems , the solutes are retained within the spherical dimensions 
 of system a , ra through harmonic flat - bottom 
 spherical restraints , while the water 
 and any other molecules in system b encompassing system a are subjected 
 to ( i ) harmonic flat - bottom spherical restraints with rb = ra + dr when system b is spherical or ( ii ) periodic boundary conditions .
vary the value of ex of solute and water , p ( p = f1 , f2 , ... , fm solutes or w , water ) by dp .
the magnitude of dp is determined by the 
 deviation of the current concentration , np ( p = f1 , f2 , ... , fm , or w ) in 
 system a from the target concentration ( nptarget ) .whereelse2with the variation of dp system - dependent , various schemes ( standard state , aqueous 
 mixture , 
 heterogeneous systems ) can be generated using eq 2 .
the gcmc portion of the simulations were run using an in - house 
 c++ code that implemented the grid - based gcmc scheme with the cavity - bias algorithm to drive solute and water exchanges between their reservoirs and 
 the aqueous systems , a. the system a is divided into a 1  grid 
 lattice .
when system a contains macromolecules , all the atoms of the 
 macromolecules are considered to be of the same size , and each atom 
 is assigned to a lattice site .
fragments within the system a are not 
 assigned lattice sites . after deleting the lattice sites assigned 
 to the macromolecule , a list of available lattice sites is maintained . 

the term fn in eq 1 is calculated as the ratio between the total number 
 of lattice sites , to the total number of available lattice sites . 
 during a gc move , for example , an insertion of a solute , one of the 
 available lattice sites is randomly chosen .
the solute molecule is 
 moved from its gas - phase reservoir such that the center of mass of 
 the solute occupies the randomly selected lattice site .
interaction 
 energy of the solute is then calculated with all of the macromolecule , 
 solutes and water atoms in the system . the updated interaction energy 

if the insertion is accepted , the solute is now a part 
 of the gc grid .
when the insertion is rejected , the solute is taken 
 back to the gas - phase reservoir .
in general , for every gc move , the 
 interaction energy of the selected molecule with every other solute 
 and solvent molecules and the macromolecule is calculated .
e , due to the attempted move is then applied to eq 1 to determine if the move is accepted or rejected .
solute empirical force - field parameters were obtained from the 
 cgenff and the tip3p water model was 
 used .
the studied solute molecules were 
 chosen to represent different chemical functionalities including apolar 
 benzene and propane , neutral polar molecules such as acetaldehyde , 
 methanol and formamide , and the negative and positively charged acetate 
 and methylammonium , respectively .
these solutes were those used in 
 our recent work on an extension of the silcs method .
as detailed previously , to prevent aggregates 
 of hydrophobic and charged solutes , thereby promoting faster convergence , 
 a repulsive energy term was introduced only between benzene : benzene , 
 benzene : propane , propane : propane , acetate : acetate , acetate : methylammonium , 
 and methylammonium : methylammonium molecular pairs .
this was achieved 
 by adding a massless particle to the center of mass of benzene and 
 the central carbon of propane , acetate , and methylammonium .
each such 
 particle does not interact with any other atoms in the system but 
 with other particles on the hydrophobic or charged molecules through 
 the lennard - jones ( lj ) force field term with parameters (  = 0.01 kcal / mol ; rmin = 12.0  ) .
all lj and coulomb interactions were calculated 
 during gcmc ( i.e. , no truncation of nonbonded interactions ) , including 
 interactions with system b. simulations are initiated with an empty 
 system a. the randomized gcmc process with the solutes and the waters 
 is initially run in multiples of 50  000 moves until the waters 
 in the system reach the bulk 55 m concentration . at each cycle in 
 this repetition of the 50  000 gcmc moves , the ex of the solutes and the water is increased by 1 kcal / mol to accelerate 
 the water and solute accumulation in system a. during this process , 
 the concentration of the solutes may also increase beyond their target 
 values .
once the bulk water concentration is attained , this is used 
 as the starting configuration for the oscillating-ex gcmc - md procedure described in steps a  d . after about 50 
 iterations of the oscillating-ex gcmc - md
in which 
 the ex values are varied , the concentration of the 
 solutes approach their target values . in each iteration ,
50  000 
 and 100  000 gcmc moves were run for schemes i and ii , respectively .
the gromacs package ( version 5.0 ) 
 was used for all md simulations . for the aqueous systems , during each 
 iteration the combined system a and b , including all the solutes and 
 the waters , were simulated with the leapfrog integrator ( gromacs integrator 
  md  ) , with a 2 fs time - step , at 300 k through a nose - hoover 
 thermostat .
the last conformation of the 
 gcmc run was used as the starting conformation for the md .
the system 
 was initially equilibrated over a period of 100 ps with reassignment 
 of velocities .
the 
 last conformation from the production md is used as the starting conformation 
 of the next mc run .
the lincs algorithm was used to constrain 
 water geometries and all covalent bonds involving a hydrogen atom . 

van der waals ( vdw ) and electrostatic interactions were switched off 
 smoothly over the range of 810  .
solutes and waters 
 were held within the spherical dimensions of system a or b by applying 
 harmonic flat - bottom restraints with 
 a force constant of 1.2 kcal / mol on the following 
 solute or water atoms : the massless particle at the geometric center 
 of benzene , propane , acetate , and methylammonium ; the carbon atoms 
 of the acetaldehyde , methanol , and formamide ; and the oxygen atom 
 of water . for the bulk aqueous systems ,
the radius , rb of the harmonic flat - bottom restraints used to define 
 system b applied to only the water was 5  larger than the radius ra = 20  defining the restraint applied 
 to the solutes in system a. for t4-l99a , pdb coordinates 181l 
 with the benzene ligand was used 
 following the deletion of the ligand .
the t4-l99a structure was inserted 
 in boxes replicating schemes i and ii systems with 1 m of benzene 
 and 0.25 m each of the different solutes ( benzene , propane , acetaldehyde , 
 methanol , formamide , acetate , and methylammonium ) , respectively .
akin 
 to the aqueous systems , 10 such boxes were built for scheme i with 
 1 m benzene and scheme ii with the multiple solutes , with the solutes 
 randomly inserted in each of the boxes .
these systems were minimized 
 for 1000 steps with the steepest descent algorithm in the presence of periodic boundary conditions ( pbc ) .
the leapfrog version of the verlet integrator with a time step of 2 fs was used for heating 
 and equilibration .
long range electrostatic interactions were handled 
 with the particle - mesh ewald method with 
 a real space cutoff of 12  , a switching function was applied to the lennard - jones interactions 
 at 12  , and a long - range isotropic correction was applied to the pressure for lennard - jones interactions 
 beyond the 12  cutoff length . during the minimization and equilibration 
 harmonic positional restraints with a force constant of 2.4 kcal / mol
were applied to protein non - hydrogen atoms . for the md in 
 the iterative protocol
, the position restraints were removed and replaced 
 by weak restraints only on the backbone c carbon atoms 
 with a force constant ( k in 1/2 kx ) of 0.12 kcal / mol .
this was done to prevent the rotation of the protein in 
 the simulation box and potential denaturation due to the presence 
 of a highly concentrated fragment solution .
the oscillating-ex gcmc - md iterations are 
 repeated 
 until convergence , typically 200 iterations for the aqueous systems , 
 yielding 100 ns of md and 10 and 20 million mc steps for schemes i 
 and ii , respectively .
furthermore , to ensure sufficient sampling and 
 convergence , 10 separate oscillating-ex gcmc - md 
 simulations are run for scheme ii and for each solute in scheme i 
 resulting in a cumulative 1 s ( 10  100 ns ) of md for 
 scheme ii and each solute in scheme i. runs differ through a randomly 
 generated seed for both gcmc and the leapfrog md integrator operations 
 in each iteration .
oscillating-ex gcmc - md of the 
 protein systems are repeated over 100 iterations , yielding a cumulative 
 500 ns of md over the 10 separate simulations from both scheme i and 
 scheme ii .
the goal of the present work 
 was the implementation of an algorithm 
 to perform equilibrium studies of solutes in aqueous solution , including 
 mixtures of solutes and in the presence of a macromolecule , with t4-l99a 
 used as the macromolecular model system in the present study . in the 
 context of scheme i ,
if the oscillating-ex gcmc - md 
 scheme is producing an approximately correct ensemble , the average 
 ex should approximate the hfe corresponding to a 
 1 m standard state aqueous system and the binding affinity of solutes 
 to the protein .
with the solute mixture in scheme ii the approach 
 would allow for determining the ex required to drive 
 sampling of the distribution of solute molecules in a heterogeneous 
 aqueous system and determine solute affinity patterns around the protein 
 site .
there are other methods available that determine ex required to maintain a target concentration of the water 
 molecules using the gcmc methodology , including one that relies on more complex proportional  integral  derivative 
 ( pid ) controller scheme , and continuously fluctuate the ex at every step of the gcmc . in comparison , the modulation 
 of the excess chemical potential in the present work is driven solely 
 by the concentrations of the solutes and the water and their target 
 concentrations or chemical potentials . in the initial calculations , 
 when the solute and solvent in the 
 finite spherical systems shared the same spherical boundary , the nonpolar 
 molecules sampled the surface of the spherical system ( si , figure s1 ) .
this is due to the favorable 
 interactions of the nonpolar solutes with the nonpolar , vacuum environment 
 outside of the spherical simulation system . to overcome this ,
the 
 radius of the simulation system for water was increased by 5  , 
 such that the  pancaking  of nonpolar molecules is avoided 
 as they stay fully hydrated .
the use of such an extended system will 
 also alleviate other potential edge effects that the other solutes 
 may encounter .
we note that the larger system b can extend to any 
 distance beyond the solute restraints defining system a , including 
 being modeled using pbc , as was performed in our calculations that 
 included t4-l99a . in the remainder of this section results using 
 both scheme i and 
 ii will be presented for aqueous systems alone , which will be followed 
 by scheme i and ii calculations on a system that contains t4-l99a . 

this lysozyme mutant was selected as a model system as it has been 
 widely used in experimental and computational studies of ligand binding 
 as well as for studies of the impact of mutations on protein structure 
 and stability .
the l99a mutation in the c - terminal domain creates a completely 
 buried , hydrophobic cavity of 150  that , 
 although inaccessible in static structures , binds small hydrophobic 
 ligands in a rapid and reversible manner . such an occluded 
 pocket offers a rigorous test of the sampling effectiveness of the 
 presented oscillating-ex gcmc - md methodology , including 
 a quantitative evaluation of the approach . for both scheme i and aqueous 
 mixture scheme ii
the concentration is the target for the variations 
 of the ex through the oscillating-ex gcmc - md iterations , with the ex of both the solute 
 and the water initially set to 0 .
the concentration is converted to 
 a number nptarget associated with the volume of system 
 a for each solute type or water and applied in eq 2 . in both cases ,
the average value of ex of both the solutes and the water converged close to their hfe and 
 the sampling of ex values approximate gaussian distributions 
 ( si , figure s2 ) .
table 1 lists both the calculated average ex values 
 and the experimental hfe ( hfe ) of the solutes .
the ex values were obtained from block averages 
 over the final five 10 cycle blocks obtained between cycles 150 and 
 200 , with the final ex calculated as the average 
 of these block averages along with the associated standard errors . 

hfes were also calculated using fep method 
 by inserting a solute molecule in a 27  cubic box with bulk 
 water and 1 m of the solute ( ie .
nine additional solute molecules 
 in addition to the solute being perturbed ) , following the simulation 
 methodology as detailed in our previous works , to account for possible limitations in the force field 
 that would yield a hfe in disagreement with the experimental data .
hfe calculated using 
 fep calculations of 1 m solutions using a previously described methodology .
figure 2 traces the progression of both 
 the concentration and ex through the gcmc - md iterations , 
 as the ex was being varied based on nptarget in 
 system a for both schemes i and ii .
both the concentration and the 
 ex at each iteration in figure 2 are presented as the average over the 10 independent simulations . 
 in both schemes
i and ii , the solutes and the solvent of system a 
 attain their target nptarget values , corresponding to 1 and 
 0.25 m of solute in schemes i and ii , respectively , and 55 m of water . 

however , for the charged fragments the acceptance rates for particle 
 insertions were low , consistent with previous studies , due to the unfavorable electrostatic interactions .
the convergence 
 took longer for these cases and simulations were hence run for 200 
 iterations .
plots of ex and concentrations of the solutes 
 and the water as a function of the oscillating-ex gcmc - md iterations .
the ex value for the solutes 
 and the water was varied every iteration based on their respective 
 concentrations in the simulation system a. the target concentration 
 of the solutes was 1 and 0.25 m in schemes i ( black ) and ii ( red ) , 
 respectively , while water was maintained at bulk concentration of 
 55 m in both the systems .
solutes are benzene ( benz ) , propane ( prpn ) , 
 acetate ( acet ) , methylammonium ( mamm ) , methanol ( meoh ) , formamide 
 ( form ) , and acetaldehyde ( aald ) .
table 1 shows the average ex and the concentration over the 10 simulations , from the 
 final 50 
 iterations .
the average ex and the concentration 
 were largely similar to those measured across iterations 50200 
 ( si , table s2 ) , consistent with the onset 
 of convergence from iteration 50 . for the hydrophobic and polar molecules , 

the largest deviation occurs with formamide with the values 
 falling between the hfe and hfe values . 

the deviation of ex and hfe from the 
 hfe for the charged fragments , acetate and methylammonium , 
 is due to the vacuum - to - solvent interface potential not being accounted 
 for in the present calculations . for 
 monovalent anions
/ cations , this contribution was calculated to be 
 about 12.5 kcal / mol , respectively , with the tip3p water model .
further , as gcmc methods with charged solutes 
 are limited by low acceptance rates for particle insertions , the agreement 
 of the ex for the charged systems was expected to 
 be limited , which is most notable with acetate .
overall , these results 
 establish that the presented oscillating-ex gcmc - md 
 methodology approximates the hfe of organic solutes via their ex in aqueous systems .
in addition the approach is indicated 
 to be suitable for more complex aqueous mixtures as evidenced by the 
 ex values obtained from scheme ii being in satisfactory 
 agreement with the experimental and hfe values . while the oscillating-ex gcmc - md simulation protocol 
 achieved the correct concentration and ex ,
the ability 
 of the method to obtain the correct spatial sampling of the solutes 
 in the finite spherical systems is important for studies of heterogeneous 
 systems .
spatial sampling was investigated via the analysis of radial 
 distribution functions ( rdf ) .
rdfs of selected solute atoms and the 
 water oxygens were calculated from the cumulative md sampling of both 
 the scheme i and ii oscillating-ex gcmc - md simulations . 

these were compared to the rdfs obtained from explicit - water 15 ns 
 pbc md simulations that maintain scheme 
 i and ii concentrations of the solutes and water in a cubic box with 
 a side of 50  and includes the explicit treatment of long - range 
 nonbond interactions via the particle mesh ewald and isotropic lj correction methods ( see si , section s1 for simulations details ) .
the rdfs from the finite spherical systems oscillating-ex gcmc - md and the explicit - water pbc md match very well ( si , figure s3 ) . this indicates the efficiency 
 of the use of oscillating-ex values to drive gcmc 
 sampling and that the treatment of the long - range nonbond interactions 
 in the presented oscillating-ex gcmc - md protocol 
 does not significantly impact the spatial sampling of the aqueous 
 systems .
some fluctuations are seen in the rdf from scheme ii pbc , 
 which are due to sampling issues : the md only simulations were run 
 for 15 ns versus a total of 100 ns md in the oscillating-ex gcmc - md protocol .
thus , the present methodology attains 
 spatial sampling consistent with that observed in unbiased md pbc 
 simulations .
we note that the ex settled close to 
 the hfe in a scheme i simulation of acetaldehyde and methanol without 
 md at the end of every iteration ( si , table 
 s3 ) .
however , since the molecules are rigid during each cycle of gcmc , 
 md is likely needed to preserve the correct conformational and spatial 
 sampling of the solutes in these bulk - phase environments ( see below ) .
subsequent calculations focused on determining if the gcmc - md approach 
 could obtain equilibrium solute sampling with known , fixed ex values . with that objective
, we ran a set of gcmc - md simulations , 
 where instead of starting with setting ex to 0 , 
 it was held fixed at the hfe throughout the iterations . shown in figure 3
are the concentration and ex for 
 the acetaldehyde system as a function of the oscillating-ex gcmc - md and fixed-ex gcmc - md iterations . 
 in the fixed-ex simulations
, we found that as the 
 number of iterations increased and waters attained bulk concentrations , 
 the number of solute exchanges decreased considerably .
similar trends 
 were seen for the other solutes ( figure s4 and s5 of the si ) .
when the moves are performed for one molecule , 
 the other molecules ( both solutes and waters ) are stationary and participate 
 only in nonbond interactions with the current molecule being exchanged , 
 thereby preventing overlapping moves of molecules into locations already 
 occupied .
thus , at the 55 m bulk - phase concentration of water , it 
 is easier for the smaller waters to fill up cavities in system a than 
 the larger solutes , leading to a drastic decrease in the gcmc exchange 
 probabilities for the solutes ( figure 3c and 
 d ) . as continuous gcmc exchanges of fragments through insertions and 
 deletions , and local relaxation through translations and rotations 
 are important to maintain chemical equilibrium between system a and 
 the coupled gas - phase fragment reservoirs this outcome was deemed 
 problematic .
( a ) concentration and ( b ) ex of acetaldehyde 
 and the average probabilities of ( c ) insertion + deletion ( pins+del ) and ( d ) translation + rotation ( ptrans+rot ) as a function of the number of oscillating-ex gcmc - md iterations from the scheme i and ii oscillating-ex gcmc - md of acetaldehyde with ex fixed 
 at the hfe ( black , green ) or fluctuated by dex ( red , 
 blue ) , respectively .
note that the number of solute exchanges with 
 the gas - phase reservoir are greater with the fluctuating ex while the average concentration and ex are approximately equivalent .
to overcome this , ex of both the solutes 
 and 
 the water were cyclically fluctuated around their respective hfes . 
 dp = np / nptarget was alternately added and subtracted to p ( p = f1 , f2 , ... , fm , w ) 
 following every three iterations of gcmc .
thus , the frequency of oscillation 
 is six cycles , and by maintaining this oscillation through the length 
 of the simulation , the average ex remains constant 
 at the values in table 1 .
we note that by periodically 
 fluctuating ex the system is likely not a formal 
 gc ensemble ; however , by maintaining the average ex constant over the length of the simulation the extent of deviation 
 is minimal , which is supported by the agreement between the oscillating-ex gcmc - md calculated average ex and fep 
 hydration free energies ( table 1 ) .
such 
 variations in ex lead to improvements in 
 sampling over the course of the gcmc - md iterations , as evidenced by 
 both the change in the number of the solute molecules and the increase 
 in the gcmc move probabilities ( pins+del and ptrans+rot for the insertions and 
 deletions and translations and rotations , respectively ) , as shown 
 in figure 3 .
similar trends were observed for 
 the apolar and the other polar solutes ( si , figure s4 , s5 ) .
thus , together with gcmc of both the solvent and 
 the solute , it is important to continue to fluctuate the ex supplied to these molecules once the target concentration , 
 as described in eq 2 , is attained to maintain 
 efficient sampling of the solutes in the simulation system a. computational chemical fragment 
 mapping studies to date have applied md simulations such 
 that the proteins that could be studied had to have binding sites 
 accessible to the bulk solvent , allowing the solutes and water to 
 diffuse in and out of the sites on the time scale of the simulations . 
 however , a larger number of biologically important proteins , such 
 as the g - protein coupled receptors ( gpcrs ) and nuclear receptors as well as other 
 macromolecules , have very deep or occluded pockets that would require 
 simulations of extensive duration . to overcome this accessibility 
 limitation ,
the present oscillating-ex gcmc - md method 
 offers great potential and we selected the well - studied t4-l99a , which 
 contains an engineered occluded binding site for benzene as a model system .
the oscillating-ex gcmc - md sampling method was identical to the aqueous systems 
 above with the exception that system b was treated as periodic with 
 the solute being studied included in system b at the target concentration 
 used to drive the gcmc sampling .
system a was a 20  sphere centered 
 on the t4-l99a binding site defined by residues ala 99 and met 102 
 ( si , figure s6 ) .
scheme i calculations 
 on t4-l99a involved only benzene as a solute at 1 m along with 55 
 m water .
this allowed for validation of the method to drive sampling 
 of benzene in the occluded binding pocket of the protein as well as 
 yield a quantitative estimate of ligand binding . across a 10  
 37.5
ns oscillating-ex gcmc - md simulation with conformations 
 saved every 10 ps , benzene was bound to the t4-l99a binding site for 
 a total of 372 ns with the pocket being empty for only 3 ns .
the average 
 benzene concentration in the simulation system a was 1.5  0.2 
 m. these results indicate that the oscillating-ex gcmc - md method can sample the occluded pocket as well as attain 
 the defined concentration in the entire system that drives the gcmc 
 sampling .
the pocket sampling also allows for estimation of 
 the binding affinity 
 of a ligand g , using3where , r is the gas constant 
 ( kcal / molk ) , t ( k ) is the temperature , [ pl ] 
 is the concentration of the bound ligand , [ l ] is the total concentration 
 of the ligand , [ p ] is the concentration of the protein , and vref is the reference volume in concentration 
 units ( 1660  per one ligand molecule for 
 1 m of ligand ) .
since the simulations 
 are maintained in equilibrium , the bound vs unbound ligand concentrations 
 can be correlated to the time fraction of ligand bound vs unbound 
 in the binding site through the simulation .
a previous technique to 
 determine binding affinity of small molecules to the protein site 
 using the gcmc methodology was limited to gas - phase simulations , using 
 generalized born / solvent accessibility ( gb / sa ) solvation calculations 
 to convert the gas - phase energies to aqueous energies . by running oscillating-ex gcmc - md 
 of both the solutes and the waters , binding affinity calculation explicitly 
 incorporates desolvation contribution of both the binding pocket and 
 of the ligand .
as the t4-l99a pocket is completely occluded , 
 the presence or the 
 absence of a solute atom at the active site is driven only by the 
 gcmc insertion / deletion moves . over the 10  37.5 ns oscillating-ex gcmc - md simulation that involves a total of 18.75 million 
 gcmc insertion / deletion attempts for the benzene , with [ l ]  
 1.5 m and a [ pl]/[p ] ratio of 99.4/0.6 ,
although 
 there is some difference from the experimental binding affinity of 
 5.19  0.16 kcal / mol , it 
 should be noted that an increase of 10 steps out of the 
 1.9  10 total steps with benzene in the pocket translates 
 to nearly 1 kcal / mol difference in the g 
 calculated .
this emphasizes the difficulty of converging the calculation 
 of a binding constant to a single site in a protein , although force 
 field effects could also impact the obtained g. scheme ii calculations involved the seven solutes along 
 with water . 
 as with the aqueous system ,
the target concentration for the solutions 
 was 0.25 m. this simulation was run for a total of 10  50 ns . 
 to facilitate analysis of the results , affinity patterns of the selected 
 atoms from the different solutes in the occluded binding pocket , called 
  grid free energy ( gfe ) fragmaps  ,
gfe fragmaps are boltzmann transformed 
 probability distributions for the solute atoms that are normalized 
 with the distributions of these molecules in aqueous solution in the 
 absence of the macromolecule .
this normalization accounts for the 
 free energy penalty of solute desolvation when calculating the gfes 
 ( see si , section s2 ) .
these maps may then 
 be visualized to qualitatively evaluate the ability of the oscillating-ex gcmc - md sampling method to reproduce the positions of different 
 ligands known to bind t4-l99a that have been subjected to experimental 
 analysis .
nine ligands were considered : ( 1 ) benzene , ( 2 ) o - xylene , ( 3 ) p - xylene , ( 4 ) ethylbenzene , ( 5 ) benzofuran , 
 ( 6 ) indene , ( 7 ) indole , ( 8) isobutylbenzene , and ( 9 ) n - butylbenzene .
figure 4 shows the aromatic 
 ( benc ) , aliphatic 
 ( prpc ) , and polar hydrogen bond donor ( meoh , forh ) and acceptor ( meoo , 
 foro , aalo ) maps along with the crystallographic orientations of the 
 ligands .
benzene occupies the aromatic fragmap ( purple ) while the 
 aliphatic moieties of the other ligands occupy the aliphatic fragmap 
 region ( green ) that protrudes away from the benzene molecule .
in addition , 
 neutral h - bond donor and acceptor maps are in the binding pocket and 
 are in the vicinity of the corresponding functional groups on benzofuran 
 and indole .
importantly , although we use the protein structure from 
 the t4-l99a  benzene complex ( pdb 181l ) as the starting conformation , the oscillating-ex gcmc - md simulations correctly identify the ability of the 
 pocket to alter its conformation to allow favorable interaction with 
 the aliphatic moieties as well as its ability to accommodate polar 
 functionality .
this is an outcome of the inclusion of md , in which 
 protein flexibility is included , in the presented methodology .
when a second scheme ii simulation was run without the md , the pocket 
 in the neighborhood of ser117 and leu118 was not available for sampling 
 by the apolar fragments ( si , figure s7 ) . 

it is possible that md sampling could be redundant when system a does 
 not contain any macromolecule , since the ex calculated 
 both with and without the md are similar ( si , table s3 ) .
however , to maintain consistency between the macromolecular 
 and the aqueous environment only sampling , md is retained during the 
 latter .
selected gfe fragmaps at the ligand binding site of the t4-lysozyme 
 l99a from a 10  50 ns oscillating-ex gcmc - md 
 simulation and the minimized crystal conformations of the 9 ligands 
 ( see text for the list ) ; protein atoms occluding the view of the binding 
 pocket were removed for clear visualization .
fragmaps are displayed 
 at a cutoff of 1.2 kcal / mol for the benc ( purple ) , prpc ( green ) 
 and 0.5 kcal / mol for aalo ( pink ) , meoo ( dark red ) , foro ( red ) , 
 meoh ( blue ) , forh ( light blue ) .
the use of gfe fragmaps also has the advantage that it allows 
 for 
 quantitative evaluation of the relative affinity of the ligands , based 
 on ligand gfe ( lgfe ) scores , as previously described .
lgfe quantifies the overlap of the atoms in the ligand 
 with the corresponding gfe fragmaps .
protein 
 orientations generated using ( i ) md sampling of the ligands bound 
 to the protein and ( ii ) mc sampling of the ligands in the field of 
 the fragmaps ( see si , section s2 ) .
as shown 
 in figure 5 , both lgfe and lgfe correlate very well with the experimental binding free energies 
 ( high r and predictive index ( pi ) ) .
importantly , the lgfes can distinguish between 
 the binding activity of both congeneric series and diverse classes 
 of ligands .
we note that the range of experimental binding affinities 
 is 2.1 kcal / mol while the lgfe scores from the protein + ligand 
 md ensemble and the mc sampling are spread over wider ranges of 4.4 
 and 6.1 kcal / mol , respectively .
this is not unexpected as 
 the lgfe scores are not true free energies of binding as numerous 
 terms that contribute binding are omitted ( e.g. , 
 the configurational entropy of the ligands ) .
correlation 
 of the experimental binding affinity g 
 ( from ref ( 37 ) ) with 
 the lgfe scores for the nine ligands considered ( see 
 text for the list ) .
the lgfe scores are obtained from md and mc conformational 
 ensembles of the ligands ( lgfe , lgfe ) and 
 the gfe fragmaps ( see main text and si , 
 section s2 ) .
overall maps have a very good correlation ( high r and predictive index , pi ) .
our novel oscillating-ex gcmc - md 
 methodology 
 therefore allows for investigations of the ex of 
 solutes in aqueous solution , including solutions containing multiple 
 solutes .
central to the approach is the use of variable ex throughout the gcmc cycles of the simulations .
this leads 
 to convergence of ex for given solute(s ) and environment 
 based on the target concentration and the maintenance of solute sampling 
 once the system has converged with respect to the ex or target concentration .
along with probing the ex required to maintain the solutes at target concentrations in aqueous 
 environments , the oscillating-ex gcmc - md approach 
 is also useful toward efficient solute sampling .
notably , the method 
 may be used to effectively sample the configurational space of an 
 occluded pocket in a macromolecule ; the engineered pocket in t4-l99a 
 in the present study .
this sampling may be performed with a complex 
 mixture of solutes to map the functional preference of the binding 
 pocket . when this is analyzed in the context of the silcs methodology
, 
 the approach is shown to qualitatively reproduce the binding orientation 
 of known ligands as well as quantitatively rank order the binding 
 affinity of the ligands .
it is anticipated that this capability will 
 be of utility for the application of computer aided drug design methods 
 to macromolecules with deep or occluded binding pockets such as those 
 in protein nuclear factors and gpcrs .
 solute 
 sampling of explicit bulk - phase aqueous environments in 
 grand canonical ( gc ) ensemble simulations suffer from poor convergence 
 due to low insertion probabilities of the solutes . to address this 
 , 
 we developed an iterative procedure involving grand canonical - like 
 monte carlo ( gcmc ) and molecular dynamics ( md ) simulations . 
 each iteration 
 involves gcmc of both the solutes and water followed by md , with the 
 excess chemical potential ( ex ) of both the solute 
 and the water oscillated to attain their target concentrations in 
 the simulation system . by periodically varying 
 the ex of the water and solutes over the gcmc - md iterations , solute exchange 
 probabilities and the spatial distributions of the solutes improved . 
 
 the utility of the oscillating-ex gcmc - md method 
 is indicated by its ability to approximate the hydration free energy 
 ( hfe ) of the individual solutes in aqueous solution as well as in 
 dilute aqueous mixtures of multiple solutes . for seven organic solutes : 
 
 benzene , propane , acetaldehyde , methanol , formamide , acetate , and 
 methylammonium , the average ex of the solutes and 
 the water converged close to their respective hfes in both 1 m standard 
 state and dilute aqueous mixture systems . 
 the oscillating-ex gcmc methodology is also able to drive solute sampling in 
 proteins in aqueous environments as shown using the occluded binding 
 pocket of the t4 lysozyme l99a mutant as a model system . 
 the approach 
 was shown to satisfactorily reproduce the free energy of binding of 
 benzene as well as sample the functional group requirements of the 
 occluded pocket consistent with the crystal structures of known ligands 
 bound to the l99a mutant as well as their relative binding affinities . 
Introduction
Method and Simulations
Details
Results and Discussion