simulating a heterogeneous solvent environment , 
 which is especially 
 important for rna , adds complexity to simulations .
intramolecular 
 interactions between dynamic rna secondary structure and tertiary 
 structure elements are typically dependent on solution conditions 
 for their formation and stability .
most 
 commonly , divalent mg ions are required for these structures 
 to form , either to effectively shield the high charge density of proximal 
 phosphates or to facilitate or stabilize a conformation required for 
 the interaction .
both secondary structure 
 elements and tertiary structure contacts can be strongly stabilized 
 by mg .
traditionally , simulations 
 have had trouble describing the highly 
 polarizable nature of divalent ions , which is not explicitly captured 
 in a classical , fixed charge model .
polarizable 
 force fields have shown mixed results for the parametrization of nucleic 
 acids , though recent parameters using 
 drude oscillator approaches to address the balance between nucleic 
 acids , monovalent ions , and solvent have shown improved results .
mg ions and their solution effects have been approximated 
 using many methods , including continuum methods , reduced 
 models , and multisite models .
reduced models have been shown to reproduce 
 both the diffuse and chelated / associated ion effects when combined 
 with explicit mg .
multisite models , which 
 shift partial charges onto external sites , have also been used to 
 capture the interactions of associated mg ions with nucleic 
 acids .
though all of these models can reproduce 
 various types of experimental data , a more fully applicable point - charge 
 model requires the ability to capture the dynamics of chelated , associated , 
 and also bulk mg effects . recently , four reparametrizations 
 for fixed charge mg have been released . in the allner 

modifications , mg parameters are optimized to 
 fit the experimental exchange rate of 
 the first shell water molecules by modifying the repulsive lennard - jones 
 ( lj ) term ( villa ) .
modifications are fit to the experimental hydration free energies 
 and the ion  oxygen distance , and then are combined into a consensus 
 set ( merz ) . though parametrized independently , 
 the resulting villa and merz 12 - 6 lj parameters are somewhat similar , 
 increasing the distance of the minimum ( rmin ) and making the well 
 depth ( epsilon ) more shallow compared to the widely used aqvist parameters 
 ( table 1 ) .
li and merz also developed a divalent model 
 based on a 12 - 6 - 4 lj functional , adding a 1/r term to capture charge induced dipole effects ( 12 - 6 - 4 ) .
panteva et al . found that these parameters resulted 
 in an overestimation of the binding free energy of mg to nonbridging phosphate oxygen atoms and an underestimation of 
 the binding free energy of mg to the n7 atom of adenosine 
 and guanine bases .
they developed modified 12 - 6 - 4 
 parameters ( m12 - 6 - 4 ) which apply different polarizabilities for these 
 atom types , tuning the interactions to experimental site specific 
 binding free energies .
following the procedure of joung 
 and cheatham , each of the preceding 
 different ion parameter sets was generated for use with specific water 
 models .
in the 12 - 6 - 4 models , this is 
 achieved by scaling the c4 terms of each atom type by the ion  water 
 oxygen atom interaction ( c4(h2o ) ) and the polarizability 
 of the specific water model , as follows:1for the simulations in this work , the c4 for 
 mg in tip3p water is 132.9 (kcal / mol ) .
the lj parameters for mg models 
 used in this work , 
 and their calculated self c12 , c6 , and c4 interactions , are reproduced 
 in table 1 .
a more 
 detailed table with the lj matrix calculated for the 12 - 6 - 4 and m12 - 6 - 4 
 mg ion interactions with all atom types , showing rmin , 
 epsilon , and c12 , c6 , and c4 coefficients , is given in supporting information tables s1 and s2 .
the dynamics and structure of the varkud satellite ribozyme stem - loop 
 v rna ( slv ) have been shown experimentally to be mg ion 
 dependent . in the absence of mg ,
the five nucleotide u - turn loop exhibits considerable dynamics 
 and adopts a wide range of structures . when mg is present , the loop is stabilized in a canonical u - turn 
 conformation bymg ions , as shown in figure 1 .
the goal of this work is to determine whether mg ion 
 models can accurately describe the dynamic structure shift from a 
 magnesium free ( mgfree ) conformation to the magnesium bound ( mgbound ) 
 conformation . using molecular dynamics simulations
, we examine the 
 differences between ion models after the rna has transitioned to the 
 mgbound form and then for all conformations sampled in the simulations . 
 by using 10 sets of independent 11.5 s md simulations 
 with varying starting ion positions for each mg ion parameter 
 set , we compare ion binding for the five models .
