Molecular Dynamics Simulations of l-RNA Involving Homo- and

Subscriber access provided by University of Newcastle, Australia

Article

Molecular dynamics simulations of L-RNA involving homo- and heterochiral complexes Marta Dudek, and Joanna Trylska J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.6b01075 • Publication Date (Web): 08 Feb 2017 Downloaded from http://pubs.acs.org on February 10, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Molecular dynamics simulations of l-RNA involving homo- and heterochiral complexes† Marta Dudek∗,‡,¶,§ and Joanna Trylska∗,‡ Centre of New Technologies, University of Warsaw, Banacha 2c, 02-097 Warsaw, Poland, Institute of Biochemistry and Biophysics, Polish Academy of Sciences, Pawińskiego 5a, 02-106 Warsaw, Poland, and Department of Hematology, Oncology and Internal Diseases, Medical University of Warsaw, Al. Żwirki i Wigury 61, 02-091 Warsaw, Poland E-mail: [email protected]; [email protected]

Abstract l-RNA is a mirror image of the naturally occurring d-RNA. We present the first simulations of l-RNA/d-RNA complexes in the AMBER ff10 force field that we modified to account for l-ribonucleotides. To validate the modifications we performed molecular dynamics simulations of several homo- and heterochiral RNA complexes. We did not observe any canonical base pairing in the l-/d- systems but we found many repeatable geometric patterns that suggest possible structural motifs between the heterochiral RNA strands. The microscopic picture of l-RNA/d-RNA interface could help formulate the rules for designing l-RNAs with desired affinity towards specific d-RNA motifs. The modified force field can be downloaded as Supporting Information and at https://github.com/martadudek/lrna-amber †

Molecular dynamics of l-RNAs To whom correspondence should be addressed ‡ Centre of New Technologies ¶ IBB PAS § Medical University ∗

1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Introduction l-RNA differs from d-RNA only in the sugar moiety that consists of β-l-ribose instead of β-d-ribose (see Fig. 1). Naturally occurring RNA is d-RNA and no l-RNA’s are known to exist in nature.

Figure 1: Structure of l- and d-guanosine. So far, Protein Data Bank (PDB) 1 contains three-dimensional structures of only three complexes involving l-RNA sequences, also referred to as Spiegelmers. The first structure is of an RNA duplex in the l-configuration composed of l-r(CUGGGCGG)/l-r(CCGCCUGG) strands. 2 The sequence is a mirror image of a naturally abundant fragment of 5S ribosomal RNA from Thermus flavus. Two other structures are of nucleic acid aptamers bound to their specific targets. One of them is a mixed l-RNA/l-DNA aptamer in the complex with a dprotein fragment, 3 a mirror image of a natural complex. Another structure is of a unnatural l-RNA aptamer bound to a chemokine termed CCL2, a protein composed of natural l-amino acids. 4 The inverse chirality of nucleic acid oligomers makes them biostable because they are not recognized and thus degraded by nucleases. In 1996 Klussmann et al. synthesized l-RNA that specifically binds d-adenosine 5 and have shown that it resists degradation in human plasma even after 60 hours of incubation at 37◦ C. Moreover, l-nucleic acids display low capacity to bind antibodies 6 so they do not induce severe self-reaction of the immune system. The affinity of Spiegelmers towards various targets has been already shown in vivo. 7–13 The 2


Page 2 of 29

Page 3 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


sequences of Spiegelmers were not rationally designed but were found experimentally with the SELEX technique (Systematic Evolution of Ligands by Exponential Enrichment). 14 The targets of these Spiegelmers varied from amino acids to hormones, proving wide applicability of l-RNA aptamers. In clinical trials Spiegelmers appear safe and well-tolerated. 15 All these properties make l-RNAs promising molecules to be tested as drugs, even though detailed mechanisms of their activity remain either unpublished or unknown. Interactions between l- and d-RNAs have not been systematically investigated and the literature on the topic is sparse. A notable example of the fundamental research in this area is by Gaubert et al. 16 who reported the synthesis of d-oligonucleotides with few randomly placed left-handed substitutions. Later, the effect of such substitutions on the affinity towards a partially complementary (in the Watson-Crick sense) d-RNA sequence was assessed. 17 Related studies typically investigate short sequences, up to 10 nucleotides, which are only partially substituted, e.g., by three l-thymidine monomers. 18 Such mixed duplexes have similar melting temperatures as pure d-DNA/d-RNA duplexes suggesting that introducing few enantiomeric nucleotides in the sequence does not disrupt the helical structures of nucleic acid hybrids. Unfortunately, the experimental data on the interactions of all-l/all-doligonucleotides are either fragmentary or inconsistent. Garbesi et al. reported that l-DNA of a random sequence and complementary (in the Watson-Crick sense) d-DNA do not bind. 19 However, other experiments have shown interactions between the d- and l-homooligomers among both deoxy- and ribonucleotides. 20 Surprisingly, most studies focused on the interactions of complementary l- and d strands 21 leaving the possible non-canonical heterochiral base pairs unexamined. As described above the research on l-oligonucleotides has been conducted by experimental techniques and to the best of our knowledge there are no reports on computational modeling of l-RNA. We present the first approach to model heterochiral RNA complexes and investigate their structural dynamics. We applied all-atom molecular dynamics 22 (MD) in explicit solvent. To perform MD simulations we had to modify and next verify the AM-

