Loading of DOX into a tetrahedral DNA nanostructure: the corner does matter.

With the rapid development of nanotechnology, various DNA nanostructures have been synthesized and widely used in drug delivery. However, the underlying mechanisms of drug molecule loading into the DNA nanostructure are still elusive. In this work, we systematically investigate the interactions of a tetrahedral DNA nanostructure (TDN) with the anti-cancer drug doxorubicin (DOX) by combining molecular docking and all-atom molecular dynamics simulations. It is found that there are five possible binding modes in the single TDN-DOX interactions, namely the outside-corner mode, the inside-corner mode, the major-groove mode, the minor-groove mode, and the intercalation mode, where the van der Waals (VDW) interaction and the electrostatic (ELE) interaction dominate in the case of unionized DOX and ionized DOX, respectively. Moreover, with the increase of the DOX number, some of the interaction modes may disappear and the inside-corner mode is the most energy-favorable mode. The present study enhances the molecular understanding of the role of TDN as the drug carrier, which may provide a useful guideline for the future design of DNA nanostructures. This journal is © The Royal Society of Chemistry.

Introduction

Nucleic acids play an irreplaceable role in biological inheritance since most forms of life use DNA as a carrier of genetic information. In nature, there are two typical forms of DNA structures, namely the single-stranded (ss) form and the double-strand (ds) form that comes from two ssDNA molecules using the hybridization principle. The ss-DNA or ds-DNA can further fold or assemble into 3D structures like nucleosomes, which determine their biological functions in living organisms. Besides the most common B-form dsDNA, the DNA molecule has a variety of other structures.

With the growing understanding of the structure of DNA molecules, recently the programmability of the DNA primary structure and richness of the secondary structure have attracted more and more attention from researchers in nanoscience and biochemistry.[1-6] The advent of various assembly methods like DNA origami or DNA tiles made it possible to compose complicated 3D structures like DNA origami nanotubes,[7,8] nanoflasks,[9] nanocages,[10,11]etc. In the meantime, the various types of DNA-based nanomaterials synthesized have been widely used in biomedical applications.[12-19]

Drug delivery could be one of the most intriguing applications, where the DNA nanostructure is used as a carrier to deliver the drug molecules in the human body.[20-23] Apart from increasing the solubility of drugs, protecting the drugs from enzyme attack as well as fascinating cellular delivery of the drugs (which other nanocarriers may also have), the good biocompatibility[24] and the low cytotoxicity[25] make the DNA nanostructure an ideal candidate for drug delivery. In the meantime, thanks to its editable characteristics, modifications for targeting can be easily applied. For example, some delivery systems or detection platforms have been developed based on DNA nanocages including tetrahedrons[26-29] and octahedrons.[30-32]

Despite the extensive research on the in vitro/in vivo delivery of drug molecules by DNA nanostructures in the experiments,[19,33-37] it is still largely unknown how the DNA structures interact with the drug molecules, which requires a deeper understanding of their interactions at the molecular level. More importantly, due to the unique topological structure of the DNA nanostructures, more possible drug-binding sites may occur, which may be quite different from the interaction of a single dsDNA chain with drug molecules.[38-43]

In this work, we take the tetrahedral DNA nanostructure (TDN) and doxorubicin[44-46] (DOX) as examples of the DNA nanostructure and drug molecule, respectively, and investigate their interactions under different conditions by combining molecular docking and all-atom molecular dynamics simulations. As shown below, there are at least five different binding modes in the TDN–DOX interaction, which are dependent on the ionization state and the concentration of the DOX. The binding energy analysis along with its energy decomposition is further used to reveal the underlying molecular mechanism of the different binding behaviors.

Model and methods

The atomistic model of TDN-17 was built using the Polygen program,[47] where 17 indicates that each side of the triangle that makes up the tetrahedron consists of 17 nucleotides including 16 hybridized nucleotides and 1 free nucleotide (Table 1). The initial conformation of DOX was downloaded from the PubChem database (PubChem CID: 31703).[48] Then the geometry optimization of the DOX molecule was calculated by using the B3LYP functional together with a 6-311G(d,p) basis set, and was carried out with the Gaussian 09 program.[49]

Sequences of nucleotides