we find that all 
 models promote the mgfree to mgbound transition in slv rna .
both sets 
 of 12 - 6 - 4 potential mg ions trend toward more occupied 
 interactions due to loss of inner shell water molecules in all simulations , 
 while the villa mg ions lost inner shell water molecules 
 in several simulations .
interestingly , the merz and aqvist mg ions give remarkably similar results as far as occupancies 
 in the folded state and the ensemble of folded and unfolded structures , 
 despite their differing lj parameters . for future simulations of slv 

rna , merz 12 - 6 mg ions will be used , but more validation 
 against experimental observables for other rna systems should be performed . slv rna 
 in the absence and presence of 
 mg .
slv rna in the absence of mg adopts 
 a loose u - turn structure ( left ; pdb i d , 1tbk ) , which shifts conformation to a canonical 
 u - turn structure upon mg binding ( right ; pdb i d , 1yn2 ) .
the positions 
 of the mg ions are based on nmr experiments and are labeled 
 14 .
the slv sequence and numbering 
 used in our simulations are shown at the bottom of the image .
the 
 hairpin loop is formed by the u6:a12 noncanonical loop - closing base 
 pair , with u7-g8-a9-c10 forming a u - turn , extruding u11 .
experimentally 
 determined hexahydrated mn binding sites are shown as 
 green spheres , numbered 14 .
the coordinates of mgfree slv rna were 
 taken from the first model 
 ( labeled as the best representative conformer ) in the 1tbk nmr ensemble .
the rna was neutralized with 8 mg ions and solvated with 3967 tip3p water molecules .
an additional four mg ions , five k ions , and 13 cl ions were added to generate approximately 
 40 mm mgcl2 and 50 mm kcl concentrations .
consensus set parameters 
 to generate the  merz  set , allner et al .
parameters to generate the  villa  set , li and merz 12 - 6 - 4 parameters to generate the 
  12 - 6 - 4  set , and panteva 
 et al . m12 - 6 - 4 parameters to generate the  m12 - 6 - 4  
 set .
amber formatted prepin and frcmod 
 files are included in the supporting information for using these parameters in amber14 ( files s1s4 ) .
additionally , 
 the add12_6_4 command in parmed.py was used to generate 
 an amber parmtop with the additional c4 terms for the 12 - 6 - 4 potential 
 simulations .
parmed.py was also used to create new atom types for 
 the n7 atoms of adenosine and guanosine ( nd and ng , respectively ) 
 and the nonbridging phosphate oxygen atoms ( o2 ) for the m12 - 6 - 4 model . 

initial ion positions were randomized 6.0  from the rna solute 
 and 4.0  from each other , using 10 different seeds for each 
 model .
this resulted in 10 starting structures with randomized ion 
 environments for each ion model tested .
simulations for the 
 12 - 6 potential systems ( aqvist , merz , and villa ) 
 were carried out using the amber14 suite of programs gpu ( cuda ) version 
 of pmemd with spfp ( a mixed single / fixed precision model on the gpu ) .
simulations for the 12 - 6 - 4 potential systems were performed using 
 the amber14 cpu version of pmemd , with the lj1264 = 1 flag set in 
 the input .
the ff12 force field for nucleic acid simulations was used , 
 combining ff99 + parmbsc0 modifications + chiol3 modifications for rna .
the particle mesh ewald method , with default parameters , 
 was used to handle electrostatic interactions with a 9  cutoff 
 for direct interactions .
bonds to hydrogen 
 were constrained using the shake algorithm , allowing use of a 2 fs 
 time step .
hydrogen mass repartitioning 
 was employed for the 12 - 6 - 4 potential simulations performed on the 
 cpu , allowing a 4 fs time step . for these 
 simulations
the masses of atoms bonded to hydrogen were decreased 
 by 3.024 da using parmed.py , and the hydrogens mass was increased 
 by the same amount .
each of the 50 total systems was minimized 
 with 500 steps steepest 
 descent and 500 steps conjugant gradient with 25 kcal mol restraints on the rna , followed by heating to 300 k over 100 ps , 
 maintaining constant volume .
temperature was regulated using the weak - coupling 
 algorithm , and a coupling time constant of 0.2 ps was used .
