Gold nanoparticle design for RNA compaction

Double-stranded RNAs (dsRNAs) play an essential, though unheralded, role in biology and medicine by regulating gene expression and activating the innate immune systems of eukaryotes.^1,2 Furthermore, advanced therapeutic methods involving RNAs hold a great promise in treating various gene-related diseases and disorders.^3,4 Short dsRNAs, microRNAs (miRNAs), and short interfering RNAs (siRNAs) regulate the expression of genes through the RNA-induced silencing complex (RISC) by destroying messenger RNA with a complementary sequence or by inhibiting translation.^4–6 This allows for targeted, transient knockdown of a specific gene, without the need for genome modification or nuclear transport. However, the development of RNA packaging is a nontrivial task due to its inherent instability, high negative charge, and significant size.^3,7

Delivery of dsRNA is typically implemented by creating carrier vectors that may be of various viral or nonviral forms. Choosing the optimal vector type and its configuration becomes a nontrivial task.⁷ For example, dsRNA viruses, such as bacteriophages φ6 and φ12, package single-stranded RNA (ssRNA) using an ATPase motor protein.^8,9 A complementary strand to the ssRNA is then created to form dsRNA within the viral capsid. Interestingly, dsRNA viral capsids seem to contain lower concentrations of polyamines than dsDNA viruses,^10–13 suggesting that the polyamine presence is not so strong a driving force for dsRNA compaction. The absence of viruses that package dsRNA in its duplex form brings to our attention the difficulty of bending and wrapping dsRNA as compared to ssRNA. Nevertheless, the final state of the packaged dsRNA shows that bending and wrapping dsRNA are possible, given an accessible folding pathway. The dsRNA is significantly stiffer than dsDNA^14–16 and is known to resist condensation in the presence of complexes of multivalent cations,¹⁷ such as cobalt-hexamine [Co-hex, Co(NH₃)₆³⁺] or polycations¹⁸ typically found at physiological pH, such as spermine. The deep binding of Co-hex in dsRNA's major groove does not reverse the negative surface charge and, therefore, does not cause condensation.¹⁹ However, dsRNAs delivered by cationic liposomes have been employed in gene therapy with recent clinical success.²⁰ 

Ligand-functionalized metal nanoparticles have been actively investigated as nonviral vectors for RNA delivery.^21–23 The most researched nanomaterial is a gold nanoparticle (AuNP),^24,25 which possesses a wide range of advantages such as biocompatibility, visibility in the infrared spectrum, capability of photothermal activation, and ability to be functionalized with various surface chemistries.^26,27 Moreover, numerous computational studies of AuNPs capped with various ligands have been reported.^28–37 In recent studies, AuNPs have been shown to deliver dsRNA to nematodes,²⁵ pests, and insects^23,24 or to transport siRNAs to tumor cancer cells;⁴ however, the NP electrochemical properties and type of surface functionalization were shown to play an important role in the delivery rate and overall effect of tissue treatment.⁴ Thus, deeper understanding of the required properties of AuNP is necessary to improve gene delivery methods, which are, therefore, needed to realize the potential of RNA therapies.

Conformation of NP-nucleic acid complexes depends on many factors: nucleic acid type, salt concentration, nanoparticle geometry, ligand chemistry, etc. All-atom molecular modeling techniques provide the opportunity to study such systems in detail.³⁷ Our previous research showed that while small, 1.5 nm, AuNPs capped with SC₁₁NH₃⁺ ligands wrapped dsDNA with the same efficiency as the nucleosome at physiological salt concentrations, the same nanoparticles fail to induce dsRNA bending.³⁸ Moreover, helical distortions in dsRNA were observed at high AuNP charge and low salt concentrations. Here, we use all-atom molecular dynamics (MD) simulations to investigate the essential properties of AuNPs, such as nanoparticle core size, ligand length, and hydrophobicity [Figs. 1(a)–1(c), Table S1],⁷³ that can facilitate the compaction of 100 base pair dsRNA. The first nanoparticle (NP1) has the core diameter of 1.5 nm with 60 S-(CH₂)₁₁-NH₃⁺ ligands attached to the surface, as in previous studies.^39–41 The second nanoparticle (NP2) has an alkyl chain ligand shorter by one carbon and a bulkier end group [S-(CH₂)₁₀-N(CH₃)₃⁺], which has the effect of increasing the hydrophobicity of the end group and reducing packing of the ligands. Third nanoparticle (NP3) has a 3.0 nm diameter of the gold core with 200 shorter ligands S-(CH₂)₆-NH₃⁺.⁴² The AuNP designs were chosen to maintain the same overall size and charge [Figs. 1(a)–1(c)]. Our simulations showed the ability of the ligand to explore the environment surrounding the nanoparticle to be the principal factor determining whether a nanoparticle could cause conformational changes to RNA.

