key: cord-294196-3gox5art
authors: Jin, Hongwei; Zheng, Suxin; Wang, Zhanli; Luo, Cheng; Shen, Jianhua; Jiang, Hualiang; Zhang, Liangren; Zhang, Lihe
title: Structural insights into the effect of isonucleosides on B-DNA duplexes using molecular-dynamics simulations
date: 2006-02-01
journal: J Mol Model
DOI: 10.1007/s00894-005-0085-8
sha: 
doc_id: 294196
cord_uid: 3gox5art

Some structural insights into the conformations of the isonucleosides containing duplexes have been provided. Unrestrained molecular-dynamics simulations on 18-mer duplexes with isonucleosides incorporated at the 3'-end or in the center of one strand have been carried out with explicit solvent under periodic boundary conditions using the AMBER force field and the particle mesh Ewald method. The Watson–Crick hydrogen-bonding patterns of the duplexes studied remained intact throughout the simulation. For the modified duplexes, the changes observed in the inter-base pair parameters and backbone torsional angles were primarily localized at the isonucleoside-inserted area. All five structures studied remained in the B-form family. The decreased stacking abilities indicated by the large changes in inter-base pair parameters and the large changes in backbones made the modified duplexes show a minor thermal destabilization in comparison with native DNA. The MM_PBSA method for estimating binding free energies on two complementary strands was used. The results showed that the binding free energies of isonucleoside-incorporated DNA duplexes were lower than the native DNA duplex, which is in good agreement with experimental observations. [Image: see text]

Synthetic oligonucleotides and their analogs have attracted considerable interest because of their potential therapeutic and diagnostic applications [1] [2] [3] [4] [5] . The natural oligonucleotides, however, are usually degraded relatively easily by cellular nucleases and cannot be used directly as therapeutics. It is therefore a challenge to develop modified antisense oligonucleotides that are degradation resistant, and can hybridize selectively to target nucleic acids with strong affinity. Over the past decades, a large number of oligonucleotides with modifications in the phosphodiester, sugar, or nucleobase have been prepared. Among them, phosphorothioate, which is the first generation of antisense compounds, shows increased resistance to degradation both in vitro and in vivo [6] [7] [8] [9] , and one drug has been admitted to market by the FDA, which represents a first indication of the feasibility of the therapeutic antisense concept. Meanwhile, many antisense phosphorothioates are in phase I/II clinical trials [10, 11] . However, phosphorothioate oligonucleotides often show decreased specificity and affinity to the target sequences [12] and inefficient cellular uptake. Thus, other analogs, such as phosphoramidates [13] , PNA [14] , LNA [15] , α-L-LNA [16] , HNA [17] , ANA [18] , CAN [19] , CeNA [20] , and INA [21] , have been studied systematically. They have proved to be capable of hybridizing with complementary DNA or RNA with increased thermal stabilities compared to the parent DNA:DNA or DNA:RNA duplexes and can be valid candidates for clinical applications.

We have reported the synthesis and enzymatic stability of oligonucleotides incorporated with isonucleosides [22] [23] [24] [25] . Isonucleosides are a new class of nucleoside analogues in which the nucleobase is linked to various positions of ribose other than C1′. Previous research has shown that a single isonucleotide-incorporated oligonucleotide antagonized the enzymatic hydrolysis of snake-venom phosphodiesterase strongly, although the binding ability of the modified sequence to the complementary sequence was slightly decreased, which indicated that there was a possibility to design a stable antisense oligonucleotide by incorporating isonucleotides in the sequence [26] . For the detailed understanding of structural information of isonucleoside-modified antisense oligonucleotides, an 18-mer of the antisense sequence 5′-AACATCTCCTGAGGGA AC -3′ incorporated with isonucleosides A 1 (6′-OH free) or A 2 (6′-OH protected by allyl group), which was complementary to the mRNA region encoding spike glycoprotein of the severe acute respiratory syndrome-associated coronavirus (SARS-CoV) (22, 415 bp) has been studied by detailed unrestrained molecular dynamics in aqueous solution (Fig. 1) . The results, which agree with the experimental data, provide an insight into the local and global structural features of the duplexes containing A 1 and A 2 . The structural and energetic analysis from the MD simulations provides some information for further research and development of isonucleotide-incorporated antisense oligonucleotides.