after heating , each system was subjected to rounds 
 of 500 steps steepest descent , 500 steps conjugate gradient minimization , 
 and 50 ps constant pressure equilibration with decreasing positional 
 restraint weights ( 5.0 , 4.0 , 3.0 , 2.0 , and 1.0 kcal mol ) . during these rounds of constant pressure equilibration , the pressure 
 relaxation time
a final equilibration was run for 
 500 ps at constant pressure with 0.5 kcal mol restraints , 
 after which restraints were removed and a 50 ps unrestrained constant 
 pressure simulation was run , with temperature coupling time constant 
 set to 5 ps and the pressure relaxation time set to 5 ps . production 
 runs began for each ion set from the 10 equilibrated 
 systems , and the properties of the rna and ions were measured from 
 this portion of the trajectory .
temperature was set to 300 k and regulated using 
 the langevin algorithm with a collision frequency of 5 ps and random seeds used for initial velocity assignments , preventing 
 synchronization artifacts .
h - bond analysis 
 was performed using a distance cutoff of 4.0  between mg ions and all atoms of the rna .
all folded rna structures 
 from all ion model runs using either a 12 - 6 lj potential or a 12 - 6 - 4 
 lj potential were used to generate an average structure for grid analysis . 

grids were normalized to a mg ion density of 9.7  
 10 mol /  .
previous simulations were performed 
 on slv rna , starting from the mgbound and mgfree forms , in the presence 
 of either 50 mm nacl or 40 mm mgcl2 .
these simulations showed that when monovalent ions are 
 present , the slv rna folds to a mgbound - like structure with an rmsd 
 of 2.0  from the mgbound loop structure , irrespective 
 of the starting conformation .
gatekeeper  
 structure as a precursor to mg ion binding , upon which 
 the rmsd to the mgbound loop shifts to a lower 1.625  value 
 and the u - turn motif is stabilized via mg ion association . 

loop ) is not 
 observed in simulations with monovalent salt , either in short simulations 
 ( hundreds of nanoseconds ) or during replica exchange simulations ( 57.6 
 s of aggregate simulation time ) . in the presence of divalent ions ,
the slv rna sampled low rmsds 
 to the mgbound reference in the simulations starting from either mgbound 
 or mgfree structures . even in short hundred nanosecond time scale 
 simulations ,
the mgfree loop folded to the mgbound state upon the 
 addition of mgcl2 . in an effort to better understand the 
 role of mg in loop folding ,
ten copies of the equilibrated system , 
 where the mg ions added were parameters from aqvist , merz , villa , 12 - 6 - 4 , or m12 - 6 - 4 , were each run for 11.5 s .
the 
 addition of mg ions to the system should result in a 
 switch from the mgfree rna structure to the mgbound rna structure .
figure 2 shows the histogram 
 of loop rmsd to the folded mgbound structure and shows that very low 
 rmsds are achieved in simulations performed in each ion model .
simulations 
 sample structures which are considered folded ( those below the green 
 line indicating an rmsd of 1.625  ) in the presence of all mg ion models .
the mgbound structure is sampled in 9/10 aqvist 
 simulations , 8/10 merz simulations , 6/10 villa simulations , 8/10 12 - 6 - 4 
 simulations , and 8/10 m12 - 6 - 4 simulations .
although not all of the 
 structures have converted to the mgbound structure in the 11.5 
 s of md simulation time , the majority of simulations in each 
 case have transitioned .
given this , the simulations were not extended 
 further ; although longer simulations may show transitioning of all 
 structures , full convergence of the structural distributions would 
 likely require application of enhanced sampling methods , such as replica 
 exchange .
mgbound structure for each of 10 simulations containingthe following : 
 ( a ) aqvist mg ions , ( b ) merz mg ions , ( c ) 
 villa mg ions , ( d ) 12 - 6 - 4 mg ions , and ( e ) 
 m12 - 6 - 4 mg ions .
time dependent data with 1000 step running 
 averages are shown in supporting information figure s1 . in all plots
, the green line represents the 1.625  
 cutoff below which structures are considered  folded  
 or mgbound . to understand the influence of 
 each ion model on the bound structure , we first combined all structures 
 ( for 10 simulations using each ion model ) below a cutoff of 1.625 
  rmsd to the mgbound loop , described previously as the cutoff 
 between gatekeeper and mgbound ( or folded ) structures .
the low rmsd - to - mgbound structures were analyzed 
 using the hbond command in cpptraj to determine the fraction occupancies 
 of mg ions to all rna atoms .