Systems consisted of 100bp RNA with sequence AUCAAUAUCCACCUGCAGAUUCUACCAAAAGUGUAUUUGGAAACUGCUCCAUCAAAAGGCAUGUUCAGCUGAAUUCAGCUGAACAUGCCUUUUGAUGGAG and one ligand functionalized AuNP in an explicit solvent box. This sequence is analogous to the first 100 base-pairs of DNA in a crystal structure of the histone octamer.⁴³ RNA was aligned along the z-axis of the simulation box. Because of the length of RNA, the solvent box was much larger in x and y dimensions, with a solvent buffer of 20–25 Å to allow rotation of the RNA. Nanoparticles were initially placed with a center of mass distance of at least 4 nm from RNA's helical axis. All systems were solvated with TIP3P water, neutralizing ions, and 0.1M NaCl salt concentration. The RNA was described by the ff99bsc0 force field.⁴⁴ Systems typically contained between 750 000 to just over 1 × 10⁶ atoms (long RNA systems) including solvent and ions. The force field for nanoparticles were the same as described in our earlier work,^39,40,45 where the nanoparticle consists of a spherical gold core functionalized with alkanethiol ligands. In this paper, we use gold nanoparticles with a 3.0 nm diameter with 200 total ligands or a 1.5 nm diameter with 60 total ligands. The number of ligands on the 1.5 nm particle is based on experimentally determined measurements showing that a NP with a size of 1.5 nm has 60 alkanethiol ligands on its surface.⁴¹ For the larger nanoparticle, we calculated the number of ligands based on the theoretical estimate of the number of ligands that can be attached to curved NP surfaces.⁴² Partial charges on nanoparticle ligands were calculated using the R.E.D Server.⁴⁶ For the 1.5 nm core NP, all ligands had a charge of +1. The ligands were S-(CH₂)₁₁-R where R indicates either an amine end group (R-NH3⁺) or a trimethyl amine end group [R-N(CH₃)₃⁺]. For the 3.0 nm core NPs, the ligands were S-(CH₂)₅-R, where R is either an ammonium end group (R-NH₃⁺) or a methyl end group (R-CH₃). Additional nanoparticle-only control simulations consisted of a single nanoparticle with a 15 Å buffer of TIP3P water, chloride ions to neutralize the nanoparticle's charge, and additional eight sodium and eight chloride ions. To confirm that 1.5 nm NPs have no tendency to cause compaction of dsRNA, we evaluated time-dependence of the radius of gyration [Fig. S5(a)],⁷³ end-to-end distance [Fig. S5(b)],⁷³ and bending angle as a function of sequence [Figs. S6(a) and S6(b)]⁷³ for dsRNA coupled with 1.5 and 3.0 nm NPs.

All computational procedures were performed using Gaussian 09 software.⁴⁷ We employed the hybrid functional of Truhlar and Zhao M06-2X⁴⁸ including some amount of dispersion correction along with the Los Alamos effective core potential with the associated basis set for all atoms, Lanl2dz.⁴⁹ Geometry optimizations (without any symmetry restrictions) and frequency calculations were performed with implicit effects of water taken into account (dielectric constant ɛ = 78.3553) using the self-consistent reaction field IEF-PCM method.^50,51 The UFF default model used in the Gaussian 09 package, with the electrostatic scaling factor α set to 1.0, was employed. The analysis of charge distribution was performed using the Natural Bond Orbital (NBO) method as implemented in Gaussian 9 program,⁵² using the M06-2X/Lanl2dz approach with implicit water effects. Frontier molecular orbitals (FMOs) also were calculated at the M06-2X/Lanl2dz level with implicit water effects. Open GL version of Molden 5.8.2 visualization program was used for the visualization of structures and FMOs of the studied NPs,⁵³ and Avogadro, version 1.1.1, was used to visualize the molecular orbitals and molecular electrostatic potential (MEP).⁵⁴ 