All MD simulations were performed with the AMBER 7 molecular simulation package [27] . The AMBER 99 force field was used to describe the DNA. To obtain molecular mechanical parameters for the isonucleosides A 1 and A 2 , ab initio quantum chemical methods were employed using the Gaussian 98 program. The geometries were fully optimized using the AM1 Hamiltonian and the electrostatic potentials around them were determined at the HF/6-31G** level. The RESP strategy [28] was used to obtain the partial atomic charges.

Starting models of DNA duplexes were built in the B canonical structures using the Insight II package [29] . All constructed oligonucleotide duplexes were solvated in TIP3P water using a rectangular box, which extended 10 Å away from any solute atom. This yielded about 4,400 water molecules used for solvation. To neutralize the negative charges of simulated molecules, Na + counterions were placed next to each phosphate group.

Molecular dynamics (MD) simulations were carried out using the SANDER module of AMBER 7. At first, the water box and counterions were subjected to a series of equilibration MD runs while holding the DNA fixed [30, 31] . The equilibration runs began with 1,000 steps of minimization with a large constraint of 500 kcal mol −1 Å −2 on the DNA atoms and were followed by 25 ps of MD, reducing the constraints to 100 kcal mol −1 Å −2 , during which the temperature was slowly raised from 100 to 300 K over 1 ps and was maintained at 300 K for the remaining 24 ps. The size of the box was allowed to change until the water density and pressure converged to the correct values. Subsequent equilibration steps, during which position constraints on the DNA molecules were gradually relaxed, as well as the final production runs, were done by using the particle mesh Ewald (PME) method to calculate electrostatic interactions [32, 33] . First, another 25 ps of MD were performed while still holding the DNA fixed to relax the solvent molecules fully and to complete the density equilibration. This was followed by a second set of 1,000 steps of minimization and 3 ps of MD, which were carried out with constraints on the DNA molecules of 25 kcal mol −1 Å −2 . Finally, five rounds of 600 steps of energy minimization were performed, in which positional constraints were reduced to 5.0 kcal mol −1 Å −2 in each round.

After this initial equilibration, the whole system was heated from 100 to 300 K for 20 ps without positional constraints. The production simulations of 1 ns for S0:S, 2 ns for S1:S, and 1.5 ns for S2:S, S3:S, or S4:S were then started at constant pressure (1 atm) and temperature (300 K). In the entire simulation, SHAKE was applied to all hydrogen atoms [34] . Periodic boundary conditions with minimum-image conventions were applied to calculate the nonbonded interactions. A cutoff of 9 Å was used for the Lennard-Jones interactions. The nonbonded pairs were updated every 30 steps. An integration time step of 2 fs was used and the trajectory was saved every 100 steps for future analysis.

The MD trajectories were analyzed using AMBER 7 and in-house software. The final structure of each duplex was produced from the 1,000 steps of minimized averaged structure of the last 800 ps of MD. Helix parameters were calculated with the program CURVES 5.3 [35] . All calculations were performed on a 2-CPU SGI Octane workstation.

Free-energy analysis was performed using the MM_PBSA scripts supplied by AMBER 7.0 [36] [37] [38] . The total free energy can be estimated from the molecular mechanical energy, the solvation free energy, as well as the solute entropies for a series of snapshots and then averaging the results. The binding free energies of all double strands were investigated using the single-trajectory approach. Singletrajectory results mean that the thermodynamics data are extracted from a single trajectory of the complex. The assumption here is that the natural DNA, the isonucleosideincorporated DNA, and the complementary DNA conformations are the same in the double strand and in the dissociated form and thus the vibrational, rotational, and translational entropy contributions can be neglected. Esti- mation of binding free energies in this manner has proven successful previously [36, 37, 39] . The binding free energy (ΔG binding ) between the two complementary strands is calculated as the difference between the electrostatic interaction energies, van der Waals interaction energies, the solvation free energies, as well as the entropies of the double strands and the dissociated system: We applied the same force field and parameter set used in the MD simulations and an electrostatic constant 1 to compute the electrostatic energy (E elec ) and the van der Waals energy (E vdw ) but without cutoff for nonbonded interactions.