average structures were 
 generated from the three 12 - 6 lj parameter models , aqvist , merz , and 
 villa , and from the 12 - 6 - 4 lj potential models , 12 - 6 - 4 and m12 - 6 - 4 , 
 and mg ion occupancies were mapped onto them . when 
 the rna is folded in the same conformation
, the 12 - 6 potential ion 
 models generally behave in the same manner .
the similarity of these 
 interactions can be seen in figure 3a  c , which shows the average folded mgbound 
 structure and localized mg densities ( 500 greater 
 than bulk mg ion densities , or 10% occupancy from hbond ) . 
 the sum of occupancies for the atoms in each residue of slv rna ( 117 ) 
 are shown below the slv rna for each ion model .
the profiles reinforce 
 the similarity between the 12 - 6 models , especially the aqvist and 
 merz models . in general , the highest occupancies center on phosphate 
 oxygen atoms . specifically , the u - turn nucleotides u7-g8-a9 coordinate 
 a mg ion , locking in the folded structure , differentiating 
 the folded structure from the rest of the ensemble .
another area of 
 high mg ion density is between c10 and u11 phosphate 
 oxygen atoms .
mg ions can be seen associating with rna 
 atoms in the major groove and at phosphate oxygen atoms when interactions 
 populated at 10% occupancy are considered ( figure s2 ) .
top : mg ion density grids at > 10% occupancy ( 250500 
 bulk occupancy ) for slv rna atoms for each ion model : ( a ) aqvist , 
 ( b ) merz , ( c ) villa , ( d ) 12 - 6 - 4 , and ( e ) m12 - 6 - 4 .
bottom : sum of fraction 
 occupancies for all residues ( 117 ) for all systems ( aqvist , 
 merz , villa , 12 - 6 - 4 , and m12 - 6 - 4 ) , where stem residues 15 
 and 1317 are shown in dark gray , loop - closing base pair residues 
 6 and 12 are shown in light gray , and loop residues 711 are 
 colored to match the bases in figure 1 and panel a. an interesting difference between ion models association 
 with the mgbound rna is the high occupancy in the 12 - 6 - 4 lj potential 
 models , shown in figure 3d , e . though the occupancy is mapped onto slv rna at the same 500 
 bulk concentration as the 12 - 6 models , there are many more highly 
 occupied sites .
the sum of the fraction occupancies per residue shown 
 below the slv rna reinforces that there is , in general , more highly 
 occupied sites than the 126 models .
this is due to the mg ions described by the 12 - 6 - 4 parameters and , to a lesser 
 extent , the m12 - 6 - 4 parameters losing an inner shell water molecule 
 and forming a chelated interaction with the rna , but there is no obvious 
 structural explanation for this chelation .
sites for mg ion association have been experimentally 
 determined by nmr chemical shift perturbation using manganese hexahydrate 
 and cobalt hexammine and are mapped onto 
 the slv structure in figure 1 .
simulations provide those parameters 
 for each mg ion association site , and table 2 shows the results for each 
 ion model .
the results show the distribution of interactions across 
 the rna atoms is almost identical for each ion model which uses a 
 12 - 6 lj potential .
the trend for more and less occupied binding site 
 interactions holds for the 12 - 6 - 4 lj potential models , but overall 
 the occupancies are higher .
this site was determined based on mn(h2o)6 proximity to the rna , which could 
 have formed d - orbital interactions with site 2 atoms , accounting for 
 the discrepancy between simulation and experiment .
however , the m12 - 6 - 4 
 simulations , which explicitly address the balance between mg interactions with n7 atoms , have a higher occupancy of g8@n7 , indicating 
 a deficit in the 12 - 6 lj potential models .
additionally , in the 12 - 6 
 lj potential simulations , k ions have a higher occupancy 
 than mg for site 2 atoms , shown in table s3 .
figure 4 shows the 
 aggregate data for all 10 simulations performed for each ion model , 
 separated into interactions between mg and phosphate 
 oxygen atoms , base oxygen and nitrogen atoms , and sugar oxygen atoms . 

it is apparent that each ion model prefers associated interactions 
 with phosphate oxygen atoms , since these interactions range from 10 
 to 70% occupied .
the highest occupancy sites in all 12 - 6 potential 
 models are in the five residue loop or its closing base pair , where 
 the phosphate backbone populates an unusual conformation involving 
 residues u7 and a9 .
