E2EDNA: Simulation Protocol for DNA Aptamers with Ligands

2.2.1. Parameterization

One step which may be required before running E2EDNA is force field parameterization. Depending on the analyte under study, one may require novel force field parameters, not already found in the literature. In the case of the polarizable AMOEBA force field which we employ here, one may generate parameters using the automated parameterization tool, POLTYPE[16], or its successor POLTYPE 2. POLTYPE is an open-source program, which was recently used e.g., to parameterize a variety of charged nucleosides[17]. With only an input structure, POLTYPE outputs a structure and accompanying quantum mechanically fit AMOEBA parameters, ready for simulation on the Tinker molecular dynamics platform.

POLTYPE 2 has made some improvements over the original POLTYPE in torsion and van der Waals parameterization. To achieve faster and more accurate parameterization, a fragmenter is implemented to derive torsion parameters. Bulky molecules often have steric interactions while sampling the conformational surface that lead to artificial barriers on the dihedral energy surface. Thus, fragmenting is important for accurately fitting torsion parameters. The fragmenter will start with a minimal fragment containing the atoms a-b-c-d plus neighbors of atoms, while not cutting any aromatic ring systems. The fragment will keep growing until a convergence criterion is met. That criteria being either hitting a maximum number of growth steps or relative Wiberg Bond Order between fragment and parent. Each fragment molecule follows the original POLTYPE procedure to derive multipoles, and torsion parameters from a rigid dihedral scan. Parameters from fragments are transferred back to the parent molecule.

If implementing this pipeline with a different MD platform or force-field, the parameterization tool would change accordingly, and similar utilities such as CGenFF and AmberTools Paramfit exist for the popular CHARMM and AMBER force fields respectively [18,19]. It is worth noting now that the MD platform, force-field, and virtually all other the individual components of the E2EDNA are modular; they could each be easily replaced if e.g., one desires e.g., a less expensive force field, or a different structure prediction protocol.

2.2.2. Secondary Structure Prediction

A crucial difficulty in aptamer-analyte binding analysis is the pre-folding of the aptamer to the correct equilibrium structure. In general, a single DNA sequence may adopt a very wide array of folded structures [20], and brute force approaches such as naïve molecular dynamics search are prohibitively computationally expensive.

Despite decades of work on this problem, due to the complexity of the underlying physical system, secondary structure prediction is not always possible with high confidence with existing tools. Indeed, in sequences with length exceeding a few dozen bases, de novo prediction without experimental input becomes extremely difficult, as the number of possible stable structures may become large, and possible structures may be difficult to discriminate [11]. Common software packages approach this issue by issuing an ensemble of predicted structures with associated probabilities, when appropriate, or by adding base-by-base pairing probabilities to overall structure predictions. One may also solicit predictions from various packages which each employ different numerical methods such as template-based models, fragment assembly, and minimization of empirical potentials [11].

E2EDNA leverages such purpose-built models to a-priori predict a structure or set of structures a given aptamer has a realistic probability of adopting, before analyzing and refining down to a final structural prediction. In the absence of experimental data, we assess the quality of predicted structures via molecular dynamics simulation with an appropriate force-field, similar to approaches which have been used in RNA structure evaluation[21]. MD sampling then allows us to compare the stability of proposed structures on the nanosecond-microsecond timescale.

In our E2EDNA implementation, we primarily use the ‘NUPACK’ and ‘seqfold’ Python packages for secondary structure prediction [22], which implement a range of modern empirical energy models to identify the minimum free energy structure at a given temperature and ionic strength. Both are remarkably fast and easy-to-use, with straightforward I/O and highly useful utility functions. One has only to supply the sequence FASTA string and associated simulation conditions.

2.2.3. 3D Structure Prediction

Given a proposed aptamer secondary structure in the form of a list of paired bases, we perform a directed simulation to fold it from an extended initial condition into the prescribed secondary structure. For this task, we use the program, MacroMoleculeBuilder (MMB) [23], which folds the structure through inexpensive simulation, with user-scalable attractive fictitious forces between the paired bases. Using such a tool, we can fold from an arbitrary DNA sequence to a 3D structure which agrees with the predicted secondary structure in a matter of minutes.

Depending on the size and complexity of a given aptamer structure, a customized user-specified annealing schedule may be required. In the SI we supply a script with a robust protocol which is efficient even for multi-hairpin structures exceeding 60 bases in length, along with a script which automatically generates MMB input files from secondary structure and sequence information.

Given the strong fictitious forces used to fold the structure in MMB, the output of this simulation is only a rough prediction of the aptamer 3D structure, which is analyzed and refined by subsequent MD simulation. Beyond refining the details, extended MD sampling is also used to verify the stability of a predicted secondary structure. For this refinement, we use the MD protocol detailed in Section 2.2.5.

2.2.4. Docking Site Prediction

Once a suitable representative 3D aptamer structure has been isolated, before binding analysis, one must complex the aptamer with the analyte. As in the case of structure of the aptamer itself, it is generally wasteful to randomly initialize the DNA and analyte initial positions, and search for an appropriate binding configuration by exhaustive molecular dynamics sampling. Instead, we a-priori predict one or more possible binding configurations and use these as the initial condition for binding trajectories. This prediction is made according to some docking protocol, which suggests likely analyte positions binding orientations, from which we initiate our sampling trajectories. Such a protocol could take the form of a docking program which predicts a set of likely binding sites using conformational and/or electrostatic features [24–29].

For sequences sampled from an experimentally obtained distribution of sequences that bind a particular analyte, a data driven approach may be taken [28,29], which takes as inputs a selection of aptamer sequences from a SELEX run on the target analyte. From such a procedure, we identify certain bases whose mutation is likely to be important for determining aptamer-analyte binding efficiency, and therefore sites where it may be profitable to attempt to dock the analyte molecule.

2.2.5. Molecular Dynamics Sampling

Given a set of likely docking configurations between the analyte and folded aptamer, we sample the complex via molecular dynamics simulation. We generate trajectories with and without the analyte and compare the equilibrium probability densities along a series of user-identified reaction coordinates, with the difference between complexed and free aptamer equilibria determining the influence/affinity of the analyte with the aptamer. The same simulation and analysis protocols enable us to analyze and evaluate aptamer 3D structures, in the absence of the analyte.

In order to accurately capture the nuanced dynamics of the high charge of the analyte in our test case (see Section 3), we employ the advanced electrostatic and polarizable force field, AMOEBA[30]. We used the Tinker9 software suite running on Compute Canada NVIDIA V100 GPUs for the bulk of our simulations, as well as for simulation setup and post-processing. Detailed simulation parameters, derived from POLTYPE 2, in the form of Tinker keyfiles are given in the SI.

2.2.6. Binding Analysis

To assess the influence of analyte binding to the aptamer, we track the time evolution of user-specified reaction coordinates, identified by observing the dynamics of free and complexed aptamers. In the test case examined in Section 3, these coordinates describe base-pairing, large-scale structural rearrangements, and the distance from the analyte to predicted binding sites. Given sufficient sampling time we may compute the equilibrium density, P(r), and accompanying free energy profile, F(r) along the reaction coordinate and interrogate the stability of the folded structure and the influence of the analyte on the aptamer complex.

The free energy profile, as a function of the reaction coordinate distance, is given by,

The probability distribution P(r) is computed via histogram of time-series analysis of the relevant reaction coordinates r, which are tracked using the Python package, MDAnalysis. Selection of reaction coordinates is done on an aptamer-by-aptamer basis, though at least partial automation of this procedure is possible.