The solvation free energy (ΔG solv ) was estimated from the electrostatic solvation energy (ΔG PB ) and the nonpolar solvation energy (ΔG nonpolar ):

The electrostatic contribution to the solvation free energy (ΔG PB ) was estimated with a Poisson-Boltzmann electrostatic continuum method using the program Delphi II [40] . The dielectric boundary is the molecular surface defined by a 1.4 Å probe sphere and by spheres centered on each atom with radii taken from the PARSE parameter set [41] . We used an interior dielectric of unity, and the outside dielectric was set to 80. For the Poisson-Boltzmann calculation, a cubic lattice with linear dimensions 80% larger than the longest dimension was applied with a 0.5 Å grid spacing; potentials at the boundaries of the finite-difference lattice were set to a sum of Debye-Huckel potentials [42] . The nonpolar contribution to the solvation free energy (ΔG nonpolar ) was estimated from a surface-area dependent term using the solvent accessible surface area (SASA) algorithm and program of Sanner [43] with a parametrization of

where γ=0.00542 kcal Å −1 and β=0.92 kcal mol −1 . Solute entropic contributions (−TΔS) were estimated based on a harmonic approximation to the normal modes and standard (quantum) formulae at 300 K [44] . Snapshots from the MD trajectories of the natural and isonucleosideincorporated DNA-DNA hybrids with water and counterions removed were considered for the binding-free-energy calculations. A total of 80 snapshots were selected at 10 ps intervals from each of the 800 ps trajectories.

The five unrestrained MD simulations of the duplexes studied led to stable trajectories of geometry and energy terms. Calculations of root-mean-square deviations (RMSDs) with respect to the starting B-form structures confirmed the stability of the trajectory. Five plots of the RMSD values as a function of the simulation time are shown in Fig. 2 .

For the native DNA duplex, a fast convergence of RMSD was observed. The RMSD rose to 2.2 Å within about 100 ps, and averaged around 2.5 Å in the final stage of the simulation. The RMSD of duplex S1:S and S3:S trajectories also showed a sharp increase to about 3.5 Å in the first 200 ps, and then remained stable for the rest of the simulation (average RMSD around 4.0 Å). The RMSD values of S2:S and S4:S duplexes showed a significantly slow increase during the first 700 ps and 500 ps, respectively, which indicated a movement and structural change of the molecules. Then the RMSD values stabilized and remained stable, oscillating around their average values up to the 1.5 ns of the MD simulation. Therefore, only conformations generated during the last 800 ps period of simulations were used to monitor the structure properties. The comparison of the time evolution of the RMSDs showed that the duplexes of S2:S and S4:S were more difficult to achieve conformational equilibrium than the other ones. For the S0:S duplex, a slightly smaller average RMSD value (2.65 Å) was obtained than for the S1:S (3.68 Å), S2:S (3.23 Å), S3:S (3.70 Å), and S4:S (3.35 Å) duplexes (Table 1) , which indicated that S0:S was closer to the started B-form conformation than the other four modified duplexes.

An important indicator of the stability of duplex structures is the length of hydrogen bonds and percentage occupancy during the MD simulations. A summary of Watson-Crick hydrogen-bond lengths and their standard deviations between the base pairs in each of the five duplexes during the last 800 ps of the MD is given in Table 2 . The hydrogen bonds in any of the five structures were well maintained and the majority was 100% occupied during the simulation. The results showed clearly that there were no obvious differences in hydrogen bonding properties among S0:S, S1:S, S2:S, S3:S, and S4:S duplexes.

The data showed that the average distances of N1(17)… H3 (20) and H61 (17)…O4(20) hydrogen bonds for S0:S were 1.98 and 2.01 Å, respectively. The corresponding values for S1:S and S2:S were 2.01 and 2.01, 2.09 and 1.98 Å, respectively. Moreover, the average distances of N1 (12)…H3(25) and H61(12)…O4(25) hydrogen bonds for S0:S, S3:S, and S4:S were 1.99 and 1.96, 1.99 and 2.05, 2.08 and 1.92 Å, respectively. The lengths of hydrogen bonds for base pairs A 1 -T and A 2 -T in the modified duplexes were quite similar to those in the native DNA. Therefore, the position change of nucleobase in isonucleosides A 1 or A 2 did not obviously affect the hydrogen bonding of the Watson-Crick base pairs in the modified DNA duplexes.