the 12 - 6 - 4 simulations show very high occupancies 
 for phosphate interactions , indiscriminately interacting with slv 
 rna throughout the combined data set , instead of showing specific 
 localization in the u - turn region .
this effect is reduced when the 
 m12 - 6 - 4 lj potential is used , but the resulting occupancies for the 
 m12 - 6 - 4 model are not as low as in the 12 - 6 models .
fraction occupancies 
 by residue ( 117 ) for mg to rna interactions for 
 all mg ion models ( left to 
 right : aqvist , merz , villa , 12 - 6 - 4 , and m12 - 6 - 4 ) .
top : phosphate atom 
 interactions ( including o3 , o5 , op1 , and op2 ) .
middle : 
 base atom interactions ( including n9 , n7 , o6 , n1 , n2 , n3 , n4 , o2 , 
 n6 , and o4 ) . bottom : sugar atom interactions ( including o4 
 and o2 ) . stem residues 15 and 1317 are shown 
 in dark gray , loop - closing base pair residues 6 and 12 are shown in 
 light gray , and loop residues 711 are colored to match the 
 bases in figures 1 and 3a .
interactions with base atoms range from 0 to 25% occupancy , 
 with 
 the most occupied interaction occurring with the u6 base o4 atom in 
 all ion models .
u6:a12 is a noncanonical loop - closing base pair , and 
 the u6@o4 atom could be associated with an mg ion .
this 
 agrees with the experimental observation that this base pair is transiently 
 formed based on the very weak experimental noe in the presence of 
 mg .
the increased occupancy for base atoms in the 12 - 6 - 4 
 model , and particularly in the m12 - 6 - 4 model , arises from interactions 
 with the n7 atoms of the adenosine and guanosine bases in the loop 
 ( residues g8 and a9 ) and indicates that the polarizability introduced 
 by the 12 - 6 - 4 model , and its subsequent tuning in the m12 - 6 - 4 model , 
 is important to adequately describe these interactions .
the 
 sugar atoms show very low occupancy , typically less than 2% 
 for all models .
the outlier occupancy of residue u7 in the m12 - 6 - 4 
 model arises due to the proximity of mg interacting with 
 the op1 phosphate of its neighboring residue g8 .
it is interesting 
 to note that the 12 - 6 - 4 lj parameter models trend toward higher occupancies 
 for all interaction types , and though this effect is ameliorated by 
 using the m12 - 6 - 4 corrections , those occupancies are still higher 
 than the 126 potential models .
the high occupancies 
 for the 12 - 6 - 4 lj potential models , and their 
 larger error bars , result from chelated interactions with the rna 
 atoms .
direct chelation involves the hexahydrated mg ion 
 losing an inner shell water molecule in favor of binding to an rna 
 atom , and these events are summarized in table 3 .
the villa simulations are the only 12 - 6 
 potential model where chelation is observed , indicated by the high 
 occupancy or high variability in the phosphate interaction with residues 
 a9 , a12 , and g16 in figure 4 .
a single simulation shows a chelated interaction between 
 an mg ion and the op1 atom of residues a9 and g16 .
half 
 of the simulations contain a chelated interaction between an mg and op1 of residue a12 . in all 10 simulations using both 
 the 12 - 6 - 4 and m12 - 6 - 4 lj parameters an mg ion loses 
 an inner shell water molecule to form at least one chelated interaction 
 to a phosphate oxygen atom . every residue s phosphate group 
 forms a chelated interaction with the 12 - 6 - 4 and m12 - 6 - 4 mg ions during the combined set of simulations . on average ,
eight mg ions are chelated to rna atoms during a simulation with 
 12 - 6 - 4 lj parameters , and this effect is decreased to four chelated 
 mg ions during a simulation using the m12 - 6 - 4 lj parameters . n.f . =
when the rna has already folded to the mgbound structure , 
 the chelating 
 ion has no effect on the rna , as shown in figure 5a .
however , in some cases the chelation occurs 
 before the loop folds to the mgbound structure , and the mg 
rna interaction locks the loop backbone conformation , prohibiting 
 reorganization to the correct u - turn fold , depicted in figure 5b .
once directly chelated with 
 rna , it would be difficult on these time scales to see the loss of 
 this interaction , and we do not ; as with the expected direct chelation 
 binding affinity , likely simulations on the order of 10s of milliseconds 
 would be required to see exchange .