All-atom molecular dynamics simulations were performed using Amber14⁵⁵ and Amber 16⁵⁶ software packages. Minimization and solvent equilibration simulations were run using CPUs, while production simulations were run using Amber's GPU-accelerated code.^57–59 Before production simulations, we used an 11-step minimization, heating, and equilibration modified routine protocol similar to that described in our prior work.⁶⁰ Briefly, systems were first minimized for 10 000 cycles using a steepest descent method, followed by several steps of constrained heating and solvent equilibration. During equilibration, the RNA and nanoparticles were held in place using gradually reduced harmonic constraints using an isothermal-isobaric ensemble to allow adjustment of system density. After solvent equilibration, all constraints were released, and simulations were switched to a constant volume and temperature (NVT) ensemble using a Berendsen thermostat.⁶¹ Simulations were run with a 2 fs timestep. Solvent molecules were constrained using the SHAKE algorithm.⁶² Particle mesh Ewald summation was used to calculate long-range interaction. Simulations continued until structural conversion occurred but for a minimum of 120 ns for all NP-RNA systems.

Simulations were analyzed using CPPTRAJ⁶³ and in-house scripts. For the analysis of RNA base pair parameters, we used the program Curves+.⁶⁴ Visualizations were performed using the Visual Molecular Dynamics (VMD) program.⁶⁵ 

The relative anisotropy index, κ², which ranges from 0 (spherical system) to 1 (line), can quantify nanoparticle ligand shell shape. Representative snapshots from simulations of the three nanoparticle designs in 0.1M NaCl solution are displayed in Fig. 1(d). Due to longer ligands, the ligand shells for NP1 and NP2 are more flexible and show greater shape deviation than the ligand shell of NP3. Measurement of the relative shape anisotropy, κ², from simulations [Fig. 1(e), Table S2]⁷³ of the three nanoparticle designs shows that although all nanoparticle ligand shells are somewhat spherical with approximately the same radius of gyration (Fig. S1),⁷³ the anisotropy is reduced by bulky end groups (NP2) or by shortening ligand length (NP3).

To characterize ligand behavior, we calculated the ligand excess free volume defined in Eq. (1), which gives the excess volume a ligand can sample due to nanoparticle curvature relative to a flat surface,⁶⁶ 

where σ specifies the ligand grafting density, l is the ligand length, and r is the NP radius. As can be seen from Eq. (1), the ligand excess free volume increases with increasing ligand length and decreasing nanoparticle diameter (Table S2).⁷³ Figure 1(f) shows that NP1 and NP2 have relatively high ligand excess free volumes, while NP3 has a significantly smaller excess free volume. Decreasing ligand length and increasing the size of the NP core reduce ligand excess free volume through both r and l parameters. Simulations of AuNPs in TIP3P-water and 0.025M NaCl showed that NP1 end groups have high mobility compared to NP2 and NP3 (Fig. S2).⁷³ Thus, with these nanoparticle designs, the effect of nanoparticle ligand shell sphericity and ligand excess free volume on NP-RNA binding varied independently of overall nanoparticle charge and size.

To accurately assess the behavior of the ligands on the surface of the nanoparticle, we performed density functional theory (DFT) simulations (for details, see Figs. S5–S9 and Tables S4 and S5).⁷³ We used a cluster of 20 gold atoms (Au₂₀ NP) as a model for AuNPs and three types of ligands (Fig. S9)⁷³ reproducing the NP designs in Fig. 1. Our results show that the behavior of the ligands depends on the ligand length and type. The longer-chain ligands S(CH₂)₁₁-NH₃⁺ tend to keep closer to each other as it can be seen in Fig. S9(a).⁷³ This happens due to attractive dispersion interactions between the ligand chains that offset the repulsive Coulombic interactions between the NH₃⁺ groups of the ligands. The situation is somewhat different for another type of the longer-chain ligands, S(CH₂)₁₀-N(CH₃)₃⁺: these ligands tend to stay somewhat further apart from each other due to their bulkier head groups [Fig. S9(b)].⁷³ However, still some similarity in the alignment might be noticed between both longer-chain ligands [Figs. S9(a) and S9(b)].⁷³ The ligand excess free volume is slightly smaller for S(CH₂)₁₀-N(CH₃)₃⁺ ligands, which is in line with the DFT results. However, for shorter ligands, S(CH₂)₆-NH₃⁺, the situation is noticeably different: not only their S-centers move away from each other on the AuNP surface, but the ligand chains move away from each other due to Coulombic repulsions between the NH₃⁺ groups being stronger than the attractive dispersion interactions between the ligand chains [Figs. S9(c) and S9(d)].⁷³ These observations are in a qualitative agreement with MD simulation results, indicating that the shorter ligands experience extra mobility due to increased repulsion.