There are three major categories of helical parameters, i.e., axis-base pair, intra-base pair and inter-base pair parameters Table 1 The RMSD values from averaging the trajectories for the five duplexes studied Table 3 . Values of the parameters slide, rise, roll, and twist along the sequence for the MD-averaged structures of the studied duplexes are plotted in Fig. 3 . 

For the five duplexes studied, the helix parameters Xdisp, Ydisp, inclination, and tip were closer to the canonical Btype DNA than to the canonical A-type DNA. There was a significant variation in the values of Ydisp and tip along the strand of modified structures S1:S, S2:S, S3:S, and S4:S. This was related to the slightly bent helix axis in these structures (Fig. 4) .

The calculated values for shear, stretch, and stagger of five DNA structures did not show much deviation from B-type DNA. The buckle changes significantly with a range of −20 to 30°along the strand. The values of propeller twist and opening all deviated from either A-type or B-type DNA. Especially, large positive values of the opening parameter were observed for all five duplexes, and indicated opening toward the major groove. A more noticeable parameter was the propeller twist (the relative rotation between the two-paired bases). Propeller twist up to −30°is native [46] to the AT base pair of DNA. The calculated values for propeller twist varied significantly along the strand with a range of −25 to 10°.

The shift, slide, tilt, and roll parameters did not deviate much from either B-or A-DNA and were within values normally observed in either A-type or B-type DNA. The rise parameter, which indicated vertical displacement of one base pair with respect to the other, remained in the range for canonical B-DNA (average value 3.38 Å) for all base pair steps of the five DNA duplexes. All were more characteristic of B-DNA than A-DNA. They remained at around 3.29, 3.46, 3.33, 3.32, and 3.46 Å for S0:S, S1:S, S2:S, S3:S, and S4:S, respectively. For all duplexes, the values of twist changed along the strands. However, for the majority of the base pairs of each duplex, values even lower than that of A-DNA were observed. The average twists over the sequence were 32±1°. The value was more characteristic of A than B. This showed particularly that the modifications did not cause the helix to unwind. In general, it could be concluded from the calculated results of all helix parameters that the overall conformations of the five duplexes remained in the B-conformations during the entire course of the simulations. It is known that the global intra-base pair parameters describe the relative displacement and orientation of the bases in the Watson-Crick base pairs. These parameters of the native S0:S duplex obtained during the simulation were very similar to those found in the four modified DNA duplexes (especially in the modified base-pair steps) according to our calculations. The parameters also explained why there were no fundamental differences in hydrogen-bonding properties among duplex S0:S, S1:S, S2:S, S3:S, and S4:S. However, for the inter-base pair parameters, several showed large differences in the modified area, which might indicate the effect of the isonucleosides on the conformation of the neighboring bases.