3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 29

BER ff10 23 force field parameters. Along with test runs we simulated the association of two equally long, separated strands of mixed chirality. From MD trajectories we identified possible hydrogen bond patterns between the l- and d-RNA sequences. All MD simulations totalled to 3.2 µs.

Theory and Computational Details Force fields: variant and invariant terms AMBER 24 is a well-established force field for MD simulations. The potential energy function used in this force field is given by the following equation:

EAMBER =

X

kl (rl − rl0 )2 +

l∈bonds

+

X m∈angles

X (q,r)∈atom pairs

X

X

p∈dihedrals

n

km (φm − φ0m )2 +

Aqr Bqr + 12 + 6 rqr rqr

X (s,t)∈atom pairs

Vpn [1 + cos(nφp − γpn )] +

1 qs q t . 4πε0 rst

The first three terms describe atom-type dependent bonded potential energy terms with their respective force constants and equilibrium values. The parameters of the third dihedral term, describing torsional rotation around a bond, are the dihedral phase γpn and periodicity n. The summations in the last two terms are over the nonbonded atom pairs separated by at least three bonds: the Lennard-Jones potential and Coulombic electrostatic term. Since d- and l-nucleotides are mirror images we focus on the symmetry of Eq. (1) upon spatial inversion, i.e., upon the transform which changes the signs of all coordinates: (x, y, z) 7→ (−x, −y, −z). Note, the inversion in three-dimensional space is equivalent to the change of sign of just one coordinate, e.g., (x, y, z) 7→ (−x, y, z), which can be interpreted as the mirror reflection of the system. For RNA systems this reflection means the transition from d-RNA to l-RNA.

4


(1)

Page 5 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The energy associated with bond lengths is invariant upon inversions (since inversion is an isometry). The same applies to Lennard-Jones and electrostatic terms. Also, the energy of planar angles remains unchanged (since inversion is a conformal map). The only energy term that may change upon inversion is the dihedral angle term because each torsion angle φp changes its sign when the system is reflected (φp 7→ −φp ). The AMBER force field parameters solely responsible for the symmetry of dihedrals (or lack of thereof) are the γpn phases. If γpn is equal to 0, the maximum of the cosine term is at φp = 0. Conversely, if γpn is equal to π, the energy minimum is at φp = 0. If all γpn are equal either 0 or π, the dihedral energy terms are symmetric under the inversion (since cosine is an even function). In the older versions of the AMBER force field ( 25 up to 26 ), the γpn phases for nucleic acids are always equal to 0 or π. As a result, the dihedral angle energy term is an even function of φp and does not change under the inversion. Therefore, the older AMBER force fields can be readily used to simulate enantiomers of nucleic acids. However, the above does not hold for the more recent 26–28 modifications of the ff99 force field, such as parmbsc0 and χOL3, which are otherwise necessary to correctly describe the dynamics of RNA on longer time scales. 24 Therefore, in silico studies of l-RNA can be performed either in the earlier (pre-2007) AMBER force field or in its newer version but with modifications. These necessary modifications are described in the next section.

l-modification of ff10 AMBER force field: reflection of parameters As already stated, the ff10 AMBER force field has many phase terms that are neither equal to 0 nor π. The consequence of this fact is the violation of parity which is illustrated for the CT-CT-CI-OS dihedral (Fig. 2). To modify AMBER ff10 we first identified eight dihedral angle terms, reported in Table 1, that required correction. The angles, for the adenine dinucleotide as an example, are shown in Figure 3. For each term we introduced its mirror reflection. Second, we defined all left-handed (l-A, l-U, l-G, l-C) ribonucleotides and named them as LA, LU, LG, LC. The new nucleotides were given new atom types 5



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Page 6 of 29

Page 7 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


for O4’ (previously OS, now LD) and O5’ (previously OH, now LH). The choice which atom types to replace was dictated by the fact that all non-symmetric dihedrals contain at least one OS or one OH atom. These two atom types LD and LH were introduced to the AMBER files, following the patterns for OS and OH, respectively. Next, the sign of γpn -s was inverted for non-symmetric dihedral terms that involved the LD and LH atoms. The modified files (leaprc.ff14SB, parm10_LRNA.dat, nucleic12_LRNA.lib, nucleic10.in) are described in Section S1 and are deposited in Supplementary Information. These files can be readily used in leap.

Figure 3: Atom labels for dihedral angles of the adenine dinucleotide α, γ and χ that were modified in our correction of the AMBER ff10 (refer to Tab. 1 for the dihedral angles).

