NAFlex: a web server for the study of nucleic acid flexibility.

We present NAFlex, a new web tool to study the flexibility of nucleic acids, either isolated or bound to other molecules. The server allows the user to incorporate structures from protein data banks, completing gaps and removing structural inconsistencies. It is also possible to define canonical (average or sequence-adapted) nucleic acid structures using a variety of predefined internal libraries, as well to create specific nucleic acid conformations from the sequence. The server offers a variety of methods to explore nucleic acid flexibility, such as a colorless wormlike-chain model, a base-pair resolution mesoscopic model and atomistic molecular dynamics simulations with a wide variety of protocols and force fields. The trajectories obtained by simulations, or imported externally, can be visualized and analyzed using a large number of tools, including standard Cartesian analysis, essential dynamics, helical analysis, local and global stiffness, energy decomposition, principal components and in silico NMR spectra. The server is accessible free of charge from the mmb.irbbarcelona.org/NAFlex webpage.

INTRODUCTION

Nucleic acids are polymorphic molecules whose conformations depend not only on the sequence but also on external factors, such as temperature, solvent, ionic environment and presence of ligands. Sequence and environment-dependent polymorphisms make the experimental description of nucleic acids extremely difficult. Our coverage of the nucleic acid structural space in the near future is not likely to approach that already available for proteins (1–3). The experimental information regarding flexibility is even more scarce, being partial and very much limited to the B-type DNA duplex (1–3). Only recently, NMR spectroscopy has started to provide dynamical information at the atomistic level for some model nucleic acid systems (4–6). Although the results are of impressive quality and impact, the technique is clearly still unable to provide a complete description of the general flexibility of nucleic acids.

In the absence of an experimental approach, simulation techniques are now widely accepted tools with which to describe nucleic acid structure and flexibility (7–11). According to PUBMED, in 2012, nearly 900 articles were published quoting the words ‘DNA’, ‘RNA’ or ‘nucleic acids’ in combination with ‘simulation’ or ‘molecular dynamics’. Given the improvements in the simulation methods and the accessibility to increasingly faster computers, greater popularization of the field is expected in the near future (9).

Theoretical approaches for the study of nucleic acids are diverse, but can be classified on the basis of two basic parameters: (i) the nature of the Hamiltonian used; and (ii) the level of resolution (7,8,11,12). The simplest simulation approaches, such as the wormlike-chain (WLC, 13) model, assume simple Hamiltonians that describe global helical properties of the DNA fiber. Such methods were developed to examine general properties of very long fibers of canonical B-DNA with great computational efficiency. Atomistic molecular dynamics (MD) takes a completely different approach. It combines a complex physical Hamiltonian with a rigorous method based on the resolution of Newton’s equations of motion to obtain trajectories. MD allows an accurate atomistic description of any nucleic acid (8,9), but at the expense of a very large computational cost. The applicability of MD for very large systems or for the analysis of processes happening over large time scales is very limited. Between WLC and MD lie a wide variety of mesoscopic simulation methods. These were designed for the study of medium-to-large DNA fibers at the base-pair or even pseudo-base level of resolution using Hamiltonians of intermediate complexity (7,8,11–16).

For many years, the simulation of nucleic acids was performed by a small number of expert groups, often the developers of the software, force fields or algorithms. As the popularity of the simulation techniques has increased, groups with limited knowledge of the theoretical tools and often of the physics behind them have become users of these approaches. This development has generated significant confusion, as the potential user now faces a plethora of simulation packages based on a wide range of physical models and dealing with different levels of resolution. Furthermore, in general, a nonexpert user can find little information on how to setup the systems and encounters significant problems when attempting to perform meaningful analysis of the trajectories collected.

Here, we present NAFlex, a web tool designed to facilitate the use of nucleic acid simulation tools for newcomers to the field. The server allows the introduction of nucleic acid structures from diverse sources, as well as automatic structure generation from the sequences. The structural models can be subjected to several simulations, these are based on: (i) a colorless WLC model (13); (ii) a base-pair resolution mesoscopic model (14,15); or (iii) atomistic MD simulations. In the latter case, MDWeb technology (17) is used to help the user during all the setup and equilibration steps, providing all the input files required to launch the simulations. Finally, the trajectory obtained (or uploaded) can be visualized and analyzed using a large set of nucleic acid-specific tools. The server is freely accessible from the mmb.irbbarcelona.org/NAFlex webpage.