The A 1 -T20/A16-T21 (0.78 Å), A 2 -T20/A16-T21 (0.97 Å) and A 1 -T25/G11-C26 (0.63 Å), A 2 -T25/G11-C26 (1.04 Å) base-pair steps in S1:S, S2:S and S3:S, S4:S structures, respectively, exhibited significantly larger slide values than the corresponding A17-T20/A16-T21 (0.12 Å) and A12-T25/G11-C26 (0.20 Å) base-pair steps in S0:S. The rise parameters showed large deviations for the A 1 , A 2 , and their neighboring base-pair steps in the S2:S, S3:S, and S4:S duplexes from corresponding native S0:S values. Wherever the isonucleoside A 1 or A 2 was located at the 3′end or in the center of the sequence, the slide and rise values of the A 2 base pair steps was larger than those of the A 1 ones. The A 1 -T20/A16-T21 (29.1°) or A 2 -T20/A16- T21 (21.2°) base pair steps in the S1:S or S2:S duplexes, respectively, showed larger roll values than in the S0:S duplex. In the S3:S and S4:S duplexes, the isonucleosides A 1 and A 2 induced changes in the roll of the steps below and above but not in the steps connected by A 1 and A 2 themselves. As shown in Fig. 3 , the twist angles of S0:S, S2:S, and S3:S did not show obvious changes along the sequence. However, for S1:S and S4:S duplexes, significant differences were observed in the twist values. The lowest twist angle was found for the A 1 -T20/A16-T21 step (12.5°) in S1:S structure at this base pair step. And the C6-G31/T5-A30, T10-A27/C9-G26 and G13-C24/A 2 -T23 steps in the S4:S structure also showed small twist values lower than 20°. These small twist values were an indicator of untwisting of the DNA helix at the corresponding sites. Such small twist values for the G13-C24/A 2 -T23 step might be caused by the incorporation of isonucleoside A 2 , but we need much more studies to understand why the T10-A27/C9-G26 and G13-C24/A 2 -T23 steps also show so small twist values. Inter-base-pair parameters were used to describe the stacking interactions of base pairings. Thus, the unusually large changes in these parameters for the isonucleoside area decreased the stacking interactions of base pairings in the modified S1:S, S2:S, S3:S, and S4:S duplexes. These decreased stacking interactions made the duplexes incorporated with isonucleosides a little less stable than the unmodified DNA molecule, in agreement with the energy comparison.

For the S0:S duplex, the average values of the torsional angles over all sequences were close to those of the standard canonical B structure. Only small differences of all torsional angles were observed for the S0 and S strands. Most of the S0:S backbone torsions were restricted to move in a relatively narrow range around the optimum value. The range of the dynamic averages and the RMSD fluctuations for the torsional angles indicated that the backbones of S0: S were quite rigid (data not shown). The presence of the isonucleosides A 1 and A 2 produced large deviations in backbone torsional angles in the modified area in comparison with those in the corresponding area of the natural DNA duplex (Table 4) , which resulted in visible kinks in the backbone but did not have a significant impact on the whole backbone conformation of the modified duplexes. The rest of the sugar-phosphate backbone remained in the range expected for the B-DNA family.

The sugar puckers for the native residues in all five duplexes possess C1'-exo conformations. Minor distortions Table 4 Averaged values of the backbone torsional angles (in degrees) for A12 and A17 in the S0 strand, A 1 in the S1 and S3 strands, and A 2 in the S2 and S4 strands a Torsional angles S0(A17) S1(A 1 ) S2(A 1 ) S0(A12) S3(A 2 ) S4(A 2 ) Table 5 .

The results indicate that the backbones of native DNA bear higher structural stability than those of modified oligonucleotides in the neighborhood of the modified sequence. The large changes in the backbone torsional angles of the isonucleoside region result in large deviations in some of the inter-base-pair parameters for the modified base pair step, and hence minimize the degree of hybridization of the modified duplexes.

The side views of the five DNA structures studied with their global axis are shown in Fig. 4 , which demonstrate that all duplexes have the features of B-DNA. The C6′-OH and C6′-OCH 2 =CH-CH 3 groups in the isonucleosides A 1 and A 2 point toward the major groove. The average structures of the isonucleosides incorporated region in the S1:S, S2:S, S3:S, and S4:S duplexes are shown in Fig. 5 .

Binding free energies of the five duplexes studied

The binding free-energy contributions of all DNA duplexes studied are listed in Table 6 . The calculations are based on single trajectories of the duplexes.

