Molecular Dynamics Simulations of RNA Systems - I.B.M.C

560

34

Molecular Dynamics Simulations of RNA Systems Pascal Auffinger and Andrea C. Vaiana 34.1

Introduction

As a result of important methodological advances, the number of molecular dynamics (MD) simulations related to RNA systems has become significant, and simulations on more than 30 different RNA, RNA–protein and RNA–ligand systems have been reported (Table 34.1). However, the simulation of the dynamics of such systems presents specific problems associated with the complex three-dimensional folds adopted by their highly charged polyanionic backbone and by the large variety of naturally modified nucleotides they comprise [1]. It is, thus, timely to present the current state of the art of MD simulations of RNA systems in a comprehensive and practically oriented manner. Several reviews on nucleic acid MD simulations [2–10] already address important methodological issues. Here, we will focus on how to setup a MD simulation of an RNA system by using explicit representation of solvent and Ewald summation methods for the treatment of the long-range electrostatic interactions. Specific features related to MD simulations of RNA systems using implicit solvation models can be found in [5]. The simulation protocols that will be described in the following refer mainly to our own experience with the AMBER program with which we are most familiar and which is used in many laboratories working on nucleic acids. They should be easily transposable to other MD packages.

34.2

MD Methods

Common force field-based simulation methods make use of an empirical potential energy function and parameter set to describe the physical properties of the system to be studied. In MD simulations, the ‘‘possible’’ time evolution of the system starting from an initial set of coordinates and velocities is calculated by integrating Newton’s equations of motion over a certain period of time [11]. Descriptions of the integration algorithms associated with such programs can be found in [12– Handbook of RNA Biochemistry. Edited by R. K. Hartmann, A. Bindereif, A. Scho¨n, E. Westhof Copyright 8 2005 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim ISBN: 3-527-30826-1

(V7(M) 9/9/04 15:12) VCH/G J-1158 Hartmann PMU: WSL(W) 07/09/2004 ScalaLF&SansLF (0).3.04.05 (170240mm) pp. 560–576 III-34_P (p. 560)

34.2 MD Methods

561

Tab. 34.1. List of recent MD simulations of RNA systems (up to December 2003) using an explicit representation of the solvent and Ewald summation methods for the treatment of the long-range electrostatic interactions (simulations using truncation [107–110] or density functional methods [39] are not listed in the table).

Starting structures

nt

Model PNA.RNA duplex HNA.RNA duplex Model helix rf(ApU)12 g2 Model helix rf(CpG)12 g2

12 16 48 48

Model helix rf(CmpGm)12 g2 r(CCAACGUUGG)2 r(CCAACGUUGG)2 r(CGCGCG)2 r(GAGUACUC)2 r(GCGAGUACUCGC)2 r(CGCGAUCGCG)2 r(CCUUUCGAAAGG)2 RNRA tetraloop hairpins

Ions

Modified nt

References

2.5 1.1 2.4 2.4

Naþ Naþ @0.2 M KCl @0.2 M KCl

PNA HNA

48 20 20 12 16

4.4 2.0 1.3 0.7 5.0

@0.2 M KCl Naþ Naþ Naþ @0.3 and @1.0 M NaCl

Cm; Gm

111 112 20 19, 20, 22, 85 22 21, 83 113 23 114

24 20 24 26 10 21 62

3.0 1.4 5.0 3.0

@0.1 M NaCl Naþ NaCl Kþ

8 12 24 24

2.0 0.7 4.0 5.0

tRNAAsp

76

0.5

Naþ Naþ @0.1 M NaCl @0.3 and @1.0 M NaCl NH4 þ

tRNAAsp anticodon hairpin

17

0.5

NH4 þ

5S rRNA loop E

24 24

10.0 11.5

HIV kissing loop complexes 16S rRNA Pseudoknot Hammerhead ribozyme

46 81 26 41 41

7.5 5.5 5.0 1.1 0.8

21

10.0

U1A hairpin Human U4 snRNA X-ray ApU and GpC steps (in crystal) r(CmGmCmGmCmGm)2 r(GGACUUCGGUCC)2 r(UAAGGAGGUGAU)2

RNA–ligand (X-ray) 16S rRNA site A/neomycin B

Length (ns)

Naþ & Mg 2þ @0.2 and @1.0 M KCl and Mg 2þ Naþ and Mg 2þ @0.1 M NaCl Naþ Naþ and Mg 2þ @0.1 M NaCl and Mg 2þ

Naþ

MOE

I

Cm; Gm

D; C; m 1 G; m 5 C; m 5 U C; m 1 G

Cþ

26 115 116 117

44 23 118 114 90, 119 38, 89, 91–94, 119 50 38, 48, 49

52 120 34 121, 122 51, 123

124


562

34 Molecular Dynamics Simulations of RNA Systems Tab. 34.1 (continued)

Starting structures

RNA–protein (X-ray) U1A RNA–protein complex

nt

Length (ns)

Ions

Modified nt

References

125 126

5.0 2.2

Naþ @0.1 and @1.0 M NaCl Naþ and Cl Kþ

12 12 22 22

0.2 2.0 2.5 2.0

Naþ Naþ NH4 þ Naþ

128 129 45 24, 25

23 29 28 29

10.0 1.6 16.0 3.0

Naþ Naþ Naþ Naþ

RNA/ligand (NMR) FMN aptamer BIV tat–TAR complex

35 28

1.7 1.2

Naþ Naþ

134 135

RNA/protein (NMR) U1A RNA–protein complex dsRBD/dsRNA

30 30

1.0 2.0

Naþ Naþ and Cl

125 136

U2 snRNA–protein complex NMR GCAA tetraloop hairpin UUCG tetraloop hairpin tRNAAla acceptor stem hairpin

HIV-1(lai) SL1 hairpin HIV-1 TAR RNA hairpin H3 kissing loop complex TAR–aptamer complex

21 21

1.0 1.8

21 23

116 127

I; 2AA; IsoC; dU; Z; M; 7DAA 2AP

130, 131 132 52 133

15]. These are common to all current simulation packages and will not be further discussed here. Various methods having the goal of pushing back currently accessible timescale limits by using multiple time-step integrators are currently under development [12, 13, 16–18].