Face	Sequence
F1	AAGGCGAGGAATTTTTCCTGAGGGTGCTACCTGTCGGGTGTTTAAAGTTAT
F2	GAAAATCTCGGCCTCTTGCCCACAAATTTCAATGCTACGAGAGAAGTGGCA
F3	CTTTTAGAGCCGGAGAGGTTACATCCAGGCTACTTTCCGCTCCTTAAAAAT
F4	GACTCCCACGATGGACTCAATGTAGGTCCGATGTGATGCTCTCTTCACCGG

The HDOCK webserver was used for molecular docking.[50,51] The optimized DOX along with a stabilized TDN-17 was submitted to the webserver (http://hdock.phys.hust.edu.cn/), where a pre-equilibrated TDN-17 was used as the receptor molecule and the DOX was used as the ligand molecule. After the docking job was completed, the server gave 100 pre-generated binding models based on the docking score and the root-mean-square deviation (RMSD) of the ligand. The top 30 predictions were examined manually based on the score ranking and cluster analysis to evaluate all the possible binding modes.

All-atom molecular dynamics (MD) simulations were further performed based on the above-selected structures. Notably, each simulation system only considered one DOX monomer/aggregate and one TDN. All MD simulations were carried out by the GROMACS 2019.6 software package,[52,53] and the Amber14sb_Parmbsc1 (ref. 54 and 55) force field was used to describe the properties of DNA and ions together with the TIP3P water model. NaCl was used to neutralize the system, and the concentration was set at 0.15 M, close to that under physiological conditions. The parameters for the DOX were generated by the general Amber force field[56,57] (GAFF) with the restrained electrostatic potential[58] (RESP) charges. The LINCS constraints[59] were used for all bonds involving hydrogen atoms. The particle mesh Ewald (PME) method was used when calculating the long-range electrostatic interactions, and the Lennard–Jones (LJ) interactions were cut off at a distance of 1.0 nm. The periodic boundary conditions were adopted in all three directions.

In the simulation, the system was firstly energy-minimized by the steepest descent method until convergence was reached. Then each system was subjected to a 500 ps pre-equilibration process in the NVT and NPT ensembles, separately, with all the heavy atoms harmonically constrained by 1000 kJ mol−1 nm−2. Finally, the constraint was released and 100 ns free NPT simulations were performed. The temperature was controlled at 300 K by the V-rescale thermostat[60] with a time constant of 0.2 ps and the pressure was kept at 1.0 bar by the Parrinello–Rahman barostat[61] with a time constant of 2.0 ps. The typical snapshots of this work were all drawn using VMD software.[62,63] The schematic diagrams of TDN–DOX interactions were drawn using PyMOL.[64]

In the binding energy calculation, the last 10 ns of each stabilized simulation trajectory was used with an interval of 100 ps, namely 100 frames for each run using the molecular mechanics Poisson–Boltzmann surface area (MMPBSA) method.[65-69] The gmx_mmpbsa script (https://jerkwin.github.io/gmxtool) together with the APBS[70,71] program were used to calculate the binding energy and the corresponding energy decomposition term, namely the electrostatic (ELE) energy, the van der Waals (VDW) energy, and the solvation energy including the polar part (PB) and the non-polar part (SA) in this work.

Results and discussion

Binding of a single unionized DOX molecule to the TDN-17

A typical dsDNA molecule is formed by two helical strands of ssDNA. The nucleobases on the two strands complement each other to form the secondary structure of DNA. Since the double-helical structure is not spatially symmetric, there are two different grooves defined by spatial distance (i.e., the major groove and the minor groove, Fig. 1a), and the dsDNA may have two different binding sites for the DOX (Fig. 1b). Moreover, the TDN is a tetrahedron composed of 6 B-DNA double helices and 6 single-base staples. Thus, the TDN has a richer spatial structure than the regular dsDNA and may provide more binding sites for the DOX, which makes the TDN a promising drug carrier. As shown in Fig. 1c–g, five typical binding modes were observed by the molecular docking, namely the inside-corner mode, the outside-corner mode, the intercalation mode, the major-groove mode, and the minor-groove mode. Notably, the intercalation mode, the major-groove mode, and the minor-groove mode have been reported in previous DNA–DOX experiments.[72,73] Moreover, a recent experimental study found that the DOX-loading ability of the TDN was much higher than that of the dsDNA with the same bps.[74] Considering that the corner is the only difference between the TDN and dsDNA for the DOX loading, the efficient loading of DOX at the corner of the TDN should be the point, which gives indirect evidence for the two new modes (i.e., inside-corner and outside-corner modes) here.

Fig. 1

(a) Schematic illustration for a single DOX binding to the TDN-17 and (b) the structural formula of the DOX. The possible binding modes observed in the molecular docking: (c) the inside-corner mode, (d) the outside-corner mode, (e) the intercalation mode, (f) the major-groove mode, and (g) the minor-groove mode.

To test the robustness of all the binding modes, we then performed all-atom molecular dynamics simulations for each complex above, respectively. As shown in Fig. 2, the root-mean-square deviation (RMSD) of TDN, the contact surface area (CSA) and the minimum distance between the DOX and TDN changed little with simulation time, indicating that the binding of DOX to TDN was relatively stable in all five modes.

Fig. 2

The root-mean-square deviation (RMSD) of TDN-17 (a), the contact surface area (CSA) between the TDN-17 and DOX (b), and the minimum distance between TDN-17 and DOX (c) as a function of time during the molecular dynamics simulations.

Having demonstrated the stability of the observed modes, we then calculated the binding energy of TDN–DOX and analyzed the binding driving force in each mode (Fig. 3a). Among the five modes, the binding energy was the lowest in the intercalation mode (−182 kJ mol−1), followed by the outside-corner mode (−169 kJ mol−1), the inside-corner mode (−162 kJ mol−1), the minor-groove mode (−150 kJ mol−1), and the major-groove mode (−102 kJ mol−1). Notably, the VDW energy was very strong in the intercalation (−212 kJ mol−1) and outside-corner modes (−190 kJ mol−1) since there was an obvious π–π stacking between the adenine of DNA and an anthraquinone ring part of DOX (Fig. 3b and c). Notably, there were some unpacked nucleotides in the outside corner of the TDN, which can provide a more flexible region for π–π stacking. As a result, the DOX with aromaticity was more likely to adsorb on the outside of the TDN. In contrast, the electrostatic interaction was strong in the inside-corner, minor-groove, and the major-groove modes since there were some hydrogen bonds formed by the nitrogen atoms (blue) in DNA and oxygen atoms (red) in DOX (Fig. 3d–f). Nevertheless, there were many polar residues at the binding site, leading to a high solvation energy result, thus the binding energy was weak in these cases. In general, the VDW interaction was the dominant driving force for the binding of single unionized DOX to the TDN.

Fig. 3

(a) The total binding energies and corresponding energy decomposition terms for the five binding modes. Structural depiction of key interfacial interactions in the five binding modes: (b) the intercalation mode, (c) the outside-corner mode, (d) the inside-corner mode, (e) the minor-groove mode, and (f) the major-groove mode. The red thick line indicates the π–π stacking and the blue dashed line represents the polar interaction.

Binding of multiple unionized DOX molecules to the TDN-17

Due to its aromatic ring, the DOX molecules can aggregate with each other in water when the concentration of DOX is not very low. Could the aggregation of DOX affect its binding to the TDN?

To answer this problem, we first simulated the self-assembly of two DOX molecules, four DOX molecules, and eight DOX molecules in water, respectively. After a 50 ns simulation, the DOX dimer, tetramer and octamer were obtained (Fig. S1†). Similar to the previous case, molecular docking was used to obtain the possible binding modes, followed by a 100 ns all-atom molecular dynamics simulation. The final snapshots for the possible binding modes are shown in Fig. 4a–c. Due to the formation of the aggregate, the size of the DOX became larger (than that of the single DOX), and there was not enough room for the DOX to insert into the intercalation space between two base pairs as well as the minor groove. As a result, the intercalation and the minor-groove modes can no longer be observed in this case. With a further increase of the DOX number, the outside-corner mode disappeared in the case of the DOX tetramer and octamer.

Fig. 4

Final snapshots illustrating the typical modes in the multiple DOX–TDN-17 interaction: (a) DOX dimer, (b) DOX tetramer, and (c) DOX octamer. (d) The total binding energies and corresponding energy decomposition terms for the above binding modes.

We then calculated the binding energies of TDN–DOX in these cases. As shown in Fig. 4d, in the case of the DOX dimer, the binding energy was the lowest in the inside-corner mode (−217 kJ mol−1), followed by the major-groove mode (−146 kJ mol−1) and the outside-corner mode (−72 kJ mol−1). Notably, the lowest binding energy here was a bit stronger than that in the intercalation mode of single DOX (−182 kJ mol−1). Moreover, the binding energy (−72 kJ mol−1) in the outside-corner mode here was even weaker than that in the case of single DOX (−169 kJ mol−1). With further increase of the DOX number, the lowest binding energy decreased gradually and was about −300 kJ mol−1 in the case of the DOX octamer. Actually, the formation of the DOX aggregates was mainly caused by the π–π stacking of the DOX molecules, which may shield some aromatic rings (from π–π stacking with the adenine in the TDN). As a result, the binding energy was not enhanced obviously. Moreover, the VDW interaction also dominated the binding process, which is similar to that in the single DOX.

Notably, the lowest binding energy of the multiple DOX–TDN interactions was observed in the inside-corner mode, which is the unique mode in the TDN (due to its 3D topological structure) and cannot be observed in the case of a single DNA chain. More importantly, the DOX was inside the TDN, which may protecte it from the undesired interactions during the in vivo cellular delivery. Notably, previous studies[75,76] indicated that the corner-attack is an efficient way for DNA nanostructures to minimize electrostatic repulsion and translocate through the negatively charged cell membrane. Here, the corner of the TDN again exhibited its merit and the TDN was indeed a good nanocarrier for drug delivery.

Binding of ionized DOX molecules to the TDN-17

The pKa of doxorubicin is about 8.0, indicating that it can be partially ionized in the physiological environment (pH ∼ 7.4), namely, the amino group at the end of the DOX molecule is protonated into a positively charged one, leading to the whole molecule being positively charged. How did the ionization of DOX affect the binding behaviors?

The possible binding modes along with their binding energies are listed in Fig. 5a and the corresponding final snapshots are shown in Fig. S2.† Similar to the case of unionized DOX, there still existed five, three, two, and two modes for the single ionized DOX, ionized DOX dimer, ionized DOX tetramer, and ionized DOX octamer, respectively. Nevertheless, the binding energies and the ranking of the binding modes became totally different. Since DNA is a highly negatively charged biomolecule and the ionized DOX carries one positive charge, there should be a strongly attractive electrostatic interaction between them. As a result, the binding energies were greatly enhanced. For example, the lowest binding energy of a single DOX–TDN was −398 kJ mol−1 and it was over twice as high as that in the unionized case (−182 kJ mol−1). Interestingly, the lowest binding energies were all observed in the inside-corner mode, which is probably due to the fact that the DOX molecules may contact more DNA bases inside the TDN.

Fig. 5

(a) The total binding energies as a function of the DOX number for the possible binding modes in the case of ionized DOX; (b) the energy decomposition terms as a function of the DOX number in the inside-corner binding mode.

We further analyzed the energy terms of the inside-corner mode under different DOX numbers (Fig. 5b). Not surprisingly, the ELE energy dominated in all the cases and increased with the DOX number in a nearly linear manner. As a result, the binding energy increased obviously with the DOX number, which was different from the slow increase of the binding energy (with the DOX number) in the case of unionized DOX.

Conclusions

In summary, we have investigated the binding of DOX molecules to TDN-17 under different conditions by combining molecular docking and all-atom molecular dynamics simulations. Due to the unique structure and extra space of the TDN, five possible binding modes in the TDN–DOX interactions are observed, namely the outside-corner mode, the inside-corner mode, the major-groove mode, the minor-groove mode, and the intercalation mode. The binding energies of the binding modes and the ranking of the modes are dependent on the number and the ionization of the DOX. Moreover, it is found that the VDW interaction dominates in the TDN–unionized DOX interaction and the ELE interaction dominates in the TDN–ionized DOX interaction. Interestingly, the inside-corner mode is found to be the most energy-favorable mode in nearly all cases, which is due to the topological structure of the TDN. In general, the present study reveals the loading mechanism of the drug molecules into the TDN at the molecular level and may shed some light on the future development of functional DNA nanostructures in biomedical applications.

Author contributions

Conceptualization: H. M. D. and Y. Q. M.; methodology: H. M. D. and Y. X.; formal analysis: Y. X.; investigation and visualization: Y. X. and S. W. H.; writing – original draft: H. M. D. and Y. X.; writing – review and editing: H. M. D. and Y. X.; funding acquisition: H. M. D. and Y. Q. M.; project administration: H. M. D. All authors approved the final version of the manuscript.

Conflicts of interest

The authors declare no competing financial interest.