Molecular dynamics simulations

Protocol: Differences in Ion-RNA Binding Modes Due to Charge Density Variations Explain the Stability of RNA in Monovalent Salts

This protocol is based on the GROMACS 5.05 package.

1. Simulation Box Preparation

1.1 RNA Structure Preparation in PyMOL

1. Prepare an RNA duplex with the sequence 5’-GGCUCUGGUUAGAC using the Nucleic Acid Builder [1] or similar tools.
2. Open the corresponding PDB file of this duplex structure in PyMOL.
3. Load the PDB entry for the crystal structure of HIV-TAR using: 

     fetch 1ANR

1. Use the 'drag coordinates' option in PyMOL to align 1ANR with the duplex structure: 
   • Ensure the 5’ end of the duplex matches the 3’ end of 1ANR and vice versa.
2. Select the sequence of both structures and copy the selection to a new object.
3. Verify that the sequence of this new object matches the full hairpin: 5′-GGCUCUGGUUAGACCAGAUCUGAGCCUGGGAGCUCUCUGGCUAACUAGGGCC-3′
4. Save the final structure as a new PDB file.
5. Preparation of the Simulation Box
6. Generate GROMACS coordinate and topology files using: gmx pdb2gmx
   • Input: PDB file from Step 1.1
   • Select the ff99-parmbsc0 force field [2] and TIP3P water model [3].
   • Modify ffnonbonded.itp for ion parameters if necessary.
7. Generate a simulation box using gmx editconf with the flags -bt triclinic -box 42.0 6.5 12.5 -center 21.0 3.25 6.75
8. Fill the box with water molecules using gmx solvate with the flags -cs spc216.gro
9. Replace water molecules with 400 mM NaCl or KCl using gmx genion
10. Perform energy minimization to eliminate high-energy contacts with the following mdp file parameters: 

integrator = steep

emtol = 100.0 

emstep = 0.001

nsteps = 50000 

nstlist = 1 

ns_type = grid 

rlist = 0.8 

coulombtype = PME 

rcoulomb = 0.8 

rvdw = 0.8 

pbc = xyz 

2. Equilibration of Ions and Water

• Heavy atoms of the RNA molecule are restrained with 1000 kJ/nm² force constant.
• The posre.itp file (generated via pdb2gmx) should be linked in topol.top: 

#ifdef POSRES

#include "posre.itp"

#endif

• Apply Periodic Boundary Conditions (PBC) in three directions.
• Long-range electrostatic interactions are computed via Particle Mesh Ewald (PME) [4].
• Set temperature to 300 K using Berendsen thermostat [5] and pressure to 1 bar with Parrinello-Rahman barostat [6].
• Bond constraints: 
  • SETTLE for water molecules.
  • LINCS for RNA covalent bonds.
• Integration scheme: Leap-Frog with 1 fs time step.

2.1 2.5-ns NPT Simulation at 300 K & 1 bar

Key mdp parameters:

title             = Eq NPT run 

; Run parameters

define =-DPOSRES 

integrator        = md            

nsteps            = 2500000       

dt                = 0.001         

; Output control

nstxout           = 1000          

nstvout           = 5000          

nstxtcout = 5000          

nstenergy = 5000          

nstlog            = 5000          

; Bond parameters

continuation      = no             

;constraint_algorithm = lincs      

;constraints      = all-bonds     

;lincs_iter       = 1             

;lincs_order      = 4             

; Neighborsearching

ns_type           = grid          

nstlist           = 25            

rlist             = 1.0           

rcoulomb  = 1.0           

rvdw              = 1.0           

; Electrostatics

coulombtype       = PME           

pme_order = 4             

fourierspacing    = 0.16  

; Temperature coupling is on

tcoupl            = V-rescale

tc-grps           = System        

tau_t             = 0.5   

ref_t             = 300   

; Pressure coupling is on

pcoupl            = Parrinello-Rahman     

pcoupltype        = isotropic     

tau_p             = 2.0   

ref_p             = 1.0   

compressibility = 4.5e-5  

refcoord_scaling = com

; Periodic boundary conditions

pbc               = xyz           

; Dispersion correction

DispCorr  = EnerPres      

; Velocity generation

gen_vel           = no

2.2 60-ns NVT Simulation

• Uses the final snapshot of the NPT simulation as input.
• Same settings as NPT but barostats are turned off.

3. Generation of the Unfolding Pathway

3.1 Definition of 3’ and 5’ Termini

• Use make_ndx to define in index.ndx : 
  • Residue 1 as Chain_A
  • Residue 52 as Chain_B

3.2 300-ns Steered MD Simulation

• Apply slow pulling along the x-direction using: 

pull                    = yes

pull-ngroups            = 2

pull-ncoords            = 1

pull_group1_name        = Chain_A 

pull_group2_name        = Chain_B 

pull_coord1_groups       = 1 2

pull_coord1_type        = umbrella 

pull_coord1_geometry    = direction-periodic

pull-coord1-vec              = 1.0 0.0 0.0

pull_coord1_dim              = Y N N

pull_coord1_rate        = 0.0001    ; in nm/ps       

   pull_coord1_k           = 1000  ; kJ/(mol×nm²)         

pull_coord1_start            = yes           

pull_nstxout                 = 100 

4. Umbrella Sampling

• Extract frames every 0.2 nm from the pulling trajectory using gmx trjconv with -dump <time>
• Generate a total of 142 gro files representing different separation lengths.
• Run 100-ns umbrella sampling simulations [7] using the same mdp parameters as in steered MD, but with: pull_coord1_rate = 0.0

5. Free Energy Landscape Calculation

• Compute the Potential of Mean Force (PMF) using Weighted Histogram Analysis Method (WHAM) [8]: 
  1. Adjust box sizes to prevent PMF calculation issues by setting a unique x-value (e.g., 60.000) in all gro files, and generate new tpr files based on this gro files using gmx grompp.
  2. Filter pull force files (pullf-umbrella.xvg) to exclude values outside the range -500 to 500.
  3. Prepare input files: 

tpr-files.dat   # List of umbrella tpr files

pullf-files.dat # List of pull force files

1. Run gmx wham: 

gmx wham -if pullf-files.dat -o profile.xvg -hist histo.xvg -bsres -nBootstrap 100 -it tpr-files.dat -unit kCal

References:

[1] T. J. Macke, D. A. Case, Am. Chem. Soc. 682, 379–393 (1997)

[2] A. Pérez, I. Marchán, D. Svozil, J. Sponer, T. E. Cheatham III, C. A. Laughton, M. Orozco, Biophys. J. 92, 3817–3829 (2007)

[3] W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, M. L. Klein, J. Chem. Phys. 79, 926–935 1983.

[4] T. Darden, D. York, L. Pedersen, J. Chem. Phys. 98, 10089–10092 1993.

[5] H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, J. R. Haak, J. Chem. Phys., 81, 3684–3690 1984.

[6] M. Parrinello, A. Rahman, J. Appl. Phys. 52, 7182–7190 1981.

[7] G. M. Torrie, J. P. Valleau, Chem. Phys. Lett. 28, 578–581 (1974)

[8] S. Kumar, J. M. Rosenberg, D. Bouzida, R. H. Swendsen, P. A. Kollman, The method. J. Comput. Chem. 13, 1011–1021 (1992).