MATERIALS AND METHODS

NAFlex is divided (see Figure 1) into three main blocks: (i) Input; (ii) Simulation engines; and (iii) Analysis. The webserver has been designed to obtain a maximum coverage of potential user needs and as such is extremely flexible at the three levels (input, simulation and analysis).

Figure 1.

Representation of the NAFlex web server overview showing its three main blocks: (i) the inputs accepted: Sequence, structure and trajectory; (ii) the simulation engines implemented: molecular dynamics, mesoscopic simulations and WLC calculations; and (iii) the set of nucleic acid-specific flexibility analysis and visualization tools offered.

Input

The user can introduce nucleic acid information in many ways depending on the nature of the problem and the calculation planned (see Figure 1). Thus, it is possible to upload a Protein Data Bank (PDB) structure (18) (of either isolated or complexed nucleic acids), or any user-derived structural model. The user can also upload a trajectory (all usual formats are accepted), which is then sent directly to the analysis modules of the server. When the user plans to work with canonical DNA and/or RNA duplexes, NAFlex allows the automatic generation of double helices using the Nucleic Acid Builder (nab) program from the Ambertools package (19). Structure generation can use parameters from various sources: (i) fiber diffraction data (20,21); (ii) sequence-dependent average X-ray structural information (1); (iii) tetramer-dependent average MD-simulation results (22–24); and (iv) user-defined helical values. The structures generated can be relaxed later to remove any distortion that may occur during model generation before being used in mesoscopic or atomistic simulations. In the case of WLC calculations, B-DNA geometry is assumed, and only sequence information is required as input.

After input, the server checks the structure to correct for potential gaps or missing atoms. Before launching the simulation, the setup procedures and a series of quality controls (see server help) are applied to the model. The server’s check list includes, among others, alternate atom/residue locations, unusual distances between consecutive bases, steric clashes and the presence of metal ions or modified nucleotides/ligands. The server warns the user about potential structural errors (not trivial to correct).

Simulation engines

NAFlex offers the nonexpert user a selected variety of simulation tools that differ in complexity and resolution and that are designed to make their use as simple as possible.

WLC model (13) was introduced using the formalism of Jian et al. (25,26), implemented into a Monte Carlo sampling procedure. Accordingly, the DNA is represented as a set of N beads, each comprising M base-pairs (typically M = 10 base-pairs, but the exact bead resolution is selected by the user), and the potential energy of the DNA is defined by a simple Hamiltonian (see the help section for additional details): where the stretching energy (ES) is computed as: where K is the stretching constant, l is the actual distance between beads and l is the optimum bead–bead distance. Based on (25,26), the value of K was set at 100 k/l02 (as it reproduces the correct DNA bond variance), with k the Boltzmann constant and T the temperature. The bending energy is determined as: where is the bending constant and is the Euler bending angle between the local coordinate systems of two consecutive beads; and the torsion potential is defined by: where K stands for the torsional rigidity, and the sum of α and γ Euler angles defines the torsion between the local coordinate systems of two consecutive beads. ϕ is a parameter that gives the mean DNA twist, and it depends on the bead–bead equilibrium distance and the helical repeat. The values of the bending and torsional rigidity constants are related to the DNA bending (P) and twist (C) persistence lengths, respectively, and the bead–bead equilibrium distance: and

The electrostatic repulsion energy (), the only nonlocal term in the Hamiltonian, is determined using Debye–Hückel potential: where ν is the salt-dependent Stigter’s effective DNA linear charge density (27), D is the dielectric constant of water, κ is the inverse Debye length and r is the distance between two beads. Alternatively, the user can define a value for the effective DNA-bead charge, qDNA, and account for the electrostatic potential as follows: where ε0 is the electric permittivity of vacuum, and ε is the dielectric constant (set to 80).