The calculated value for binding of the S0:S hybrid is −23.6 kcal mol −1 , more stable by 2.1, 2.0, 7.7, and 7.9 kcal mol −1 , respectively, than the isonucleoside-modified S1:S, S2:S, S3:S, and S4:S hybrids. The binding free energy of S1:S (−21.5 kcal mol −1 ) is close to that of the S2:S (−21.6 kcal mol −1 ). Those of S3:S (−15.9 kcal mol −1 ) and S4:S (−15.7 kcal mol −1 ) agree similarly. However, the incorporation of A 1 and A 2 modifications at the 3′-end of the strand (duplexes S1:S and S2:S) create lower binding free energies by 5.6 and 5.9 kcal mol −1 , respectively, compared to the duplexes with isonucleosides in the center of the sequence (duplexes S3:S and S4:S). UV thermal stability measurements [26] show that the oligonucleotide consisting of isonucleoside A 1 or A 2 decreases the affinity Fig. 5 Views of the average structures of the isonucleosideincorporated region. The isonucleosides A 1 and A 2 are indicated in red. The A16-A 1 -C18/G19-T20-T21 in the S1:S (left, the structure of base C18 distorted because of the end base pair opening) and the A16-A 2 -C18/ G19-T20-T21 in the S2:S (right) are shown at the top, the G11-A 1 -G13/C24-T25-C26 in the S3:S (left) and the G11-A 2 -G13/C24-T25-C26 in the S4:S (right) are shown at the bottom Values from [26] for the complementary DNA. A 2 -containing oligonucleotides do not show obvious affinity differences to the complementary DNA compared to corresponding A 1containing sequences. However, the duplex in which A 1 (A 2 ) is at the 3′-end of the sequence has higher thermal stability than the duplex in which A 1 (A 2 ) is in the center of the sequence. The calculated results and experiment are in good agreement.

Comparing the van der Waals/nonpolar (ΔE vdw /ΔG nonpolar ), the electrostatic (ΔE elec +ΔG PB ) and the solute entropic contributions (−TΔS), we find that the association between the two strands in the five duplexes studied is mainly driven by more favorable van der Waals interactions. As shown in Table 6 , the van der Waals, nonpolar and solute entropic contributions in the natural DNA hybrid are close to those in the other four isonucleoside-incorporated hybrids. In the S0:S duplex, the electrostatic contribution ÁG binding elec decreases by 3.5, 2.2, 11.2, and 7.7 kcal mol −1 , respectively, in comparison with those in the S1:S, S2:S, S3:S, and S4:S duplexes, which causes a l repulsive force and thus further leads to the highest binding affinity in the natural S0:S duplex. The electrostatic contributions are lower for the S1: S and S2:S than for the S3:S and S4:S by 7.7 and 5.5 kcal mol −1 , respectively. The unfavorable strand-strand electrostatic repulsive interactions make the S3 and S4 strands produce worse affinity with the S strand, in comparison with S1 and S2 stands.

Overall, the binding free-energy differences observed between the five DNA duplexes studied are mainly driven by the different electrostatic interactions between complementary strands. For the S1:S, S2:S, S3:S, and S4:S duplexes, more unregularity of the structures on the modified region (compared with the natural DNA) make some atoms approach each other closely, causing a large repulsive force and thus further contributing to the unfavorable thermal stabilization.

Unrestrained molecular-dynamics simulations were performed on a 1.0-2.0 ns scale for an 18-mer DNA duplex S0:S and isonucleoside-incorporated duplexes S1:S, S2:S, S3:S, or S4:S, in water with full periodic boundary conditions and Na + counterions using the particle mesh Ewald method. The results show that the S2:S and S4:S duplexes achieve conformational equilibrium with difficulty. Watson-Crick hydrogen bonds remained intact throughout the simulations and there were no fundamental differences in the hydrogen-bonding properties of the S0:S, S1:S, S2:S, S3:S, and S4:S duplexes.

The backbone torsional angles and the helicoidal parameters support the view that the overall structures of all five duplexes are quite close to the canonical B-type DNA. For the four modified duplexes S1:S, S2:S, S3:S, and S4:S, the large changes in the inter-base-pair parameters and backbone torsional angles for the modified area in comparison with natural S0:S demonstrate that dramatic alterations of structures occur because of the incorporation of the two isonucleosides. These alterations produce decreased stacking abilities and large electrostatic repulsive forces and thus lead to unfavorable thermal stabilization for the modified DNA duplexes.

The MM_PBSA method was used to compute binding free energies of two complementary strands of the DNA duplexes studied. The calculated values showed that the natural DNA duplex has a lower binding free energy than the duplexs containing isonucleosides and the duplex in which A 1 (A 2 ) is at the 3'-end of the sequence has a lower binding free energy than the duplex in which A 1 (A 2 ) is in the center of the sequence, which were in good agreement with experimental results. The differences in binding free energies between five duplexes are mainly driven by the different electrostatic interactions between complementary strands.

Antisense research and applications

Nucleosides and nucleotides as antitumor and antiviral agents

Statistical mechanics