34.3

Simulation Setups 34.3.1

Choosing the Starting Structure

Starting structures are mainly derived from (1) model built systems, (2) low- to high-resolution X-ray data and (3) NMR data.


34.3 Simulation Setups

563

Model Built Structures Such structures can be used to study precise issues related to (1) the conformation and solvation of simple RNA building blocks such as Watson–Crick base pairs and base pair steps [19, 20], (2) the differences between RNA and DNA structures of similar sequence [19–21] or (3) the effects associated with the introduction of natural [22, 23] and non-natural modified nucleotides [24, 25]. In order to generate larger systems, model built fragments can also be assembled with motifs derived from X-ray or NMR structures [26]. 34.3.1.1

X-ray Structures When derived from moderate- to high-resolution data, X-ray experiments certainly provide the most accurate source of starting structures since, besides precise coordinates for the solute, they also supply important information related to its hydration and ionic environment [27, 28]. The main sources for X-ray structures are the NDB (Nucleic acid Database; http://ndbserver.rutgers.edu/ [29]) and the PDB (Protein Data Bank; http://www.rcsb.org/pdb/ [30]). 34.3.1.2

NMR Structures NMR structures are similarly deposited in the NDB and PDB. Due to the lower amount of detail contained in these structures compared to those derived from Xray data, they are generally chosen only when no X-ray data for a given system are available. Furthermore, NMR cannot provide data related to the hydration and ionic environment of the solute. 34.3.1.3

34.3.2

Checking the Starting Structure

The imprecision associated with model built and NMR structures is usually well characterized. Whereas X-ray structures are generally very precise, these may nonetheless contain local errors which, if not noticed and corrected, may seriously impair the quality of the calculated trajectories which greatly depends on the correctness of the starting structure. Conformational Checks For proteins, the precise orientation of several residues (essentially His, Asn and Glu) cannot be correctly assessed by a single examination of the electron density maps obtained at resolutions above around 1.0 A˚. Some groups have addressed this problem [31] and have proposed various automated solutions based on an examination of ‘‘contact maps’’ (see, e.g. MolProbity; http://kinemage.biochem.duke. edu/ [32]) for detecting incorrect or suspicious orientations of amino acid residues in existing experimental structures. Checking the orientation of such residues is therefore essential if one desires to perform accurate simulations of protein–RNA complexes. In nucleic acids, the orientation of the natural nucleotides can generally be determined without ambiguity even at low resolution. However, this is not 34.3.2.1


564

34 Molecular Dynamics Simulations of RNA Systems

the case for the orientation of the pseudouridine or C residue which is a naturally modified nucleotide [33] and possibly other side-chains in modified nucleotides. Protonation Issues It is now recognized that some nucleic acid residues such as adenines or cytosines can be found in their protonated state [34]. Yet, X-ray structures rarely provide direct evidence for such protonation states. These have to be derived from chemical intuition by inspecting the structural environment of the involved residues. In cases where it is impossible to decide which protonation state is the most probable, it is advisable to perform simulations of both the charged and uncharged states [34]. Evidently, protonation states of charged amino acids and of drugs such as aminoglycosides [35, 36] have also to be clearly defined before starting a simulation. 34.3.2.2

Solvent The interpretation of solvent density maps derived from X-ray experiments is often difficult. A recent survey of the NDB has revealed that it is not unusual for water molecules to be assigned to sites with very large electron densities (above those expected for a water molecule); this may be indicative of the presence of a metal ion at those locations [37]. It has also recently been emphasized that some density spots resulting from the presence of anions (SO4 2 , Cl , etc.) in the vicinity of nucleic acids have been wrongly assigned to divalent cations (Mg 2þ or Mn 2þ ) [38]. Incorrectly assigned metal ions may represent an important source of errors in MD simulations. 34.3.2.3

34.3.3

Adding Hydrogen Atoms

Most of the structures deposited in the NDB or PDB are devoid of hydrogen atoms. Hence, before starting a simulation, it is essential to add them either by using the tools provided by the main MD packages or by using external programs such as MolProbity (http://kinemage.biochem.duke.edu/ [32]). 34.3.4

Choosing the Environment (Crystal, Liquid) and Ions

Two types of environment may be chosen: the crystal and the liquid phase. Simulations of the crystal phase are rare for nucleic acids [39–43] and only the r(CpG) and the r(ApU) dinucleotides have been simulated in this environment [39, 44] even though such simulations probably provide the most accurate confrontation between experimental and theoretical data [3]. Simulations of RNA systems in the aqueous phase are the most common. The ions used to neutralize the nucleic acid can be of various types (Naþ , Kþ , NH4 þ , Mg 2þ ) and different ionic conditions ranging from a number of ions sufficient to neutralize the charge carried by the nucleotidic backbone to a number of ions that would appropriately represent



565

a solution containing 0.1–1.0 M of added salt can be selected (Table 34.1). In this last case, co-ions (Cl , SO4 2 , etc.) must be placed in the simulation box as well. 34.3.5

Setting the Box Size and Placing the Ions Box Size For simulations of the crystal phase, the box size is determined by the size of the crystal cell. In solution, the box size is usually allowed to extend by approximately 8–12 A˚ around the solute, thus ensuring three to four hydration layers in each direction. 34.3.5.1

Monovalent Ions Monovalent ions are then placed in the simulation box by replacing the water molecules with the lowest (for cations) and highest (for anions) electrostatic potential [19, 45]. If one wants to avoid biases due to the initial placement of ions, it is also essential to set limits ensuring that no ions are closer than a certain distance, usually 5–6 A˚, to any nucleic acid atom or to another ion. In this manner, no ions are initially found in the nucleic acid grooves, in direct proximity of the phosphate groups or forming contact ion pairs. 34.3.5.2

Divalent Ions The limited length of the simulations does not allow reproduction of dehydration processes of Mg 2þ ions and formation of direct Mg 2þ O contacts with the RNA which take microseconds to occur [46, 47]. Thus, divalent ions cannot be placed in the simulation box as described above. Such systems can only be realistically simulated if the Mg 2þ ions occupy their crystallographic positions already at the beginning of the MD run [48–52]. It has to be noted that, although the Mg 2þ ions occupy their crystallographic positions, the water molecules located in the first hydration shell of the ions find their equilibrium positions during the equilibration phase or, in other words, the crystallographic position of these water molecules is not taken into account. 34.3.5.3

34.3.6

Choosing the Program and Force Field Programs Two very widely used MD simulation packages are AMBER (Assisted Model Building with Energy Refinement; http://amber.scripps.edu/ [53]) and CHARMM (Chemistry at HARward Molecular Mechanics; http://yuri.harvard.edu/). Each of these is associated with its own force field. Other more recent and non-commercial MD packages such as GROMOS (http://www.igc.ethz.ch/gromos/ [54]) and NAMD (http://www.ks.uiuc.edu/Research/namd/ [55]) are also available. Package choice is a very subjective matter and is generally dictated by the lab history. However, one 34.3.6.1


566


should be attentive to the fact that some desired features may not be present in all of these packages. Force Fields AMBER [56] and CHARMM [57, 58] are the most commonly used force fields for simulating the dynamics of nucleic acids and proteins. They co-exist with the less frequently used BMS [59] and OPLS force fields [60]. However, it is essential to be aware of their limitations. Although they are all able to reproduce quite well most of the important structural and dynamical features of biological systems, they do differ in some specific aspects. Comparisons of the AMBER and CHARMM force fields performed on DNA duplexes [61, 62] have led to the conclusion that in some instances a ‘‘force field dependent polymorphism’’ [3] is observed. This polymorphism contributes to an artifactual drift of the duplexes toward an A- or a B-form that is independent of the starting conformation. Some of these problems have been corrected in recent versions of the AMBER and CHARMM force fields [63, 64]. Hence, since these force fields are constantly evolving, it is recommended to check for the latest available updates. Nevertheless, it can be concluded that up to now, none of the widely used force fields is able to reproduce fully all of the structural features experimentally observed for nucleic acids in solution [8, 64]. New, polarizable force fields [65] such as those recently implemented in AMBER are expected to bring significant improvement in the informative and predictive power of MD simulations (although at the expense of the use of more computer resources). These force fields are still at an initial stage of development and will have first to be thoroughly tested on small systems before using them for simulating solvated macromolecules. 34.3.6.2

Parameterization of Modified Nucleotides and Ligands Most empirical force fields used in common MD packages are equipped with parameters for modeling all canonical nucleotides and amino acids, but not for modified nucleotides and ligands. For simulation of such systems, new parameter sets must be developed. Although some general sets of procedures for automated parameterization based on fitting to both experimental and quantum chemical data have been proposed, force field parameterization remains a matter for experts. A classical review on force field parameterization has been written by Dinur and Hagler [66]. Fortunately some tools and methods for generating reasonable to high-quality force field parameters are available [67, 68]. Among them it is worth mentioning the AMBER module Antechamber which greatly simplifies the generation of new parameters for the AMBER force field and the Automated Frequency Matching Methodology (AFMM) developed for the CHARMM force field but extendible to all other atom based force fields [69]. AFMM allows development of parameter sets for small to medium sized molecules using high quality quantum chemical calculations as reference data. A new set of scripts (RED) has also been made available to facilitate the derivation of partial charge sets from electrostatic potential calculations (http://www.u-picardie.fr/labo/lbpd/RED/). In addition to those for nucleic acids and amino acids, some force fields have been developed 34.3.6.3



567

for specific classes of molecules, such as carbohydrates (OPLS-AA [70]; GLYCAM http://glycam.ccrc.uga.edu/), sulfates and sulfamates [71] or polyphosphates [72]. Parameters for the ions which are compatible with the chosen force field can be easily retrieved from the literature. Water Models Among the available water models, the TIP3P (Transferable Interatomic Potential – 3 Point [73]) model is generally associated with the AMBER force field. The more computationally ‘‘expensive’’ four (TIP4P [73]) and five (TIP5P [74]) point models have been developed and tested on various systems. SPC (Single Point Charge [75]) and SPC/E (SPC/Extended [76]) models are also very popular. For DNA, it has been shown that the use of TIP3P and SPC/E water models resulted in comparable hydration patterns although the water densities associated with the water model with the highest diffusion rate (TIP3P) appear to be more blurred than those calculated with the SPC/E model that has calculated diffusion rates in better agreement with experimental data [10]. Besides, a MD simulation of a protein in a crystal environment has led to the conclusion that results obtained with the SPC/E model are in much better agreement with neutron-scattering data than those collected with the TIP3P model [77]. Thus, a particular water model can significantly alter the results from MD simulations and more studies are needed to precisely evaluate their influence. 34.3.6.4

34.3.7

Treatment of Electrostatic Interactions

The ensemble of parameters in the above referenced force fields is usually considered to be the most important factor in determining the quality of a MD trajectory. Yet, if the empirical equation that models in a condensed way the intermolecular forces at play in biomolecular systems is not accurate enough or if one of its terms is not correctly evaluated then even the most precisely developed set of parameters cannot realistically account for the dynamics of these systems. Indeed, numerous studies have revealed that severe artifacts emerge when the long-range electrostatic contributions are neglected. In other words, when only a part of the Coulomb term qi q j =4pe0 rij is estimated [78–80]. It has been shown that an accurate treatment of the long-range electrostatic interactions is of paramount importance for generating realistic trajectories of nucleic acid systems. The particle mesh Ewald (PME) summation method [81, 82] has proven to lead to very stable nucleic acid trajectories [80, 83]. The program default parameters are generally adequate for using the PME method. These parameters are such that a cubic interpolation scheme and a 105 A˚ tolerance for the direct space sum cutoff is chosen. To speed up the fast Fourier transform in the calculation of the reciprocal sum, the size of the PME charge grid is chosen to be an integer power of 2, 3 and 5, and to be slightly larger than the size of the periodic box. This leads to a grid spacing of around 1 A˚ or less.


568


34.3.8

Other Simulation Parameters Thermodynamic Ensemble MD simulations can be conducted in a microcanonical ensemble or (N,V,E) for constant number of particles, volume and energy. Yet, the isothermal–isobaric or (N,P,T) ensemble for constant number of particles, pressure and temperature is more commonly chosen. Working at constant pressure is particularly important during the equilibration stage since during the building stage of the system, the water molecule positions are not optimized at the solute/solvent interface. This can lead to the generation of holes in he simulation box when constant volume options are selected. 34.3.8.1

Temperature and Pressure In order to maintain a constant temperature and pressure, the Berendsen temperature coupling scheme and an isotropic molecule based constant pressure scaling with a time constant of 0.2–0.5 ps for both are often used [84], although other coupling schemes and time constants are found in the literature. The target temperature and pressure values are usually set to 1 atm and 298 K or 25 C, also called ‘‘room temperature’’. At this point, it is interesting to note that the choice of this target temperature is dictated by early studies which used experimental data obtained at room temperature in order to parameterize the force fields. This is close to the temperature at which many in vitro experiments are performed in a laboratory. Nonetheless, simulations using different target temperatures (e.g. 5, 25 and 37 C) have been used to reveal some differences in the hydration of RNA systems [85]. 34.3.8.2

Shake, Time Steps and Update of the Non-bonded Pair List A tolerance of 0.0005 A˚ is generally used for the SHAKE algorithm [86] that allows to ‘‘safely’’ utilize a 2-ps time step (instead of 1 ps) by artificially freezing the most rapid vibration motions (CaH, OaH, NaH, etc., elongations) which can be observed in biomolecular systems. A strong tolerance value has been found to reduce the so-called ‘‘flying ice cube’’ phenomena described below [87, 88]. Another time saving device is associated with the use of a pair list for calculating non-bonded interactions between all atoms in the system. This pair list is not calculated at each time step, which would be computationally expensive, but only once every nth step (usually every 10 steps). 34.3.8.3

The Flying Ice Cube Problem Since the non-bonded pair list is not updated every step (in order to save CPU resources) and constant pressure algorithms as well as uniform scaling of velocities by the Berendsen coupling scheme are used, some small energy drain during the simulation can occur and the center of mass velocity can slowly grow. Therefore, periodically in the simulation this center of mass velocity has to be removed otherwise the ‘‘flying ice cube’’ syndrome may appear [87, 88]. It is also advised to update the non-bonded pair list more frequently [87, 88]. 34.3.8.4



569

34.3.9

Equilibration

The main aim of equilibration procedures is to alleviate tensions created in the system during the early building stages. These tensions would, if neglected, lead to unrealistic trajectories. In the following, a typical equilibration protocol used by us will be described. First, 500 steps of steepest descent minimization are applied to the entire system with periodic boundary conditions in the (N,P,T) ensemble. This is followed by 25 ps of simulation in the same thermodynamic ensemble where only the water molecules and the hydrogen atoms of the solute are allowed to move (the heavy atoms of the solute and the ions are frozen by using the BELLY option in AMBER). During the next 25 ps, the constraints placed on the monovalent ions are released so that they can start to equilibrate around the RNA system while the BELLY option is still used to freeze the heavy atoms of the solute. Then, several rounds of 50-ps MD are performed during which positional constraints of 10, 5, 2, 1, 0.5, 0.1, 0.01 and 0.001 kcal/mol A˚ 2 are applied to the heavy atoms of the solute, yielding a partially constrained 450 ps MD trajectory. Then the production phase can take place. However, in order to allow a better sampling of the conformational space by the solute and the ions, the next 0.5–1.0 ns are generally included in the equilibration phase [19]. The length of the equilibration phase can be extended or shortened at will. However, it has to be noted that the phase during which constraints are applied ensures that the solvent is appropriately equilibrated in the vicinity of the solute, but also throughout the simulation box. The part during which no constraints are applied allows then for the solute to relax in the presence of relatively well-equilibrated solvent environment. Other equilibration protocols start at low temperature (50 K); then, the temperature is raised to the target temperature at discrete time intervals [51]. An apparent advantage of the former procedure is that the sampling of the conformational space by the solvent is conducted at the target temperature and, hence, probably more efficient. 34.3.10

Sampling How Long Should a Simulation Be? Ideally, MD simulations should be long enough to sample all the conformational transitions occurring in a particular biomolecular system. Given present and predictable computational means, this goal will certainly not be achieved in the next decade. It is, thus, impossible to address in the near future biochemical processes with very long relaxation times by using classical MD simulations with explicit solvent representations. Up to now, the longest MD simulations of RNA systems have barely overcome the 10-ns timescale (Table 34.1). Nevertheless, many processes associated with limited conformational perturbations lie within reach of present methods. Among them, those that involve the binding of the smallest ligands to nucleic acids (monovalent ions and water molecules [19, 20, 48]) or the orientation 34.3.10.1


570


of the 2 0 -OH hydroxyl groups [89] can trustfully be evaluated on the current nanosecond timescale. When to Stop a Simulation Another important and often overlooked issue is: when to stop a MD simulation or when does it start to be a waste of time to continue sampling? Indisputably, a very large number of approximations are included in MD simulations. These approximations are, among others, related to: (1) imprecision in the force field parameters, (2) incomplete evaluation of the intermolecular forces, (3) neglect of polarization and charge transfer effects and (4) an incomplete knowledge of the starting conditions. This last point by itself can easily lead to the generation of partially or totally unrealistic simulations. For example, an instability of the tertiary core structure of the yeast tRNAAsp molecule characterized by a reordering of several base triples has been observed during a 500-ps MD simulation [90]. This has first been attributed to the absence of Mg 2þ ions in the model. Later on, it was proposed that an adenine involved in the above-mentioned base triples is protonated [34]. Probably both effects led to the calculated transitions and further sampling would, in this case, not have contributed to a better understanding of the structure and the dynamics of this molecule. Hence, there is a certain risk for long MD simulations to oversample regions of the ‘‘configurational space’’ attached to the current force field and MD package that are not overlapping with the ‘‘true’’ configurational space explored by the investigated system. 34.3.10.2

Multiple MD (MMD) Simulations An alternative to long MD simulations is to generate an ensemble of several ‘‘shorter’’ trajectories by using the MMD simulation technique. This method introduces slight perturbations (such as different initial velocities) in the starting conditions and, thus, exploits the chaotic nature of MD simulations in order to generate several uncorrelated trajectories [91–97]. It has been stated that several trajectories generated from similar but slightly different starting conditions may provide, on a statistical basis, more information than a single long trajectory. A parallel view developed by other authors consists of generating and comparing several trajectories obtained by starting from different initial configurations [96]. In short, MMD methods allow a statistical evaluation of MD simulation that cannot be achieved from the analysis of single trajectories [80]. 34.3.10.3

34.4

Analysis 34.4.1

Evaluating the Quality of the Trajectories

Before trying to extract information from MD simulations, it is imperative to evaluate the quality of the generated trajectories by checking their internal consistency and confronting them to all available experimental data [98].


34.4 Analysis

571

Consistency Checks Internal consistency checks are related to the detection of unphysical behaviors such as those resulting from the use of truncation methods for the evaluation of the electrostatic interactions [78, 79] or ‘‘flying ice cube’’ syndromes. When such issues come to the foreground, it is clear that the simulation can and should no longer be used to derive biologically relevant information but that the focus of the study must shift toward the correction of the detected artifacts. By doing so, the importance of the neglected contributions is generally brought to light. For example, the stabilizing role of ‘‘hydration forces’’ could be assessed by comparing trajectories issued from simulations including or discarding long-range electrostatic interactions [3, 80, 92]. 34.4.1.1

Comparison with Experimental Data Confrontation with available experimental data (that are most frequently of a structural type) is in all cases mandatory. For instance, if important tertiary interactions (Watson–Crick base pairing, etc.) present in the starting structure break during the simulation over short timescales it is probable that the simulation protocols are not adequate and the reasons underlying such behavior must be uncovered [3, 79]. Early simulations performed in ‘‘in vacuo’’ conditions or with truncation of the electrostatic interactions where disruptions of important tertiary interactions were observed clearly illustrate this point [79, 99]. Moreover, if the structure of the investigated system diverges significantly from the initial structure, in part or entirely, during the time course of the simulation, it must be considered that this observation may result from force field inaccuracies. Only comparisons with experimental data on the same or related systems can help to resolve such issues. 34.4.1.2

Visualization Visual checks are also mandatory and can be performed by using the programs that are delivered with each main MD package of by using programs such as VMD (Visual Molecular Dynamics; http://www.ks.uiuc.edu/Research/vmd/) that can read most of the available MD outputs. In rare cases, the Open BABEL utility (http://openbabel.sourceforge.net/) is needed in order to convert data written in different formats. 34.4.1.3

34.4.2

Convergence Issues

Afterwards, the convergence of the simulation with respect of the investigated properties has to be assessed. For this, it is necessary for the simulation length to exceed the average relaxation time of the investigated property which is not necessarily the longest relaxation time of the system that can largely exceed the currently accessible nanosecond timescales [98]. For example, it is possible to estimate on a statistical basis the residence times of water molecules (between 10 and 500 ps) around a regular RNA helix from nanosecond MD simulations, although the statistics will obviously be less good for the water molecules with the longest residence times [19, 20]. For less regular RNA systems, water molecules may be trapped into


572


specific pockets for much longer times and only lower bounds for the residence times can be obtained [100]. A similar problem is associated with the estimation of the binding properties of monovalent and divalent ions. While it is feasible to estimate the exchange time of monovalent ions located in RNA grooves from nanosecond MD simulations [49], it is impossible on such timescales to estimate the same values for divalent ions [48]. 34.4.3

Conformational Parameters

Most of the analysis programs delivered with major MD packages (such as CARNAL for AMBER) can monitor nearly all of the main conformational parameters such as interatomic distances, angles and dihedral angles, and also parameters that are specific to nucleic acids such as sugar puckers or user defined variables. Other programs such as CURVES (http://www.ibpc.fr/UPR9080/Curindex.html) [101] or 3DNA (http://rutchem.rutgers.edu/%7Exiangjun/3DNA/) [102] that allow to estimate parameters specific to nucleic acids such as bending, shift, slide, rise, tilt, roll, twist, etc., can be also very useful in the analysis process. The majority of these parameters are now attached to the experimental structures deposited at the NDB. 34.4.4

Solvent Analysis

Solvent analysis is an important aspect of the evaluation of MD trajectories. It is related to the characterization of water and ion binding sites and to the estimation of their residence times. Several reviews have already addressed these issues for nucleic acid systems [100, 103–106] and we will not discuss them in further details.

34.5

Perspectives

The current simulation protocols are now able to deliver stable trajectories of RNA systems ranging from small duplexes to large RNA–protein complexes over nanosecond timescales. Important conformational transition and solvation processes can be, thus, studied in detail. The current trend is to head toward longer timescales. This would allow us to address, among others, folding and unfolding issues. Nevertheless, one should be aware that by extending the currently accessible timescales, new artifacts that remain up to now hidden will appear and necessitate the development of more refined force fields and simulation methods that will most evidently incorporate polarizability and charge transfer effects. The detection and correction of artifacts that will manifest themselves on the longer timescales must be regarded as an important duty of modelers.


References

573

Acknowledgments

A. C. V. acknowledges support from the European contract HPRN-CT2002-00190 (CARBONA).

References 1 H. Grosjean, R. Benne, Modification

2

3 4

5 6 7 8 9

10

11

12

13

and Editing of RNA, ASM Press, Washington, DC, 1998. P. Auffinger, E. Westhof, in: Encyclopedia of Computational Chemistry, P. v. R. Schleyer, N. L. Allinger, T. Clark, J. Gasteiger, P. A. Kollman, H. F. Schaefer III, P. R. Schreiner (eds), Wiley, Chichester, 1998, pp. 1629–1640. P. Auffinger, E. Westhof, Curr. Opin. Struct. Biol. 1998, 8, 227–236. P. Auffinger, E. Westhof, in: Perspective in Structural Biology, M. Vijayan, N. Yathindra, A. S. Kolaskar (eds), University Press, Hyderabad, 1999, pp. 545–555. M. Zacharias, Curr. Opin. Struct. Biol. 2000, 10, 311–317. J. Norberg, L. Nilsson, Acc. Chem. Res. 2002, 35, 465–472. E. Giudice, R. Lavery, Acc. Chem. Res. 2002, 35, 350–357. T. E. Cheatham, 3rd, M. A. Young, Biopolymers 2000, 56, 232–256. T. E. Cheatham, P. A. Kollman, Annu. Rev. Phys. Chem. 2000, 51, 435– 471. T. E. Cheatham, P. A. Kollman. Structure, Motion, Interaction and Expression of Biological Macromolecules: Proceedings of the Tenth Conversation, R. H. Sarma, M. H. Sarma (eds), Adenine Press, New York, 1998, pp. 99–116. W. F. van Gunsteren, H. J. C. Berendsen, Angew. Chem. Int. Ed. 1990, 29, 992–1023. M. P. Allen, D. J. Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford, 1987. D. C. Rapaport, The Art of Molecular Dynamics Simulation, Cambridge University Press, Cambridge, 1995.

14 J. A. McCammon, S. C. Harvey, in:

15

16

17 18 19 20 21 22 23

24

25

26 27 28 29

30

Dynamics of Proteins and Nucleic Acids, J. A. McCammon, S. C. Harvey (eds), Cambridge University Press, Cambridge, 1987. A. R. Leach, Molecular Modelling: Principles and Applications, Pearson Education, Harlow, 2001. T. Schlick, E. Barth, M. Mandziuk, Annu. Rev. Biophys. Biomol. Struct. 1997, 26, 181–222. P. F. Batcho, D. A. Case, T. Schlick, J. Chem. Phys. 2001, 115, 4003–4018. T. Schlick, Structure 2001, 9, R45– R53. P. Auffinger, E. Westhof, J. Mol. Biol. 2000, 300, 1113–1131. P. Auffinger, E. Westhof, J. Mol. Biol. 2001, 305, 1057–1072. T. E. Cheatham, P. A. Kollman, J. Am. Chem. Soc. 1997, 119, 4805–4825. P. Auffinger, E. Westhof, Angew. Chem. Int Ed. 2001, 40, 4648–4650. K. Kulinska, T. Kulinski, A. Lyubartsev, A. Laaksonen, R. W. Adamiak, Comp. Chem. 2000, 24, 451– 457. M. C. Nagan, P. Beuning, K. Musier-Forsyth, C. J. Cramer, Nucleic Acids Res. 2000, 28, 2527–2534. M. C. Nagan, S. S. Kerimo, K. Musierforsyth, C. J. Cramer, J. Am. Chem. Soc. 1999, 121, 7310–7317. C. Schneider, J. Suhnel, J. Biomol. Struct. Dyn. 2000, 18, 345–352. M. Sundaralingam, B. Pan, Biophys. Chem. 2002, 95, 273–282. M. Egli, Chem. Biol. 2002, 9, 277–286. H. M. Berman, J. Westbrook, Z. Feng, L. Iype, B. Schneider, C. Zardecki, Acta Crystallogr. D 2002, 58, 889–898. H. M. Berman, T. Battistuz, T. N. Bhat, W. F. Bluhm, P. E. Bourne,


574


31 32

33

34

35

36

37

38 39

40

41

42

43

44 45 46 47 48

K. Burkhardt, Z. Feng, G. L. Gilliland, L. Iype, S. Jain et al., Acta Crystallogr. D 2002, 58, 899–907. K. McDonald, J. M. Thornton, Prot. Eng. 1995, 8, 217–224. S. C. Lovell, I. W. Davis, W. B. Arendall, 3rd, P. I. de Bakker, J. M. Word, M. G. Prisant, J. S. Richardson, D. C. Richardson, Proteins 2003, 50, 437–450. P. Auffinger, E. Westhof, in: Modification and Editing of RNA, H. Grosjean, R. Benne (eds), American Society for Microbiology, Washington, DC, 1998, pp. 103–112. Csaszar, N. Spackova, R. Stefl, J. Sponer, N. B. Leontis, J. Mol. Biol. 2001, 313, 1073–1091. M. Kaul, C. M. Barbieri, J. E. Kerrigan, D. S. Pilch, J. Mol. Biol. 2003, 326, 1373–1387. D. S. Pilch, M. Kaul, C. M. Barbieri, J. E. Kerrigan, Biopolymers 2003, 70, 58–79. U. Das, S. Chen, M. Fuxreiter, A. A. Vaguine, J. Richelle, H. M. Berman, S. J. Wodak, Acta Crystallogr. D 2001, 57, 813–828. P. Auffinger, L. Bielecki, E. Westhof, Structure 2004, 12, 379–388. Hutter, P. Carloni, M. Parrinello, J. Am. Chem. Soc. 1996, 118, 8710– 8712. D. R. Bevan, L. Li, L. G. Pedersen, T. A. Darden, Biophys. J. 2000, 78, 668–682. H. Lee, T. A. Darden, L. G. Pedersen, J. Chem. Phys. 1995, 102, 3830–3834. D. M. York, W. Yang, H. Lee, T. Darden, L. G. Pedersen, J. Am. Chem. Soc. 1995, 117, 5001–5002. N. Korolev, A. P. Lyubartsev, A. Laaksonen, L. Nordenskiold, Biophys. J. 2002, 82, 2860–2875. H. Lee, T. Darden, L. Pedersen, Chem. Phys. Lett. 1995, 243, 229–235. Nina, T. Simonson, J. Phys. Chem. B 2002, 106, 3696–3705. H. Ohtaki, Monatsh. Chem. 2001, 132, 1237–1268. H. Ohtaki, T. Radnai, Chem. Rev. 1993, 93, 1157–1204. P. Auffinger, L. Bielecki, E.

49

50

51

52

53

54

55

56

57 58 59 60

61 62 63

64

Westhof, J. Mol. Biol. 2004, 335, 555– 571. P. Auffinger, L. Bielecki, E. Westhof, Chem. Biol. 2003, 10, 551– 561. K. Reblova, N. Spackova, R. Stefl, K. Csaszar, J. Koca, N. B. Leontis, J. Sponer, Biophys. J. 2003, 84, 3564– 3582. T. Hermann, P. Auffinger, E. Westhof, Eur. J. Biophys. 1998, 27, 153–165. K. Reblova, N. Spackova, J. E. Sponer, J. Koca, J. Sponer, Nucleic Acids Res 2003, 31, 6942–6952. D. A. Pearlman, D. A. Case, J. W. Caldwell, W. S. Ross, T. E. Cheatham, S. DeBolt, D. Ferguson, G. Seibel, P. Kollman, Comp. Phys. Commun. 1995, 91, 1–41. ¨ nenberger, W. R. P. Scott, P. H. Hu I. G. Tironi, A. E. Mark, S. R. Billeter, J. Fennen, A. E. Torda, T. ¨ger, W. F. van Huber, P. Kru Gunsteren, J. Phys. Chem. A 1999, 103, 3596–3607. Kale, R. Skeel, M. Bhandarkar, R. Brunner, A. Gursoy, N. Krawetz, J. Phillips, A. Shinozaki, K. Varadarajan, K. Schulten, J. Comput. Phys. 1999, 151, 283–312. W. D. Cornell, P. Cieplak, C. I. Bayly, I. R. Gould, K. M. Merz, D. M. Ferguson, D. C. Spellmeyer, T. Fox, J. W. Caldwel, P. A. Kollman, J. Am. Chem. Soc. 1995, 117, 5179–5197. Foloppe, A. D. MacKerell, J. Comput. Chem. 2000, 21, 88–104. A. D. MacKerell, N. Banavali, J. Comput. Chem. 2000, 21, 105–120. D. R. Langley, J. Biomol. Struct. Dyn. 1998, 16, 487–509. G. A. Kaminski, R. A. Friesner, J. Tirado-Rives, W. L. Jorgensen, J. Phys. Chem. B 2001, 105, 6474–6487. M. Feig, B. M. Pettitt, J. Chem. Phys. B 1997, 101, 7361–7363. M. Feig, B. M. Pettitt, Biophys. J. 1998, 75, 134–149. T. E. Cheatham, P. Cieplak, P. A. Kollman, J. Biomol. Struct. Dyn. 1999, 16, 845–862. S. Y. Reddy, F. Leclerc, M. Karplus, Biophys. J. 2003, 84, 1421–1449.


References 65 T. A. Halgren, W. Damm, Curr. Opin. 66

67 68

69

70

71 72

73

74 75

76

77 78 79

80

81 82 83

Struct. Biol. 2001, 11, 236–242. U. Dinur, A. T. Hagler, in: Reviews in Computational Chemistry, K. B. Lipkowitz, D. B. Boyd (eds), VCH, New York, 1991, pp. 99–164. J. Wang, P. Kollman, J. Comput. Chem. 2001, 22, 1219–1228. S. L. Njo, W. F. van Gunsteren, F. Muller-Plathe, J. Chem. Phys. 1995, 102, 6199–6207. A. C. Vaiana, A. Schulz, J. Wolfrum, M. Sauer, J. C. Smith, J. Comput. Chem. 2003, 24, 632–639. D. Kony, W. Damm, S. Stoll, W. F. Van Gunsteren, J. Comput. Chem. 2002, 23, 1416–1429. C. J. M. Huige, C. Altona, J. Comput. Chem. 1995, 16, 56–79. K. L. Meagher, L. T. Redman, H. A. Carlson, J. Comput. Chem. 2003, 24, 1016–1025. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, M. L. Klein, J. Chem. Phys. 1983, 79, 926– 936. M. W. Mahoney, W. L. Jorgensen, J. Chem. Phys. 2000, 112, 8910–8922. H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, J. Hermans, in: Intermolecular Forces, B. Pullman (ed.), Reidel, Dordrecht, 1981, pp. 331–342. H. J. C. Berendsen, J. R. Grigera, T. P. Straatsma, J. Phys. Chem. 1987, 97, 6269–6271. M. Tarek, D. J. Tobias, Biophys. J. 2000, 79, 3244–3257. Auffinger, D. L. Beveridge, Chem. Phys. Lett. 1995, 234, 413–415. G. Ravishanker, P. Auffinger, D. R. Langley, B. Jayaram, M. A. Young, D. L. Beveridge, in: Reviews in Computational Chemistry, K. B. Lipkowitz, D. B. Boyd (eds), VCH, New York, 1997, pp. 317–372. S. Louise-May, P. Auffinger, E. Westhof, Curr. Opin. Struct. Biol. 1996, 6, 289–298. C. Sagui, T. A. Darden, Annu. Rev. Biophys. Struct. 1999, 28, 155–179. T. Darden, L. Perera, P. Li, L. Pedersen, Structure 1999, 7, R55–R60. T. E. Cheatham, J. L. Miller, T. Fox,

84

85 86

87

88

89 90 91

92

93 94 95

96 97 98 99

100 101 102 103

575

T. A. Darden, P. A. Kollman, J. Am. Chem. Soc. 1995, 117, 4193–4194. H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, J. Chem. Phys. 1984, 81, 3684–3690. Auffinger, E. Westhof, Biophys. Chem. 2002, 95, 203–210. J. P. Ryckaert, G. Ciccotti, H. J. C. Berendsen, J. Comput. Phys. 1977, 23, 327–336. S. C. Harvey, R. K. Z. Tan, T. E. Cheatham, J. Comput. Chem. 1998, 19, 726–740. S. W. Chiu, M. Clark, S. Subramaniam, E. Jakobsson, J. Computat. Chem. 2000, 21, 121–131. Auffinger, E. Westhof, J. Mol. Biol. 1997, 274, 54–63. P. Auffinger, S. Louise-May, E. Westhof, Biophys. J. 1999, 76, 50–64. P. Auffinger, S. Louise-May, E. Westhof, J. Am. Chem. Soc. 1995, 117, 6720–6726. P. Auffinger, S. Louise-May, E. Westhof, J. Am. Chem. Soc. 1996, 118, 1181–1189. P. Auffinger, E. Westhof, Biophys. J. 1996, 71, 940–954. P. Auffinger, E. Westhof, J. Mol. Biol. 1997, 269, 326–341. M. Braxenthaler, R. Unger, D. Auerbach, J. A. Given, J. Moult, Proteins 1997, 29, 417–425. L. S. D. Caves, J. D. Evanseck, M. Karplus, Prot. Sci. 1998, 7, 649–666. H.-B. Zhou, L. Wang, J. Phys. Chem. 1996, 100, 8101–8105. W. F. van Gunsteren, A. E. Mark, J. Chem. Phys. 1998, 108, 6109–6116. Swaminathan, G. Ravishanker, D. L. Beveridge, J. Am. Chem. Soc. 1991, 111, 5027–5040. P. Auffinger, E. Westhof, Biopolymers 2000, 56, 266–274. Lavery, H. Sklenar, J. Biomol. Struct. Dyn. 1989, 6, 655–667. X. J. Lu, W. K. Olson, Nucleic Acids Res. 2003, 31, 5108–5121. P. Auffinger, E. Westhof, in: Water Management in the Design and Distribution of Quality Food, Y. H. Ross, R. B. Leslie, P. J. Lillford (eds), Technomic, Basel, 1999, pp. 165–198.


576

34 Molecular Dynamics Simulations of RNA Systems 104 M. Feig, B. M. Pettitt, Biopolymers 105

106

107 108 109

110 111

112

113

114 115

116

117

118

119

1998, 48, 199–209. V. A. Makarov, B. M. Pettitt, M. Feig, Acc. Chem. Res. 2002, 35, 376– 384. P. Auffinger, B. Masquida, E. Westhof, in: Computational Methods for Macromolecules: Challenges and Applications, T. Schlick, H. H. Gan (eds), Springer Verlag, Heidelberg, 2002, pp. 61–70. Y. Tang, L. Nilsson, Biophys. J. 1999, 77, 1284–1305. A. Lahiri, L. Nilsson, Biophys. J. 2000, 79, 2276–2289. J. Sarzynska, T. Kulinski, L. Nilsson, Biophys. J. 2000, 79, 1213– 1227. J. Norberg, L. Nilsson, J. Chem. Phys. 1996, 104, 6052–6057. R. Soliva, E. Sherer, F. J. Luque, C. A. Laughton, M. Orozco, J. Am. Chem. Soc. 2000, 122, 5997–6008. H. De Winter, E. Lescrinier, A. Van Aerschot, P. Herdewijn, J. Am. Chem. Soc. 1998, 120, 5381–5394. K. E. Lind, V. Mohan, M. Manoharan, D. M. Ferguson, Nucleic Acids Res. 1998, 26, 3694– 3799. Y. Pan, A. D. MacKerell, Nucleic Acids Res. 2003, 31, 7131–7140. J. Sarzynska, L. Nilsson, T. Kulinski, Biophys J 2003, 85, 3445– 3459. F. Pitici, D. L. Beveridge, A. M. Baranger, Biopolymers 2002, 65, 424– 435. J. Guo, I. Daizadeh, W. H. Gmeiner, J. Biomol. Struct. Dyn. 2000, 18, 335– 344. C. Schneider, M. Brandl, J. ¨ hnel, J. Mol. Biol. 2001, 305, 659– Su 667. P. Auffinger, S. Louise-May, E. Westhof, Faraday Discuss. 1996, 103, 151–174.

120 W. Li, B. Ma, B. A. Shapiro, Nucleic

Acids Res. 2003, 31, 629–638. 121 R. A. Torres, T. C. Bruice, Proc. Natl

122 123

124

125 126 127 128 129 130

131 132

133

134 135

136

Acad. Sci. USA 1998, 95, 11077– 11082. R. A. Torres, T. C. Bruice, J. Am Chem. Soc. 2000, 122, 781–791. Hermann, P. Auffinger, W. G. Scott, E. Westhof, Nucleic Acids Res. 1997, 25, 3421–3427. J. L. Asensio, A. Hidalgo, I. Cuesta, C. Gonzalez, J. Canada, C. Vicent, J. L. Chiara, G. Cuevas, J. JimenezBarbero, Chemistry 2002, 8, 5228– 5240. C. M. Reyes, P. A. Kollman, RNA 1999, 5, 235–244. Hermann, E. Westhof, Nat. Struct. Biol. 1999, 6, 540–544. J. Guo, W. H. Gmeiner, Biophys. J. 2001, 81, 630–642. D. A. Zichi, J. Am. Chem. Soc. 1995, 117, 2957–2969. J. Miller, P. A. Kollman, J. Mol. Biol. 1997, 270, 436–450. F. Kieken, E. Arnoult, F. Barbault, F. Paquet, T. Huynh-Dinh, J. Paoletti, D. Genest, G. Lancelot, Eur. Biophys. J. 2002, 31, 521–531. G. La Penna, D. Genest, A. Perico, Biopolymers 2003, 69, 1–14. T. Kulinski, M. Olejniczak, H. Huthoff, L. Bielecki, K. PachulskaWieczorek, A. T. Das, B. Berkhout, R. W. Adamiak, J. Biol. Chem. 2003, 278, 38892–38901. F. Beaurain, C. Di Primo, J. J. Toulme, M. Laguerre, Nucleic Acids Res. 2003, 31, 4275–4284. ¨hnel, Biopolymers C. Schneider, J. Su 1999, 50, 287–302. C. M. Reyes, R. Nifosi, A. D. Frankel, P. A. Kollman, Biophys. J. 2001, 80, 2833–2842. T. Castrignano, G. Chillemi, G. Varani, A. Desideri, Biophys. J. 2002, 83, 3542–3552.