Molecular dynamics simulations Nucleic Acid Builder 29 was used to build ideal A-form helices composed of two 12 nucleotide long strands. Six distinct starting structures were prepared: two ideal anti-parallel A-form 7



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table 1: Non-symmetric dihedral angles from the parmbsc0 and χOL3 corrections with the new atom labels introduced by us. no. dihedral angle 1) OS-P-OS-CI 2) OH-P-OS-CI 3) CT-CT-CI-OS 4) CT-CT-CI-OH 5) OS-CT-N*-C5 6) OS-CT-N*-CP 7) OS-CT-N*-C4 8) OS-CT-N*-CS

conventional name α α γ γ χ χ χ χ

ff amendment bsc0 bsc0 bsc0 bsc0 chiOL chiOL chiOL chiOL

atom labels for l-RNA LD-P-LD-CI LH-P-LD-CI CT-CT-CI-LD CT-CT-CI-LH LD-CT-N*-C5 LD-CT-N*-CP LD-CT-N*-C4 LD-CT-N*-CS

helices (d-A12 /d-U12 , l-A12 /l-U12 ), two pairs of anti-parallel separated strands (d-A12 ..dU12 , l-A12 ..l-U12 ) and two pairs of parallel separated strands (l-A12 ..d-U12 , l-U12 ..d-A12 ). The separation was achieved in Chimera 30 by moving the centroid of one strand away from the other by 23 Å. The reflection of strands was performed in a text editor by changing for each atom the sign of one coordinate. Note that we expected similar results for some starting structures: (d-A12 /d-U12 and l-A12 /l-U12 ) and (d-A12 ..d-U12 and l-A12 ..l-U12 ). Such redundant computations helped to assess the correctness of the modified force field. Next, the hydrogen atoms and counterions (to neutralize the negative charge of the system) were added. Then with leap each structure was immersed in a rectangular box of TIP3P water molecules, ensuring that the minimum distance from the water box edges to the solute was at least 12 Å. The starting coordinates for the l-..l- (homochiral) systems were obtained by reflecting the entire d-..d- setup after adding the water molecules. This was because leap adds water molecules randomly, so hydrating d-..d- and l-..l- independently would break the symmetry between systems, which we initially did not want. Also, inverting the solute with fixed water positions would result in clashes. For heterochiral systems that were simulated from separated strands, the initial coordinates of the l-..d- structure were obtained by reflection of the corresponding d-..l- structure prior to solvating the systems. The purpose of this approach was to introduce extra diversity into simulations of heterochiral complexes which should allow for better sampling of their structural patterns. However,

8


Page 8 of 29

Page 9 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


while simulating RNA enantiomers, one can either invert the entire system (including water molecules) or solvate strands after performing the inversion. Ionic strength of ∼300 mM NaCl was achieved by further addition of sodium cations and chloride anions. High ionic strength improves the stability of heterochiral complexes, as indicated by our UV-monitored thermal melting experiments (data not shown). Energy minimizations and MD simulations were performed with the AMBER 23 and NAMD 31 packages, respectively. First, energy minimizations were carried out with sander; 8500 steps with steepest descents followed by 1500 steps with the conjugate gradient methods. For homochiral left- and right-handed systems the initial single point energies were identical. Electrostatic interactions were treated by Particle Mesh Ewald (PME) summation, 32 in periodic boundary conditions. The cutoff distance for evaluation of short-range non-bonded interactions was set to 12 Å. The temperature and pressure were kept constant with the Langevin thermostat and Berendsen piston, 33 respectively. Bonds involving hydrogen atoms were restrained using the SHAKE algorithm, 34 and the simulation time step was set to 2 fs. In the subsequent 0.5 ns-long thermalization stage, the temperature was gradually raised from 30 to 310 K and the constraints on heavy atoms of RNA were weakened to k = 50 kcal/(mol·Å2 ) for the first 85 ps and k = 25 kcal/(mol·Å2 ) for the next 35 ps. Following our earlier protocol 35 the equilibration stage was divided into two parts. First, the constraints were gradually released in 6 rounds of 50 ps each. Next, in the 600 ps long equilibration stage the restraints were relaxed, keeping only four terminal phosphate atoms restrained with small harmonic potential of 0.35 kcal/(mol·Å2 ) in order to preserve the structural integrity of the duplex strands. The homochiral helices were simulated for 200 ns in three independent production runs and for the initially separated systems each production simulation lasted 500 ns. ptraj was used to cluster the trajectory frames into 10 groups using the hierarchical agglomerative average-linkage clustering algorithm with RMS as a metric. The number of hydrogen bonds between base pairs in the RNA helix was determined with 3DNA. 36 Stacking interactions,

9



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

occupancies of inter- and intra-strand hydrogen bonds and base pair types and configurations were analyzed with MINT 37 and visualized with VARNA. 38 The following criteria were applied for detection of hydrogen bonds (including the water-mediated ones): the maximal distance between the donor and acceptor of 3.5 Å and the angle between the donor-hydrogen and acceptor not lower than 150 degrees. VMD 39 was used to calculate ion and water densities around tertiary structures. Plots were prepared in R and Grace.