The elastic mesoscopic model provides an intermediate level of resolution and potential energy complexity (7,8,11,12,14–16). It assumes that DNA deformations can be approximated as the addition of harmonic distortions of equilibrium base-pair step geometries. Three rotational (twist, roll and tilt) and three translational (slide, shift and rise) degrees of freedom are considered, thereby allowing us to define the Hamiltonian as: where k is the Boltzman constant, T is the absolute temperature, E is the energy associated with the deformation ΔX and k stands for the different stiffness constants defined by the 36 elements of the stiffness matrix (Ξ) [twist (w), roll (r), tilt (t), rise (s), slide (l) and shift (f)]. The Ξ can be calculated (see equation 9) by inversion of the covariance matrix obtained from either analysis of MD trajectories (at dinucleotide or tetranucleotide level, 22–24) or from the analysis of dinucleotide step variability in the crystal structures of DNAs and DNA–protein complexes (1,15). In all the cases, sampling is obtained using Monte Carlo simulations in the helical coordinate space following the DNALive (28) protocol.

Atomistic molecular dynamics (MD) is the most flexible, universal and accurate simulation engine incorporated to NAFlex. Unfortunately, it is also the most expensive, thus limiting its practical use to medium-sized nucleic acid structures (7–10). After the quality control check, and before running a MD simulation, the system is neutralized, immersed in the desired solution, minimized, thermalized and equilibrated in a few steps (Supplementary Figure S1). Once the equilibration is finished, the server provides NAMD (29), AMBER (30) or GROMACS (31) adapted input files, which can be then used to perform the production run on the user’s local systems (see server help). Output trajectories can be uploaded back to the server for analysis.

Analysis tools

NAFlex integrates a variety of analysis packages for mining nucleic acid trajectories. Analysis results are presented in several numerical and graphical formats. Available analyses include the following: standard Cartesian and helical analysis (17,32); principal component analysis (PCA) (33,34) to determine the essential deformation movements and their associated stiffness; and helical stiffness analysis, performed following Olson–Lankas’s approach (14,15). A basic energy analysis is conducted to determine stacking and hydrogen bonding interaction energies along the trajectory. Finally, NAFlex also provides estimates of a number of NMR observables (see help for details), which can be used for either trajectory validation or for ab initio spectral predictions.

SERVER IMPLEMENTATION AND USAGE