these long residence times are 
 reflected in the high occupancies of the chelated interactions which 
 are shown in table 3 .
neither the aqvist nor merz models mg ions 
 lose an inner shell water molecule in favor of directly chelating 
 an rna atom , though this may be an issue of simulation length . a complete 
 list of chelated mg  rna interactions , separated 
 by ion model , simulation , and atom ,
( a ) 
 ion chelation ( blue ) to g8 phosphate after slv rna loop folding does 
 not disrupt the folded loop structure or association of other ions 
 ( silver ) .
( b ) ion chelation ( blue ) to the same g8 phosphate oxygen 
 disrupts the phosphate backbone , and loop folding does not occur .
in 
 our previous work we looked at the dynamics of slv rna , which 
 has been experimentally shown to adopt a mg dependent 
 u - turn structure .
our re - refinement and solution state simulations 
 of the mgfree loop showed that , in the presence of monovalent ions , 
 the slv rna adopts a fold that is similar to the mgbound loop .
when 
 mg is present , the loop becomes more compact , stabilizing 
 slv as a scaffold for substrate binding and neutralizing the charges 
 of the u - turn region to promote binding . as a proof of concept to 
 show the structure is dependent on the ion environment , we ran short 
 simulations ( hundreds of nanoseconds in length ) to show that the mgbound 
 structure is only adopted in divalent salt . in this work , we 
 extended the simulations of the mgfree loop in 
 the presence of mg ions and compare the effectiveness 
 of the aqvist , merz , and villa 12 - 6 potential mg ion models , the li 12 - 6 - 4 potential mg ion 
 model , and the panteva et al .
modifications 
 to the 12 - 6 - 4 potential model in reproducing 
 the conformational change of slv rna from mgfree to mgbound .
all five 
 models induce the conformational change , to some extent , shown by 
 low rmsds to the mgbound nmr structure in figure 2 , as well as time dependent rmsd from the 
 mgbound nmr structure in figure s1 .
though 
 we can not converge the rna structure distribution on this time scale , 
 we can look at the effect of the ions on rna structure over the ensemble 
 of structures generated in the same amount of simulation time , as 
 well as structures which have transitioned to the  folded  
 or mgbound conformation .
for the folded slv rna , this yielded at least 
 175 k individual frames per ion model , and in some cases 500 k frames , 
 i.e. , a non - negligible number of structures .
when the rna adopts 
 the same conformation , specifically a loop 
 rmsd to the mgbound reference lower than a 1.625  cutoff , the 
 mg ions associate with mgbound rna in almost identical 
 sites in all models .
the most similar are the aqvist and merz models , 
 which show mg ion association in three out of the four 
 previously described binding sites .
the subtle difference between 
 these two models can be seen in the fraction occupancies reported 
 in table 2 , with aqvist 
 ions trending toward higher fractions of occupied interactions .
the 
 villa mg ion model follows the same trend in occupancies 
 for the folded slv rna , with the exception that direct chelation occurs 
 between mg and the 3 phosphate group of residue 
 12 , leading to seemingly artificial high occupancy for this region . 

both the 12 - 6 - 4 and m12 - 6 - 4 models , which use an additional c4 potential 
 term to approximate the effect of polarizability , have higher occupancies 
 in general when compared to the 12 - 6 lj potential models .
the overestimation 
 in the binding free energy of mg to rna phosphate nonbridging 
 oxygen atoms was noted previously , and 
 corrected in the m12 - 6 - 4 model .
though 
 the occupancies are lowered in the m12 - 6 - 4 simulations compared to 
 the 12 - 6 - 4 model , they remain higher than the 126 lj potential 
 models .
overall , when slv rna is in its folded conformation , binding 
 sites 1 , 3 , and 4 show consistent occupancy in all ion models , while 
 binding site 2 shows systematically lower , but nonzero , occupancy . differences between models
arise when all data across the 10 simulations 
 performed for each model are considered .
figure 4 shows the fraction occupancies for mg associating with phosphate , base , or sugar atoms .
the models 
 trend toward the 12 - 6 - 4 having the highest occupancies , followed by 
 m12 - 6 - 4 , aqvist , mer , z and villa .
the exception to this trend can 
 be seen in specific interactions where villa , 12 - 6 - 4 , and m12 - 6 - 4 
 ions have directly chelated to the rna .