Results and Discussion of Molecular Dynamics Simulations Validation of the force field modification for homochiral duplexes The force field modification was first tested and evaluated in MD simulations of homochiral helices (d-A12 /d-U12 and its mirror-image l-A12 /l-U12 ; both starting from a perfect A-form helix, see Methods). Note that the dynamics of these systems will not be equally symmetric, even though the reflections of systems’ coordinates are exact, due to numerical errors in MD simulations arising from the integration algorithm and different random starting velocities. However, the global dynamical picture should be similar. The representative structures from clustering of trajectories are presented in Fig. 4. For d-A12 /d-U12 the most common clusters differ only in the position of the terminal U13, which was flipped out in 20% of frames. The average Root Mean Square Deviation (RMSD) confirms the overall stable dynamics of both homochiral duplexes (Fig. S1). RMS fluctuations (RMSF) are about 1 Å for the majority of nucleotides in the helix and, as expected, only for the terminal nucleotides increase up to 4 Å (Fig. 5). Helical parameters of these homochiral duplexes display similar absolute values and are shown in Tab. S4. Other selected parameters such as Twist, Opening and Rise are presented in Figures S2, S3, S4. The hydrogen bond types and their occurrence in homochiral duplexes is listed in Tab. 2

10


Page 10 of 29

Page 11 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 4: a) l-A12 /l-U12 and b) d-A12 /d-U12 representative structures obtained from trajectory clustering (numbers denote the percentage of cluster occupancy).

and shows that all base pairs (except the terminal ones) were present during more than 90% of the simulation time. This is further confirmed by a similar number of hydrogen bonds per each base pair shown in Fig. S5 for the l-A12 /l-U12 and d-A12 /d-U12 systems. In Fig. 6 we present the secondary structures of homochiral duplexes with nucleotides colored according to the number of hydrogen bonds they form via the WC and non-WC edges, as well as their van der Waals energies. Overall, the simulations for the homochiral systems starting from the A-form helix indicate that our modification to the AMBER force field correctly handles the l-ribonucleotides.

Simulations of initially separated strands As the next step we simulated the following homo- and heterochiral systems: l-A12 ..dU12 , l-U12 ..d-A12 , d-A12 ..d-U12 , and l-A12 ..l-U12 starting from the RNA strands initially separated (see Methods). The choice was limited to these two nucleotide types but our force field modification includes also l-guanosine and l-cytidine monophosphates, and can 11



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 5: a) Heavy-atom RMSF averaged per nucleotide, computed for six MD production simulations (3 for d-A12 /d-U12 and 3 for l-A12 /l-U12 ). b) Histograms of heavy-atom RMSF comparing the occupancies of RMSF values in each simulation (colors are as in the legend of a)).

be easily extended to modified nucleotides. Fig. 7 presents the evolution of RMSD with respect to the last frame of each simulation, showing that the structures equilibrated at particular conformations. For the homochiral strands the forming of typical A-form helices was not completed in the 500 ns time scale but the strands started to form bound states, and some of the 12


Page 12 of 29

Page 13 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 2: Base pairing types with their occurrence in the simulations (percentage of trajectory frames) for d-A12 /d-U12 and l-A12 /l-U12 systems. WC stands for the Watson-Crick and HG for the Hoogsteen edge, according to the classification of Leontis and Westhof. 40,41 The asterisk indicates that it was impossible to assign the WC or HG edge uniquely; this happens if only one hydrogen bond is formed between the bases and involves the indistinguishable hydrogen of amine groups of either A or U. end

nucleotides

end

5’

A1 | U24 A1 | U24 A1 | U24 A2 | U23 A3 | U22 A4 | U21 A5 | U20 A6 | U19 A7 | U18 A8 | U17 A9 | U16 A10 | U15 A11 | U14 A12 | U13

3’

3’

5’

base-pair type WC/WC WC*HG/WC*HG WC*HG/WC*HG WC/WC WC/WC WC/WC WC/WC WC/WC WC/WC WC/WC WC/WC WC/WC WC/WC WC/WC

pair configuration Cis Trans Cis Cis Cis Cis Cis Cis Cis Cis Cis Cis Cis Cis

% of frames d-A12 /d-U12 29.1 20.7 18.7 90.2 97.7 97.6 98.2 98.4 98.4 98.2 98.1 98.0 97.0 61.1

% of frames l-A12 /l-U12 31.7 13.4 16.6 90.8 97.8 98.0 98.3 97.1 98.2 96.7 97.5 97.6 97.1 73.0