NAFlex is a web portal implemented in PHP and MySQL. It provides users with a personal workspace where intermediate data, trajectories and results of analysis can be stored. Access to the server is free; however, registration is required to maintain a permanent workspace. Nonregistered users can use the full functionality of the server, although they need to download the results at the end of the session. In the present implementation, registered users are provided with 2 Gb of storage space. Storage for specific projects can be increased on demand. NAFlex is powered by the MDWeb platform, which was designed for automated MD simulations of proteins (17). MDWeb provides the necessary modules for data management, structure checking, MD simulation setup and standard Cartesian trajectory analysis (see http://mmb.irbbarcelona.org/MDWeb/). NaFlex incorporates a variety of tools specifically oriented for nucleic acids, including input facilities, simulation engines and analysis tools (see ‘Materials and Methods’ section).

Workspace

The starting points for the analysis could be single structures, uploaded as PDB files or derived from a given sequence, or previously simulated trajectories (see above). The input structure or trajectory is used as initial step to initiate a project in the personal workspace. There is no limit to the number of initiated projects while the permitted storage capacity is not exceeded. For each project, the workspace holds all intermediate structures and results, organized as a tree view (Supplementary Figure S2), thereby allowing the user to track the history of operations performed. For each entry of the tree, a number of operations, selected on the basis of the type of data, are offered to the user. For all intermediate results, a download tool and a series of JMol-powered (http://www.jmol.org) visualization tools are available. Complete projects can also be downloaded for local storage. Downloaded projects can be restored in NAFlex again at a later time, thus recovering the original workspace.

Simulations

Coarse-grained simulations from initial structures can be performed within the NAFlex environment. When atomistic MD simulations are done, a series of preparation tools is required. Specific setup procedures adapted to nucleic acids are available, although experienced users may design and incorporate their own protocols into the workflow. Running short MD simulations for testing purposes is allowed on the server; however, production simulations should be run locally in the user’s own facilities using NAFlex-generated input files (see above). The user can then upload trajectories into his/her workspace for analysis.

NAFlex provides a wide repertoire of specific tools for the analysis of nucleic acids, some of them local, others incorporated from external sources. All the tools are offered in a common interface, which requires no additional expertise for its use. The user has complete control over the level of resolution in the analysis through an interactive duplex viewer (Figure 2a). He/she can also analyze time course or average values and study in detail the resulting structural models (see Figure 2b and c). When available, standard values obtained from the literature are included for comparison purposes (Figure 2b). Plots can be opened in separate windows to compare the structural or mechanical behavior of different regions of the molecules. In all cases, raw data can be downloaded for local analysis. It should be noted that full analyses are done upon the initial request. This allows the user to browse the results even for a large structure without noticeable delay. Trajectories can be manipulated to extract individual snapshots for further inspection (Figure 2d).

Figure 2.

Screenshots of the analysis results interface. (a) Interactive duplex viewer allowing the level of resolution where the analysis is focussed: base-pairs (|), base-pair steps (X) and nucleotide pair steps (-). (b) Sequence distribution of time-averaged Roll, compared with MD- and X-ray averaged values. (c) Time course of the Roll at 6-GGCC base-pair step, with an accessory histogram (mean and standard deviations are also displayed). (d) JMol interactive representation of a nucleic acid structure. Graphical visualization is offered for structures and also for trajectories.

Please note that trajectory upload is limited by network bandwidth. The current limit for NAFlex-uploaded trajectories is 100 Mb (as indicated in the corresponding help pages). However, in our experience, once the trajectory has been stripped of solvent and ions, and taking a representative ensemble of snapshots, the analysis of the resulting trajectory with NAFlex can provide useful information about the general flexibility of the molecule. See Example 4 in the next section as an illustration of a real analysis.

EXAMPLES OF USE

NAFlex is an extremely flexible server that offers a large number of options, and it does not require deep knowledge of either simulation engines or the physics of nucleic acids. A few examples of use are included here as references. These examples can be also accessed online on the NAFlex server using the ‘demo’ account.

Example 1: Atomistic MD from a nucleotide sequence

The files required to run an atomistic MD simulation can be obtained following a robust protocol. First of all, the user should define the type of input (Supplementary Figure S3a), in this case, DNA/RNA Simulation From sequence. The only inputs required are the nucleotide sequence (either typed or from a FASTA-formatted file) and the desired set of DNA/RNA helical parameters (Supplementary Figure S3b). For instance, to generate a DNA molecule with a higher value for twist than the canonical B-DNA, User-defined DNA should be chosen at DNA/RNA type selector and the Twist value should be modified accordingly (in the example from 36.0° to 38.0°). The structure will then be generated using the Nab program from the Ambertools package (19). A new window will appear showing the newly created structure in a JMol applet, together with the results of a set of structural checks (Supplementary Figure S3c and d). After validation of the original structure, a new workspace for the project is prepared (see Supplementary Figure S3e). The interface allows the user to select the required operation, in this case a complete setup for an AMBER MD simulation (see Supplementary Figure S3e) using the ff99SB (35) and PARMBSC0 (36) force fields with recent corrections [ildn (37) and OL3 (38)]. These are the recommended settings for a wide set of structures, including protein–DNA complexes and RNA; however, users have the choice of other force-field combinations. The evolution of the running workflow can be visualized in real time (Supplementary Figure S3f), and the final files can be downloaded for local production runs (Supplementary Figure S3g).

Example 2: Mesoscopic simulation of DNA

This example illustrates a mesoscopic simulation of a long protein-free DNA molecule at the base-pair level of resolution. As in the previous example, the first steps are the definition of the type of input (DNA/RNA Simulation From sequence), the selection of a title for the project and the introduction of the sequence (see Supplementary Figure S4a and b). In this case, DNA/RNA type should be set to Coarse-Grained DNA Model (Base Step Level), to select the desired coarse-grained resolution level. The prepared structure can be visualized using the JMol applet (Supplementary Figure S4c). To launch the simulation, the user should select the desired number of ensemble snapshots (500 in this example), while the type of mesoscopic Hamiltonian (Coarse-grained DNA Elastic Mesoscopic Model) is already defined by the resolution level (see Supplementary Figure S4d). After the Monte Carlo ensemble has been collected, the snapshots can be further analyzed to trace the flexibility of the fiber being studied. In this particular example, the probability of cross-talk between distant segments of the DNA is analyzed. This can be done by simple inspection of the ensemble (Supplementary Figure S4e) and from the averaged nucleotide-contact maps (see Supplementary Figure S4f). The appearance of short contacts (red) off-diagonal is indicative of long-range contacts of potential biological relevance.

Example 3: Nucleic acid flexibility analysis

The example illustrates the analysis of a 100-ns MD simulation of a DNA dodecamer (dCGCGAGGACGCG). The trajectory in this case was uploaded from an external source. Analysis is initiated by selecting the Nucleic Acid Flexibility Analysis operation in the Analysis Tools box (Supplementary Figure S5a). The first results to be inspected are the helical characteristics of the duplex, obtained from a CURVES (32) analysis. Results (Supplementary Figure S5b and c) show the presence of a subtle Twist bimodality (Time-course Twist plot, Supplementary Figure S5b) and a marked Roll increase (taking X-ray as reference) in the first CpG step (Averaged Roll plot, Supplementary Figure S5c). This observation suggests that this step can be especially flexible and can display spontaneous curvature. Additional analysis shows that the DNA is highly robust in terms of helical structure, with strong hydrogen bonds and stacking interactions (see Supplementary Figure S5d); however, both ends show fraying effects (see, Supplementary Figure S5e), which are expected also to impact Nuclear Overhauser Effects (NOEs) (see Supplementary Figure S5f).

Example 4: Protein–DNA complex, human mitochondrial transcription factor A as a case study

The example shows the analysis of an important protein–DNA complex. The structural experimental information available (39–41) suggests that flexibility is crucial for the functionality of the complex. Analysis is initiated from a user-provided 100-ns trajectory obtained with GROMACS 4 (31). In the example, a 1-ns time window is retained for analysis.

A JMol visualizer provides initial structural views of the trajectory. Thus, in Figure 3a, the human mitochondrial transcription factor A (TFAM) protein can be clearly identified in dark pink cartoon representation, with the two high mobility group (HMG) protein domains joined by the inter-domain linker helix. The 22-base-pair oligonucleotide mitochondrial light strand promoter (LSP) attached to the TFAM protein is shown in purple. Once uploaded, the trajectory is automatically split into protein and DNA, allowing separate analysis. An extensive flexibility analysis of the protein counterpart (TFAM) can be performed using FlexServ platform (42), also available through the NAFlex workspace.

Figure 3.

Protein–DNA complex TFAM–LSP as a case study. (a) JMol representation of the TFAM–LSP complex. To facilitate inspection, the protein part is represented by a dark pink cartoon, while the nucleic part is shown by a purple cartoon. (b) Average inter-base-pair helical parameters Roll and Twist computed on the LSP oligonucleotide plotted together with standard values obtained from the literature for comparison purposes. (c) Snapshot of the animation of the first mode of the DNA obtained from PCA. (d) B-factor analysis of the TFAM protein.

Basic helical analysis (see Figure 3b) shows that roll and twist angles are completely distorted along the trajectory in the regions where the protein residues intercalate into DNA base-pairs (Figure 3a). PCA was used to obtain a clear picture of the essential movements of the oligonucleotide in the complex. In this example, the first three eigenvectors explain most of the variance (Figure 3c and interactive JMol online). The first two modes involve the kinking of the complete LSP along the direction of the linker helix, bringing both DNA ends closer together, thus adopting a U-turn. The third mode shows the twisting of the nucleic acid caused by the motion of the two TFAM domain boxes HMG1 and HMG2. Domain motions can be visualized using the available JMol Applet. In accordance with the hypothesis that the HMG domains induce an overall DNA U-turn stabilized by the helix linker (40,41), FlexServ analysis of the protein moiety (Figure 3d) reveals the presence of flexible segments near the domain boxes, while the linker helix appears to be relatively rigid.

CONCLUSIONS

The field of nucleic acid simulation has reached maturity. It has moved away from the times when only ‘proof of concept’ simulations were done by a very small number of highly specialized groups toward a ‘full production’ situation, where nucleic acid simulation tools are used by many groups, often not experts in the physics of nucleic acids or the theory behind the simulation package. NAFlex is a bioinformatics tool created to facilitate the use of simulation tools for nonexpert users interested in gaining insight into the dynamics of nucleic acid systems.

SUPPLEMENTARY DATA

Supplementary Data are available at NAR Online: Supplementary Figures 1–5.

FUNDING

Funding for open access charge: [BIO2012-3286]; European Research Council Advanced Grant (ERC); ; Consolider E-Science Project; EU-Scalalife project; Fundación Marcelino Botín; European Union’s Seventh Framework Programme (FP7/2007–2013) [275096 to R.C.-G. and M.O.].

Conflict of interest statement. None declared.