We then performed simulations of AuNPs’ binding and interactions with the 100-bp dsRNA in 0.1M NaCl and explicit water. Representative snapshots of final system conformations from our simulations are shown in Figs. 2(a)–2(c). We observed that while the charge and size of each NP are about the same, the final RNA conformation differs depending on NP properties. NP1 and NP2 cannot compact the dsRNA; however, RNA bends and wraps around NP3. Measurement of the RNA radius of gyration in each of the three simulations [Fig. 2(d)] reveals that NP1 and NP2 do not cause dsRNA conformational changes (RNA remains linear), but NP3 causes a reduction in the radius of gyration through bending of the RNA around the nanoparticle.

Our simulations show that while the overall NP size and charge are the same, the changes in NP designs result in differing RNA final conformation. Simple theory models, however, would not predict this difference. For instance, in the wormlike chain model, atomic details are neglected, and the parameters determining bending energy are bending angle, segment length, and persistence length.^67,68 In previously published theoretical calculations, parameters such as the Bjerrum length and charge separation along the molecule were used to calculate that 6% minimum charge neutralization of DNA, which is required for bending on a superhelical ramp similar to the nucleosome.⁶⁹ Based solely on these types of descriptors (persistence length, linear charge density, etc.), the theory would not predict any difference in the interaction of nucleic acids with three nanoparticles in this study. Thus, the subtle details of the molecular structure of the nanoparticles or nucleic acid molecules play a significant role in their interaction and behavior.

We next examine the atomic details of ligand-RNA binding. The snapshots in Fig. 2 show differences in end-group binding locations; the charged groups of NP1 bind exclusively to RNA's major groove [Fig. 2(a)], while NP2 and NP3 make contact with both grooves [Figs. 2(b) and 2(c)]. Measurements of the preferred binding locations of NP ligand end groups, shown in Fig. 3, demonstrate this quantitatively. Figure 3 shows the binding location of the nitrogen on the ligand end group for each nanoparticle as a radial distribution function (a)–(c) along with a heatmap looking down the z-axis of the RNA helix [(e)–(g) for ligands attached to NPs and (i)–(l) for ligands not attached to NPs). Figure 3(h) shows a schematic of major and minor groove locations on an RNA base pair. In radial distribution plots, the red line at 0 represents the RNA's helical axis, while the black line indicates the outermost positions of the RNA's helix.

We also performed control molecular dynamics simulations of ligand molecules, not attached to a nanoparticle surface, with shorter lengths (30 bp) of dsRNA. When ligands are free and not constrained in their movement by the nanoparticle, as in control simulations (Table S1),⁷³ the binding distribution of ligands is dependent on ligand length and hydrophobicity [Figs. 3(i)–3(l)]. Packing of the alkyl backbone prevents longer ligands from binding deep within major groove in contrast to shorter ligands. The hydrophobicity of RNA's minor groove leads to an increase in end-group binding for ligands with trimethylated end groups [N(CH₃)₃⁺] and the short alkanethiol ligands. One of the possible conclusions based on these data may suggest that an increase in ligand excess free volume allows for a more energetically favorable ligand distribution by alteration of the shape of the ligand corona, as opposed to a conformational change in the dsRNA.

Binding of NP1 to RNA is consistent with our earlier simulations, where the NH₃⁺ end groups bound inside the nucleic acid helix.³⁹ For both smaller core NP1 and NP2, the ligands bind primarily inside RNA, as shown in radial distribution plots. However, for NP3 with a smaller ligand excess free volume, the movement of the ligands was restricted (Fig. S2)⁷³ and the majority of charged ligands bound to the outside of RNA. This binding to the outside of the RNA helix was associated with RNA bending. In RNA-NP systems where ligands bound inside the major groove, no RNA bending occurred.

Interestingly, from the radial distribution plot, NP3 still appears to have ligands bound inside the RNA's helix [as it can be seen in Fig. 3(c)]. Examining the location of ligand binding in major and minor [Figs. 2(g), 2(h), 3(g), 3(k) and 3(l)] grooves provides more insight into this observation. Short unbound ligands travel deeper into RNA major grooves, residing closer to the RNA's helical axis. Additionally, bound ligands generally show a greater preference for the minor groove compared to free ligands.

For NP3, which is the only nanoparticle that induced RNA bending, charged ligand groups do not go into RNA's major groove, as can be seen by lack of color on the heatmap on Fig. 3(g). In contrast, for NP1 and NP2 that did not induce dsRNA bending [Figs. 3(e) and 3(f)], the ligands bind inside the RNA's major groove. For NP2 and NP3, which are the more spherical nanoparticles, some ligands also contact RNA's minor groove. Also, for NP2 and NP3, where ligands bind inside the minor groove, this may be due also to increased hydrophobicity of the NP end-groups, as the RNA minor groove is more hydrophobic. Overall, RNA bending with nanoparticles is associated with the absence of ligands inside the major groove of RNA.