nucleotides formed base pairs via the WC edges. Their starting conformations and the most representative clusters are shown in Figures S6 and S7, and the average hydrogen bond and van der Waals interaction patterns in Figure S8. For the d-A12 ..d-U12 system part of the stem formed hydrogen bonds via the WC-edges. The formation of the duplexes did not complete in this time scale also because the A/U oligomers are expected to have multiple secondary structures, e.g., with the WC-edge base pairs shifted by one or more nucleotides, which we indeed observed in these systems (Figs S6 and S7). The heterochiral strands also associated but their interaction patterns did not show hydrogen bonds formed via the WC edges of two opposite nucleotides in mirrored strands. The average number of hydrogen bonds created via WC and non-WC edges mapped on the strand sequences is shown in Fig. 8 confirming the absence of canonical pairing. However, the fraction of non-WC edge hydrogen bonds in these heterochiral complexes is substantial – contrary to homochiral helices (compare Fig. 8 with Fig. 6). 13



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 6: Secondary structures of a) d-A12 /d-U12 and b) l-A12 /l-U12 colored according to various descriptors calculated from the trajectories of the duplexes as averages per nucleotide: the number of hydrogen bonds formed via WC edges (Hbonds WC) and other edges (Hbonds non-WC), and van der Waals energy for the stacked bases (VDW Stacking, in kcal/mol). Nucleotides colored in red are characterized by the largest number of hydrogen bonds (both inter- and intra-strand) and the lowest van der Waals energy.

Stacking is another non-covalent interaction that stabilizes both double helix and single stranded structures of nucleic acids. 42 The calculated van der Waals contribution to stacking energies is presented in Fig. 8. These energy terms show that in the heterochiral complexes less bases stack than in homochiral duplexes. Overall, in heterochiral complexes the favorable interactions observed in the trajectories are non-canonical hydrogen bonds. Stacking is stronger and involves more bases in the homochiral complexes. Hydrogen bond patterns in heterochiral complexes MD simulations of heterochiral strands l-A12 ..d-U12 and l-U12 ..d-A12 indicate attractive interactions between l- and d-ribonucleotides since in each simulation some stable alignment was created. Representative conformations from clustering of the last 100 ns of the trajectories are shown in Fig. 9. The initially separated strands associated to form various types of interactions: base–base, base–sugar, and sugar–sugar. Also, we observed the interactions 14


Page 14 of 29

Page 15 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 7: RMSD in MD production runs for l-A12 ..d-U12 (blue), l-U12 ..d-A12 (red), lA12 ..l-U12 (cyan) and d-A12 ..d-U12 (green) simulated from separated strands. RMSD was calculated for heavy atoms and with the last simulation frame as a reference.

15



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 8: Nucleotides of a) l-A12 ..d-U12 and b) l-U12 ..d-A12 colored according to various descriptors averaged per nucleotide from the trajectories of initially separated heterochiral strands: the number of hydrogen bonds (both inter- and intra-strand) formed via WC edges (Hbonds WC), other edges (Hbonds non-WC) and van der Waals energy for the stacked bases (Stacking VDW, in kcal/mol).

16


Page 16 of 29

Page 17 of 29


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 10: Representative interaction patterns between heterochiral ribonucleotides obtained from the simulation of the l-A12 ..d-U12 strands. Hydrogen bonds are marked by green lines and their heavy atom distances are in [Å]. The numbering of patterns corresponds to the motif numbers listed in Tab. S5.

18


Page 18 of 29

Page 19 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 11: Representative interaction patterns between heterochiral ribonucleotides obtained from the simulation of the l-U12 ..d-A12 strands. Hydrogen bonds are marked by green lines and their heavy atom distances are in [Å]. The numbering of patterns corresponds to the motif numbers listed in Tab. S5.

19



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Hydration and ions around RNA complexes The distribution of water molecules and ions has proven important for the formation of RNA–ligand complexes (e.g. 43,44 ) and stability of nucleic acids. 45 In our trajectories we did not observe any water-mediated hydrogen bonds between the l-RNA and d-RNA oligomers. Nevertheless, even though the water distribution varies in the simulations, we detected some high water density areas, especially around the phosphate groups. Some high ion occupancies were also detected. Figure 12 shows high density areas of water and ions obtained in the simulations for the initially separated heterochiral systems (Figure S9 shows these densities for the corresponding homochiral systems). For example, we found high Na+ density in between the l-A12 and d-U12 strands. Also, some ions occupy the space in the middle of the homochiral complexes, probably precluding their further association in MD simulations. However, since the structural diversity of the geometric patterns is quite large, there are many possibilities of transient ion occupancies.