this happens at three atoms 
 with the villa parameters , but at many more rna atoms in the 12 - 6 - 4 
 potential simulations  in fact , when the entire set of 12 - 6 - 4 
 simulations is considered , there is not a single phosphate group where 
 chelation does not occur .
if 
 an ion chelates after the loop is folded , it does not necessarily 
 alter the conformation ( figure 5a ) .
however , ion chelation can prohibit folding to the correct 
 structure by trapping phosphates in alternate backbone conformations 
 precluding ion association in the experimentally determined sites 
 ( figure 5b ) .
it could 
 be bad luck that we see direct chelation in the villa simulations , 
 and if all 12 - 6 lj potential simulations were extended , we would see 
 it occur in other ion models , too .
however , on this time scale , neither 
 the aqvist nor merz models mg lose an inner shell water 
 molecule .
the chelation effect in the 12 - 6 - 4 lj potential simulations 
 is comprehensive .
the simulations performed in this study aggregate 
 15 s for 
 each ion model and converge the mg ion distribution when 
 there is only association of mg ions with slv rna binding 
 sites ( figure s3 ) .
the convergence of the mg ion distribution 
 when those ions chelate directly to rna would take simulations on 
 the order of 10s of milliseconds to observe multiple binding and unbinding 
 events for mg ions to overcome 1213.5 kcal / mol 
 free energy barriers to desolvation , and are fairly 
 prohibitive .
this is reflected in the disagreement of occupancies 
 for binding sites when all simulation data are considered ( table s5 ) .
the structures where mg is not associated 
 with rna vary between models and are not dependent on the ion models 
 themselves ( figure s4 ) . where mg is chelated to rna atoms , some simulations become trapped in structures 
 which are artificially stabilized by the chelating ion .
this means 
 that the observation that different numbers of simulations sample 
 the mgbound state is not necessarily related to the ion model used 
 in those simulations .
since this work has shown the time dependent 
 folding of slv rna , an enhanced sampling technique might be useful 
 to converge the slv rna structure distribution when mg ions are present .
the objective of this 
 study was to reproduce the mg ion dependent conformational 
 change for vsr slv rna , from the mgfree 
 to mgbound loop .
we tested five mg ion parameters which 
 used two different functional forms of the lj potential .
the 12 - 6 
 lj potential was tested with the widely used aqvist parameters , and 
 the more recent merz and villa modifications .
the 12 - 6 - 4 potential 
 li parameters , and their modifications to balance interactions between 
 mg and rna , were also tested . with all five models , we 
 observed switching from the mgfree rna to the mgbound rna structure . 

the mg ions consistently associated with common sites 
 in the folded slv rna and agreed with experimental data .
the extent 
 of folding differed between ion models and is likely a time scale 
 issue .
however , chelated interactions interfered with folding in the 
 12 - 6 - 4 lj potential simulations with both the 12 - 6 - 4 and m12 - 6 - 4 models , 
 and to a lesser extent in the 12 - 6 lj potential villa simulations . 
 in the combined data set , where mgfree and mgbound structures were 
 present , the 12 - 6 - 4 ions had the highest fraction occupancies overall . 
 of the 126 potential models , the aqvist ions had higher fraction 
 occupancies , followed by merz and then villa , following the trends 
 in increasing van der waals radii and reduced well depth .
 we compare the performance of five 
 magnesium ( mg2 + ) 
 ion models in simulations of an rna stem loop which has an experimentally 
 determined divalent ion dependent conformational shift . 
 we show that 
 despite their differences in parametrization and resulting van der 
 waals terms , including differences in the functional form of the nonbonded 
 potential , when the rna adopts its folded conformation , all models 
 behave similarly across ten independent microsecond length simulations 
 with each ion model . 
 however , when the entire structure ensemble is 
 accounted for , chelation of mg2 + to rna is seen in three 
 of the five models , most egregiously and likely artifactual in simulations 
 using a 12 - 6 - 4 model for the lennard - jones potential . despite the 
 simple nature of the fixed point - charge and van der waals sphere models 
 employed , and with the exception of the likely 
 oversampled directed 
 chelation of the 12 - 6 - 4 potential models , rna  
 mg2 + interactions via first shell water molecules are surprisingly well 
 described by modern parameters , allowing us to observe the spontaneous 
 conformational shift from mg2 + free rna to mg2 + associated rna structure in unrestrained molecular dynamics simulations . 
Introduction
Material and Methods
Results
Discussion
Conclusions