We next examine atomic details of dsRNA. The orientation of two nucleic acid base pairs relative to one another may be described using six parameters: three are translational (shift, slide, and rise) and three are rotational (roll, tilt, and twist). The roll base pair parameter, which describes groove width and depth, has been shown to be particularly important for DNA bending.^39,70 Bending of nucleic acids can be achieved through periodic variations in the roll and tilt parameters or a combination of both. Figure S3 and Table S3⁷³ show models of RNA bending that have been achieved through the application of periodic variation of roll and/or tilt. We built these models using the program 3DNA^71,72 and performed an energy minimization in implicit solvent in order to obtain energies associated with each type of bending for RNA and DNA. As shown in Table S3,⁷³ the bending energies for RNA and DNA do not differ significantly. The energy penalty per base-pair is the highest for changes in the tilt base-pair parameter (4.3 and 5.2 kcal/mol for RNA and DNA, respectively) and lower for roll or combination roll-tilt bending (∼2.1–2.8 kcal/mol for all samples).

Figure 4 shows axial bend, major groove width, and roll and tilt base pair parameters for dsRNA measured from each RNA-NP system simulation. For each parameter, we did a three-point smoothing average across the RNA helix in each snapshot to reduce noise from base pair movement before taking a time-average for plot values. Figure 4(a) shows that NP1 and NP2 do not cause significant bending. However, bending is higher for NP3. This bending is accompanied by a periodic variation on the major groove width [Fig. 4(b)]. Tilt showed no significant variation for any of the samples, indicating that it does not contribute greatly to RNA bending. However, Fig. 4(c) shows that bending of RNA with NP3 was associated with a periodic variation on RNA's roll parameter (bottom panel). This periodic variation was not observed for RNA with NP1 and NP2, indicating that RNA bending occurs through periodic changes in the roll parameter. Variations in roll parameter show direct correlation with RNA's major groove width. The width of the major groove is smallest for high roll angles. At lowest values, the major groove width is close to 0 Å, suggesting that RNA bending by nanoparticles can only occur when ligands are not inside the major groove. In contrast, the minor groove shows much smaller variation in width (Fig. S4).⁷³ 

Herein, we have examined the ability of three designs of ligand functionalized nanoparticles to cause dsRNA bending and compaction using large-scale all-atom molecular dynamic simulations. Results show that the effective nanoparticle should have a low ligand free volume that effectively inhibited cationic ligands from binding inside the major groove and caused bending and wrapping of RNA around the NP. NPs with long cationic ligands bind to RNA but did not cause bending. Characterization of RNA's base-pair parameters showed that RNA bending was associated with periodic variations in RNA's roll parameters and, consequently, major groove width. The minor groove was affected to a lesser extent. DFT studies on the model AuNP-ligand systems supported the simulations and characterization of NPs used in this study. Overall, results suggest the design rules for ligand functionalized metal NP design that can cause compaction of ds RNA.

This study was supported by the National Science Foundation (NSF) (Nos. CBET-1403871, DMR-1121107, CMMI-1150682, DGE-0946818, and DMR-2203979). J.A.N. was supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE-0946818. A.E.K. acknowledges the supercomputer resources of the Department of Chemistry, Institute of Technology of Aeronautics (ITA), São Jose dos Carlos, Brazil.

The authors have no conflicts to disclose.

Ethics approval is not required.

J.A.N. and Y.G.Y. designed research; J.A.N. performed research; J.A.N. and M.D.M. analyzed data; A.E.K. performed DFT studies; and J.A.N., M.D.M., A.V.G., and Y.G.Y. wrote the paper.

Jessica A. Nash: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (lead); Methodology (equal); Validation (equal); Visualization (equal); Writing – original draft (equal). Matthew D. Manning: Data curation (equal); Formal analysis (equal); Methodology (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Alexey V. Gulyuk: Data curation (lead); Methodology (equal); Software (equal); Writing – review & editing (equal). Aleksey E. Kuznetsov: Formal analysis (equal); Methodology (equal); Software (equal); Visualization (equal); Writing – review & editing (equal). Yaroslava G. Yingling: Conceptualization (lead); Data curation (equal); Funding acquisition (lead); Methodology (equal); Project administration (lead); Resources (lead); Supervision (lead); Writing – review & editing (lead).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

AuNP

 = gold nanoparticle

bp

 = base pair

DFT

 = density functional theory

GPU

 = graphics processing unit

MD

 = molecular dynamics

NP

 = nanoparticle

ns

 = nanosecond