Conclusion We modified the AMBER ff10 for the simulations of l-ribonucleotides because in this force field the dihedral angles are not symmetric under mirror reflection. We described these modifications and included the ready to use files, along with installation scripts and examples, in the Supporting Information and GitHub repository (https://github.com/martadudek/lrnaamber/archive/master.zip). After testing the force field modifications for homochiral canonical duplexes, we performed MD simulations of l- and d-RNA oligonucleotides to gain first insight into their interactions. Trajectories pointed to some preferred hydrogen bond patterns between the strands of mixed chirality but these patterns did not involve canonical base pairing. Despite the standpoint presented in 21 that asserts certain kind of canonical pairing between l- and d-RNA, we did not observe such interactions between the heterochiral strands. The 20


Page 20 of 29

Page 21 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 12: Areas of high density of water molecules (light blue) and Na+ ions (yellow) derived from the last 100 ns of MD trajectories of the heterochiral strands. Densities are projected on the last-frame structures of: (a) l-A12 ..d-U12 and (b) l-U12 ..d-A12 . Light blue color corresponds to observed water density of ≥0.036 water oxygen atoms per Å3 and yellow ≥0.0036 Na+ ions per Å3 .

involvement of classical WC nucleotide edges 41 was marginal in the formation of heterochiral complexes. Therefore, we believe that structural motifs of l-/d- complexes require further careful examination. The trajectories cannot be considered exhaustive, i.e., there exist far more conformational patterns than we could observe in a 0.5 µs long series of MD simulations. The presented cases obviously did not cover all kinds of secondary motifs for l-U/d-A and vice versa. Solid resampling of conformational space, preferably with enhanced sampling techniques, is necessary, especially for systems containing guanines and cytosines, as well as mixed sequences. Also, accounting for different starting conformations of the separated strands is necessary. Further, simulations of the interactions of l-ribonucleotides with other biomolecules such as proteins would be valuable. We plan to continue the studies on intermolecular interactions of l-RNAs, also experimentally. Such microscopic picture of l-RNA/d-RNA interactions 21



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

might help in the future to formulate the rules for design of l-RNAs with desired affinity toward various d-RNA motifs.

Acknowledgement The authors thank National Science Centre for support (2013/11/N/NZ1/02384) to MD and (DEC-2012/05/B/NZ1/00035 and DEC-2014/12/W/ST5/00589) to JT. Simulations were performed at the Interdisciplinary Centre for Mathematical and Computational Modelling University of Warsaw (grant no G31-4 and GB65-28) and Centre of New Technologies UW.

Supporting Information Available • SI.pdf: section S1 with Tables S1–S3 describe modification of the Amber ff10 force field. Table S4 shows local base-pair parameters for homochiral helices. Table S5 gives hydrogen bonds distances for patterns observed in heterochiral complexes. Figures S1–S9 show data from analyzed trajectories. • leaprc.ff14SB: leap file that indicates the force field to load • parm10_LRNA.dat: leap file corrected for l-RNA • nucleic12_LRNA.lib: new leap library for l-RNA • nucleic10.in: leap file with definitions of l-nucleotides This material is available free of charge via the Internet at http://pubs.acs.org/.

References (1) Berman, H.; Battistuz, T.; Bhat, T.; Bluhm, W.; Bourne, P.; Burkhardt, K.; Feng, Z.; Gilliland, G.; Iype, L.; Jain, S.; Fagan, P.; Marvin, J.; Padilla, D.; Ravichandran, V.; Schneider, B.; Thanki, N.; Weissig, H.; Westbrook, J.; Zardecki, C. The Protein Data Bank. Acta Crystallogr., Sect. D: Biol. Crystallogr. 2002, 58, 1937–1949. 22


Page 22 of 29

Page 23 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(2) Vallazza, M.; Perbandt, M.; Klussmann, S.; Rypniewski, W.; Einspahr, H.; Erdmann, V.; Betzel, C. First look at RNA in L-configuration. Acta Crystallogr., Sect. D 2004, 12, 1–7. (3) Yatime, L.; Maasch, C.; Hoehlig, K.; Klussmann, S.; Andersen, G.; Vater, A. Structural basis for the targeting of complement anaphylatoxin C5a using a mixed L-RNA/L-DNA aptamer. Nat. Commun. 2015, 6, 6481. (4) Oberthür, D.; Achenbach, J.; Gabdulkhakov, A.; Buchner, K.; Maasch, C.; Falke, S.; Rehders, D.; Klussmann, S.; Betzel, C. Crystal structure of a mirror-image L-RNA aptamer (Spiegelmer) in complex with the natural L-protein target CCL2. Nat. Commun. 2015, 6, 6923. (5) Klussmann, S.; Nolte, A.; Bald, R.; Erdmann, V.; Fürste, J. Mirror-image RNA that binds D-adenosine. Nat. Biotechnol. 1996, 14, 1112–1115. (6) Wlotzka, B.; Leva, S.; Eschgfäller, B.; Burmeister, J.; Kleinjung, F.; Kaduk, C.; Muhn, P.; Hess-Stumpp, H.; Klussmann, S. In vivo properties of an anti-GnRH Spiegelmer: An example of an oligonucleotide-based therapeutic substance class. Proc. Nat. Acad. Sci. USA 2002, 99, 8898–8902. (7) Nolte, A.; Klussmann, S.; Bald, R.; Erdmann, V.; Fürste, J. Mirror-design of loligonucleotide ligands binding to l-arginine. Nat. Biotechnol. 1996, 14, 1116–1119. (8) Vater, A.; Jarosch, F.; Buchner, K.; Klussmann, S. Short bioactive Spiegelmers to migraine-associated calcitonin gene-related peptide rapidly identified by a novel approach: Tailored-SELEX. Nucleic Acids Res. 2003, 31, 1–7. (9) Helmling, S.; Maasch, C.; Eulberg, D.; Buchner, K.; Schröder,; Lange, C.; Vonhoff, S.; Wlotzka, B.; Tschöp, M.; Rosewicz, S.; Klussmann, S. Inhibition of ghrelin action in vitro and in vivo by an RNA-Spiegelmer. Proc. Natl. Acad. Sci. USA 2004, 101, 13174– 13179. 23



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(10) Leva, S.; Lichte, A.; Burmeister, J.; Muhn, P.; Jahnke, B.; Fesser, D.; Erfurth, J.; Burgstaller, P.; Klussmann, S. GnRH binding RNA and DNA Spiegelmers: a novel approach toward GnRH antagonism. Chem. Biol. 2002, 9, 351–359. (11) Faulhammer, D.; Eschgäfller, B.; Stark, S.; Burgstaller, P.; Englberger, W.; Erfurth, J.; Kleinjung, F.; Rupp, J.; Dan Vulcu, S.; Schröder, W.; Vonhoff, S.; Nawrath, H.; Gillen, C.; Klussmann, S. Biostable aptamers with antagonistic properties to the neuropeptide nociceptin/orphanin FQ. RNA 2004, 10, 516–627. (12) Eulberg, D.; Buchner, K.; Maasch, C.; Klussmann, S. Development of an automated in vitro selection protocol to obtain RNA-based aptamers: identification of a biostable substance P antagonist. Nucleic Acids Res. 2005, 33, 1–10. (13) Williams, K.; Liu, X.; Schumacher, T.; Lin, H.; Ausiello, D.; Kim, P.; Bartel, D. Bioactive and nuclease-resistant L-DNA ligand of vasopressin. Proc. Natl. Acad. Sci. USA 1997, 94, 11285–1290. (14) Ellington, A.; Szostak, J. In vitro selection of RNA molecules that bind specific ligands. Nature 1990, 346, 818–822. (15) Vater, A.; Klussmann, S. Turning mirror-image oligonucleotides into drugs: the evolution of Spiegelmer therapeutics. Drug Discov Today 2015, 20, 147–155. (16) Gaubert, G.; Babu, B.; Vogel, S.; Bryld, T.; Vester, B.; Wengel, J. Synthesis and RNAselective hybridization of alpha-L-ribo- and beta-D-lyxo-configured oligonucleotides. Tetrahedron 2005, 62, 2278–2294. (17) Jensen, T.; Madsen, A.; Pasternak, A.; Petersen, M.; Wengel, J. Synthesis and structural characterization of 2’-fluoro-alpha-L-RNA-modified oligonucleotides. ChemBioChem 2011, 12, 1904–1911.

24


Page 24 of 29

Page 25 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(18) Keinicke, L.; Sørensen, M.; Wengel, J. Alpha-L-RNA (alpha-L-ribo configured RNA): synthesis and RNA-selective hybridization of alpha-L-RNA/alpha-L-LNA chimera. Bioorg. Med. Chem. Lett. 2002, 12, 593–596. (19) Garbesi, A.; Capobianco, M.; Colonna, F.; Tondelli, L.; Arcamone, F.; Manzini, G.; Hilbers, C.; Aelen, J.; Blommers, M. L-DNAs as potential antimessenger oligonucleotides: a reassessment. Nucleic Acids Res. 1993, 21, 4159–4165. (20) Ashley, G. Modeling, Synthesis, and Hybridization Properties of (l-Ribonucleic Acid). J. Am. Chem. Soc. 1992, 114, 9731–9736. (21) Szabat, M.; Gudanic, D.; Kotkowiak, W.; Gdaniec, Z.; Kierzek, R.; Pasternak, A. Thermodynamic Features of Structural Motifs Formed by β-L-RNA. PLoS ONE 2016, 11, 1–11. (22) Schlick, T., Ed. Molecular Modeling and Simulation: An Interdisciplinary Guide; Springer, 1997. (23) Case, D.; Babin, V.; Berryman, J.; Betz, R.; Cai, Q.; Cerutti, D.; Cheatham, I. T.; Darden, T.; Duke, R.; Gohlke, H.; Goetz, A.; Gusarov, S.; Homeyer, N.; Janowski, P.; Kaus, J.; Kolossváry, I.; Kovalenko, A.; Lee, T.; LeGrand, S.; Luchko, T.; Luo, R.; Madej, B.; Merz, K.; Paesani, F.; Roe, D.; Roitberg, A.; Sagui, C.; Salomon-Ferrer, R.; Seabra, G.; Simmerling, C.; Smith, W.; Swails, J.; Walker, R.; Wang, J.; Wolf, R.; Wu, X.; ; Kollman, P. AMBER 14. University of California, San Fransisco. (24) Banáš, P.; Hollas, D.; Zgarbová, M.; Jurečka, P.; Orozco, M.; Cheatham, I. T.; Sponer, J.; Otyepka, M. Performance of Molecular Mechanics Force Fields for RNA Simulations: Stability of UUCG and GNRA Hairpins. J. Chem. Theory Comput. 2010, 6, 3836–3849. (25) Wang, J.; Cieplak, P.; Kollman, P. How Well Does a Restrained Electrostatic Poten-

25



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

tial (RESP) Model Perform in Calculating Conformational Energies of Organic and Biological Molecules? J. Comp. Chem. 2000, 21, 1049–1074. (26) Pérez, A.; Marchán, I.; Svozil, D.; Sponer, J.; Cheatham, I. T.; Laughton, C.; Orozco, M. Refinement of the AMBER force field for nucleic acids: improving the description of alpha/gamma conformers. Biophys. J. 2007, 92, 3817–3829. (27) Yildirim, I.; Stern, H.; Kennedy, S.; Tubbs, J.; Turner, D. Reparameterization of RNA chi Torsion Parameters for the AMBER Force Field and Comparison to NMR Spectra for Cytidine and Uridine. J. Chem. Theory Comput. 2010, 5, 1520–1531. (28) Zgarbová, M.; Otyepka, M.; Sponer, J.; Mládek, A.; Banáš, P.; Cheatham, I. T.; Jurečka, P. Refinement of the Cornell et al. Nucleic Acids Force Field Based on Reference Quantum Chemical Calculations of Glycosidic Torsion Profiles. J. Chem. Theory Comput. 2011, 9, 2886–2902. (29) Macke, T., Case, D., Eds. Modeling Unusual Nucleic Acid Structures; Am. Chem. Soc., 1997; pp 379–393. (30) Pettersen, E.; Goddard, T.; Huang, C.; Couch, G.; Greenblatt, D.; Meng, E.; Ferrin, T. UCSF Chimera–a visualization system for exploratory research and analysis. J. Comput. Chem. 2004, 13, 1605–1612. (31) Phillips, J.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R.; Kalé, L.; Schulten, K. Scalable molecular dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781–1802. (32) Ewald, P, 1921, Die Berechnung optischer und elektrostatischer Gitterpotentiale. (33) Berendsen, H.; Postma, J.; van Gunsteren, W.; DiNola, A.; Haak, J. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 1984, 81, 3684–3690.

26


Page 26 of 29

Page 27 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(34) Kräutler, V.; van Gunsteren, W.; Hünenberger, P. A fast SHAKE algorithm to solve distance constraint equations for small molecules in molecular dynamics simulations. J. Comput. Chem. 2001, 22, 501–508. (35) Dudek, M.; Romanowska, J.; Wituła,; Trylska, J. Interactions of amikacin with the RNA model of the ribosomal A-site: Computational, spectroscopic and calorimetric studies. Biochimie 2014, 102, 188–202. (36) Xiang-Jun, L.; Olson, W. 3DNA: a software package for the analysis, rebuilding and visualization of three–dimensional nucleic acid structures. Nucleic Acids Res. 2003, 31, 5108–5121. (37) Górska, A.; Jasiński, M.; Trylska, J. MINT: software to identify motifs and shortrangeinteractions in trajectories of nucleic acids. Nucl. Acids Res. 2015, 17, 1–10. (38) Darty, K.; Denise, A.; Ponty, Y. VARNA: Interactive drawing and editing of the RNA secondary structure. Bioinformatics 2009, 25, 1974–1975. (39) Humphrey, W.; Dalke, A.; Schulten, K. VMD: visual molecular dynamics. J. Mol. Graph. 1996, 14, 33–38. (40) Leontis, N.; Westhof, E. Geometric nomenclature and classification of RNA base pairs. RNA 2001, 7, 499–512. (41) Leontis, N.; Stombaugh, J.; Westhof, E. The non-Watson-Crick base pairs and their associated isostericity matrices. Nucleic Acids Res. 2002, 30, 3497–3531. (42) Hobza, P. Stacking interactions. Phys. Chem. Chem. Phys. 2008, 10, 2581–2583. (43) Vaiana, A.; Westhof, E.; Auffinger, P. A molecular dynamics simulation study of an aminoglycosides/A-site RNA complex: conformational and hydration patterns. Biochimie 2006, 88, 1061–1073.

27



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(44) Kulik, M.; Goral, A. M.; Jasiński, M.; Dominiak, P. M.; Trylska, J. Electrostatic Interactions in Aminoglycoside-RNA Complexes. Biophys. J. 2015, 108, 655–665. (45) Auffinger, P.; Hashem, Y. Nucleic acid solvation: from outside to insight. Curr. Opin. Struct. Biol. 2007, 325–333.

28


Page 28 of 29

Page 29 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Graphical TOC Entry

29


Molecular Dynamics Simulations of l-RNA Involving Homo- and

Recommend Documents