Here, we report the results of all-atom molecular dynamics simulations of RNA translocation through carbon nanotube membranes in explicit solvent. These simulations allow us to study membrane translocation at atomic detail, extending earlier studies of coarse-grained polymer translocation models (20–23). By including detailed descriptions of water and ions, the simulations capture electrostatic and hydrophobic solvation effects on the translocation processes and permit a detailed study of sequence dependences. As we will show, hydrophobic interactions of the bases with the nanotube pores can transiently trap RNA at the pore walls. To analyze the resulting multiphasic kinetics, we use a trap-diffusion model that combines diffusive RNA motion along the pore with random trapping. We will show that the model motivated by our simulations can also account for detailed experimental measurements of DNA translocation through α-hemolysin pores (5). By varying the strength of the transmembrane potential, we will study translocation of RNA at different driving forces. The simulations also give us access to the RNA structure, which will allow us to study the role of conformational dynamics in RNA translocation. Finally, we will discuss the effects of RNA sequence on the overall translocation rate and implications on molecular separation devices.

Simulations. Molecular dynamics simulations at a constant pressure of 1 bar and a temperature of T = 300 K were performed with the namd program (24), the assisted model building with energy refinement 94 (AMBER 94) force field (25) for the RNA, and the transferable intermolecular potential 3 point (TIP3P) model of water (26). Potential parameters for single-wall carbon nanotubes are based on sp² carbon parameters in the AMBER 94 force field (12, 25).

In our molecular dynamics simulations, single-stranded RNA molecules A₆ and U₆ with 6 bases of adenine and uracil, respectively, are driven across membranes of hexagonally packed single-wall carbon nanotubes by electric fields of up to 0.5 V/nm (Fig. 1). To form a nanotube membrane (13) under periodic boundary conditions, four (14, 14) “armchair”-type nanotubes (1.33-nm length and 1.87-nm carbon-to-carbon diameter with a pore size of ≈1.5 nm) were packed into a hexagonal lattice in the x–y plane, with their tube axes aligned with the z axis of the simulation cell. The pore diameter of the nanotubes is comparable to the constriction region of the protein channel α-hemolysin (7) and to the diameter reported for carbon nanotube bundles used as fluid-flow sensors (27).

Fig. 1. Schematic of simulations. (A) Translocation of negatively charged single-stranded RNA through carbon nanotube membranes is driven by the electric field E along the z axis normal to the membranes. Water molecules and counterions are not shown for clarity. (more ...)

Five initial configurations each of A₆ and U₆ RNA were taken from two 5-ns simulations in free solution, starting from single strands in A-form RNA duplexes (28). The RNA molecules, together with five neutralizing K⁺ counterions (29), were then inserted into the water phase of equilibrated nanotube/water interfacial systems. The resulting RNA/nanotube/water systems with 2,267 water molecules were further equilibrated for 1 ns without electric field, and then 100 ps with field. The simulation boxes were ≈4.4, 3.8, and 5.0 nm long in x, y, and z directions, respectively. We performed 30 RNA/nanotube/water simulations of 20–28 ns duration at electric fields of E₀ = 0.3, 0.4, and 0.5 V/nm (5 simulations each for A₆ and U₆ RNA) and 10 simulations without field covering between 11 and 26 ns.

Electric Fields. Particle-mesh Ewald summation (30) with conducting boundary conditions was used to calculate long-range electrostatic interactions. The electric field averaged over the simulation cell is then equal to the external electric field E₀ (31), but the local electric field in the bulk water phase and inside the nanotube pores deviates significantly from E₀. To estimate the local electric field across the water/nanotube interface, we calculated the electrostatic potential along the z direction from simulations of water/nanotube interfaces without RNA. The potential Φ(z) averaged parallel to the membrane is calculated from the charge density ρ_q(z) through Gauss' law, , with ε₀ the vacuum permittivity and L_z the length of the simulation cell in the z direction. The integration constant E′ is chosen to enforce the boundary condition Φ(0) - Φ(L_z) = E₀L_z for Ewald summation under conducting boundary conditions (31, 32). From the slopes of the electrostatic potentials, we estimate average electric fields of 0.15, 0.17, and 0.19 V/nm in the bulk phase, and 0.56, 0.83, and 1.09 V/nm inside the nanotube membrane at external fields of E₀ = 0.3, 0.4, and 0.5 V/nm, respectively. The electric field in the bulk water region calculated from the slope of Φ(z) and the orientation of water dipoles are consistent with the previously determined relationship between the polarization and the electric field of TIP3P bulk water (32).

The electric field in typical nucleic acid translocation experiments with membrane-inserted α-hemolysin channels is in the range of 0.02–0.04 V/nm (120 mV across 3- to 5-nm membrane thickness) (2, 5), which is ≈20 times smaller than the smallest estimated electric field in our nanotube membrane simulations. Use of relatively high electric fields was necessary to induce at least nanometer-scale motion on the nanosecond simulation time scales. We recently (32) showed that the hydrodynamic and electrophoretic properties (33) obtained from simulations of single-stranded RNA homopolymers in free solution at similar electric fields are in accord with experiment. Also, in a recent microsecond-electrophoresis experiment with micrometer-separation paths, electric fields exceeding 0.01 V/nm were achieved without any significant alteration of mobilities because of Joule heating (34). Solid-state nanopore systems like the nanotubes used in our simulations are expected to be more stable than biological membranes (9, 35). Indeed, in all our simulations under high electric fields, the carbon nanotube membrane models were found to remain intact, held together by strong hydrophobic interactions between neighboring nanotubes, even though the nanotubes were flexible and could move relative to each other and the solvent (12, 13).

RNA Translocation Kinetics. In the analysis of RNA translocation kinetics, the duration of entry events is defined as the time between the first RNA atom entering into the membrane region (within 1.1 nm of the membrane center), and either the last RNA atom entering or the first atom exiting the membrane region. The duration of exit events is defined as the time between the completion of the entry event and the last atom leaving the membrane region. To correct for incomplete entry or exit events (i.e., the simulation was terminated before the RNA molecule finished entering or exiting the membrane region, respectively), we calculated average entry and exit times by using a maximum-likelihood estimate for randomly censored data, , where N is the total number of attempts, n is the number of successful events, and t_i is the duration of an event, if completed, and otherwise the time between the start of an event and the end of the simulation. For cumulative distributions, we used the Kaplan–Meier nonparametric estimator (36).

Trap-Diffusion Model. To describe the exit-time statistics, we consider a model in which RNA molecules can become transiently trapped as they pass through the pore. Similar models have been used previously in modeling molecular motors (37) and polymer translocation (21). The exit-time probability distribution can be obtained by solving the coupled equations describing the exchange between diffusively drifting and bound RNA populations,

[1]

[2]

In the trap-diffusion model, the exit time is given by the arrival time of molecules that start at a fixed distance Δz (here, 2.2 nm) upstream from the exit. We obtain the Laplace transform of the exit-time distribution by calculating the outgoing flux from the Green's functions for Eqs. 1 and 2 on a semiinfinite interval with absorbing boundary conditions at the exit:

[3]

[4]

[5]

[6]

[7]

[8]

[9]

RNA Translocation. During the combined simulation time of over 800 ns, we observe 173 translocations of RNA and 14,435 translocations of K⁺ counterions across the nanotube membranes. RNA transport across the membranes is driven by electric fields acting on the charged phosphate groups, but hydrophobic interactions between the bases of the RNA and the carbon nanotubes strongly affect the translocation process, as shown in Figs. 1B and 2 A–C.

Fig. 2. Representative configurations of RNA during translocation through the nanotube membrane. (A) Head-on entry. (B) Translocation. (C) Exit with single trapped base. The close-up shows a tight hydrophobic interaction between an unstacked base and the nanotube (more ...)

Fig. 2 shows snapshots of typical translocation events. Rapid entry occurs if a terminal phosphate group enters first into a nanotube pore (Fig. 2 A). Less frequently, one of the central phosphate groups is pulled into a pore by the electric field (Fig. 2D). The resulting confinement distorts the RNA structure, forcing it into a tightly packed U-shape. Nevertheless, despite the severe geometric restraints, we find that the small RNA molecules studied here translocate through the membrane, driven by the high electric fields inside the pore, without disrupting the nanotube membrane.

To determine the orientation of RNA before, during, and after translocations, we calculate the cosine of the angle θ between the long axis of RNA, as determined from the inertia tensor, and the z axis as a function of the z position of RNA (Fig. 3A). A decrease of cosθ near z ≈ 1.1 nm indicates that RNA molecules initially dock onto the membrane surface and, typically, spread out on it under the combined influence of hydrophobic and electrostatic forces. A large increase of cosθ around z = 0 indicates that RNAs are oriented along the z axis during the translocation events. The drop in cosθ just outside the membrane (z ≈ -1.0 nm) is caused by RNA molecules that remain hydrophobically attached to the membrane upon exit and spread out on it. The subsequent increase of cosθ for z < -1.0 nm is caused in part by the elongation of RNA during the exit process if one end remains in hydrophobic contact with the nanotube while the other end is pulled away from the nanotube by the electric field (Fig. 2C). The orientation of the RNA in the presence of electric fields is also evident in the distributions of base and phosphate groups along the z direction normal to the membrane plane (Fig. 3B). On the entrance side (z > 0), phosphate groups penetrate into the nanotube membranes, pulled by the high electric fields inside the pores, whereas base atoms are built up along the water/nanotube interface. On the exit side (z < 0), we find a high density of base atoms inside the carbon nanotubes, originating from bases forming long-lived hydrophobic interactions with the nanotube walls. In contrast, phosphate groups concentrate outside the nanotube membrane near the exit. The results shown in Fig. 3 are for an average over A₆ RNA trajectories, but the main features are present also in data from individual runs. Similar results are obtained for U₆ (data not shown).

Fig. 3. Orientational and density profiles of RNA. (A) Cosine of the angle θ between the long axis of A6 RNA, as determined from the inertia tensor, and the z axis as a function of RNA position, averaged over all simulations with applied electric fields. (more ...)

Kinetics of RNA Exit from Pores. During RNA exit from the membrane, the charged phosphate groups usually proceed rapidly into the high-dielectric, low-field solvent, as shown in Fig. 3. However, attractive hydrophobic interactions between individual bases and the nanotube channel can delay the release of RNA into the solvent (Fig. 2C). We observe structures with up to five RNA bases simultaneously attached to the nanotube pore wall (Fig. 2E). At a field of 0.3 V/nm, hydrophobic trapping prevents the release of RNA from the nanotube in a significant fraction of the ≈20-ns simulations. In the absence of an external electric field driving translocation, RNA slides into the nanotube pores but remains trapped by strong hydrophobic interactions between the bases and the channel walls. In some cases, RNA can become trapped simultaneously in two neighboring pores. Such trapping results in tight hydrophobic contacts between the bases and the nanotube channels but leaves much of the charged sugarphosphate backbone fully exposed to the solvent outside the membrane (Fig. 2F).

Results for the translocation kinetics of RNA through carbon nanotube membranes are shown in Figs. 4 and 5. The exit-time distribution of A₆ RNA at 0.5 V/nm (Fig. 4A) has a lag time, followed by a sharp peak of fast exit events at t ≈ 150 ps, and then a long tail extending to ≈3 ns. This distribution of exit times is described well by the trap-diffusion model of Eqs. 1 and 2, which combines free RNA diffusion with intermittent random trapping inside the pore. In Fig. 4A, we have included results calculated with a drift distance of Δz = 2.2 nm across the membrane, an RNA diffusion coefficient of D = 0.7 × 10^-5 cm²·s^-1, and a drift velocity v = 14 nm·ns^-1 for untrapped RNA at the highest field. These values are similar to those determined previously for free diffusion (32). The estimated rate coefficients for attaching to and dissociating from the pore wall are k_on = 7 ns^-1 and k_off = 2.5 ns^-1, respectively. With those values, the roughly 30% of RNA molecules that exit the tube without becoming trapped produce the sharp peak in the exit times. Subsequently, arriving RNA molecules were trapped in the pore and result in the slow tail found in the distribution. The amplitude and decay of the exit-time tail are accurately reproduced by the approximate analytic expression for slow diffusion, Eq. 9.

Fig. 4. Kinetics of RNA exit from the membrane pores. (A) Cumulative probability distributions Pexit(t) of exit times t of A6 RNA (blue line) at an electric field of 0.5 V/nm compared with the trap-diffusion model (solid black line). (Inset) The corresponding (more ...)

Fig. 5. Kinetics of RNA entry into the membrane pores. (A) Cumulative probability distributions Pentry(t) for the duration t of entry events of A6 (blue; right-hand scale) and U6 (red; left-hand scale) RNA at electric fields of 0.3 V/nm (long dashed lines), 0.4 (more ...)

As expected, exit slows down with decreasing electric field. This slow down is the result of a reduced drift velocity and, more importantly, a slower dissociation of the hydrophobically trapped RNA molecules from the pore walls. With some uncertainty, caused by the small number of events at low fields, we find that the characteristic time of the slow decay in the exit-time distributions seems to decrease exponentially with the applied electric field. This result is consistent with the experimental observations of Bates et al. (5). Based on Eq. 8, for slow diffusion, the characteristic time for exit is dominated by the rate of dissociation from the pore wall, k_off. An exponential field dependence of k_off is, indeed, expected from phenomenological theories of force-induced ligand dissociation (39, 40).

Exit rates (Fig. 4B) depend more strongly on electric fields than rates of entry into the pores (Fig. 5B), and drop to <1/(10 ns) at a field of E₀ = 0.3 V/nm, below which exit is limiting the rate of translocation. This sharp reduction in exit rates with decreasing electric fields (Fig. 4B) explains why we have not observed any exit of RNA from the nanotubes without external electric field on the simulation time scale. Indeed, if we assume an exponential dependence of the exit times on the electric field inside the pores, we estimate exit times of roughly 100 μs to 1 ms for a 120-mV voltage drop across the membrane. This result is on a time scale similar to that measured for DNA translocation through α-hemolysin (5).

Strong hydrophobic interactions between RNA bases and the nonpolar nanotubes (Fig. 2C) trap RNA inside the membrane pores, despite the cost of partially desolvating the phosphatebackbone charges and the entropic cost of confining a flexible polymer. At all electric fields, we find longer average exit times of A₆ compared with U₆ RNA (Fig. 4B). Based on their respective exit-time distributions, the difference in exit rates seems to be caused by a slightly higher probability of A₆ to become trapped, which we attribute to stronger hydrophobic interactions of the larger A₆ purine bases with the nanotubes, compared with the U₆ pyrimidine bases.

Kinetics of RNA Entry into Pores. Fig. 5A shows the cumulative distributions of entry times of RNA into the membrane pores for different electric fields and RNA sequences. The entry-time distributions show two distinct phases. A fast exponential phase is independent of the external electric field but depends on the RNA sequence, with time constants of ≈1 ns for U₆ and 1.3 ns for A₆ RNA, respectively. A second slow phase, although negligible at 0.5 V/nm, has an amplitude that grows strongly with decreasing electric field. We associate this slow phase with RNA molecules that became transiently trapped upon entry because of unfavorable geometry (e.g., sideways entry; Fig. 2D). To explain the fast, field-independent phase, we point out that, in addition to electrostatic and hydrophobic effects, the conformational dynamics of RNA molecules will affect their rate of translocation through the carbon nanotube membranes. Conformational dynamics should be particularly relevant during RNA entry into the pores. To explore the hypothesis that the initial fast phase in the entry-time distributions of Fig. 5A is caused by conformational dynamics on a nanosecond timescale, we measure the rate at which RNA conformations evolve in time. We quantify the RNA structural change by calculating a time-dependent rms distance (rmsd). In Fig. 7, which is published as supporting information on the PNAS web site, we plot the average rmsd between two RNA conformations separated by fixed time differences Δt along trajectories in free solution (32) as a function of Δt. A simple linear scaling of time by 1.3 and amplitude by 0.9 brings the rmsd-vs.-time curve of U₆ RNA to overlap with A₆. The ≈30% faster configurational dynamics of U₆ compared with A₆ RNA is consistent with the 30% difference in the time constants of the fast phase for RNA entry into the pores observed in Fig. 5A. This result suggests that nanosecond configurational dynamics is indeed a dominant factor in the fast phase of RNA entry that becomes rate limiting at high electric fields and accounts for shorter average entry times of U₆ RNA into the nanotube membrane (Fig. 5B). Although the entry rates in Fig. 5B are almost identical at electric fields of 0.4 and 0.5 V/nm, trapping of RNA at the membrane interface becomes relevant at a lower electric field of 0.3 V/nm, resulting in a slow tail in the entry-time distributions (Fig. 5A) and a sharp increase in the average duration of entry events.

Concluding Remarks. Exit-time distributions similar to those obtained here in the simulations and predicted by the trap-diffusion model of Eqs. 1 and 2 were observed in recent single-molecule measurements of single-stranded DNA exiting from membrane-inserted α-hemolysin channels (5). The approximately exponential tail after a rapid initial burst in the measured exit times was attributed to pore-binding DNA populations, consistent with our observations and model. In Fig. 6, we show that the trap-diffusion model suggested by our simulations can indeed account for the experimental translocation data. The model accurately captures the multiphasic kinetics measured for DNA pore blockage both with and without applied electric field. In particular, we recover the slow “tails” in the translocation and exit-time distributions without invoking more elaborate models (38).

Fig. 6. DNA translocation experiment of Bates et al. (5) modeled by the trap-diffusion model. (A) Probability distributions of poly(dA)60 translocation time t across α-hemolysin channels at 120 mV. The experimental results (bars and open circles for probability (more ...)

In our molecular dynamics simulations, we found that the translocation kinetics of RNA through the nanotube membranes is sequence-dependent. At first glance, this finding may seem surprising, considering that the nanotube pores are atomically smooth and featureless, compared with more structured protein pores (2, 4–6). However, we showed that the RNA molecules differ in their affinity for the pore wall and their conformational dynamics. These differences affect the rates of exit and entry, respectively, and thus the overall rate of RNA translocation, and highlight that RNA is not just a simple polyanion. Our results suggest that sequence information can be extracted from synthetic-pore translocation data (35, 41, 42) and that simulations could prove useful in the experimental analysis.

The number of U₆ translocations through the nanotube membranes is significantly larger than that of A₆ RNA (52, 42, and 7 vs. 45, 25, and 2 at electric fields of 0.5, 0.4, and 0.3 V/nm, respectively, during simulations covering ≈110 ns for each RNA and field). By averaging over multiple trajectories of individual RNA molecules through the periodic nanotube membranes, we effectively recreate an ensemble measurement from our single-molecule simulation data. Differences in entry and exit rates result in sequence-dependent electrophoretic mobilities of RNA, and, as shown in Fig. 1B, the “ensemble” of U₆ RNA gradually separates from that of A₆. On a time scale of multiple translocations, the simulation system of periodically spaced nanotube membranes thus acts similarly to a porous gel in an electrophoretic measurement. We expect the underlying differences in RNA conformational flexibility and affinity for the pore walls to result in sequence-dependent rates of membrane translocation also for longer polymers in larger pores (10, 11), and in nanometer-size pores of other materials (9, 35, 41), as long as frequent biopolymer-pore interactions occur during translocation. Functionalization of the carbon nanotube rims, for instance by attaching complementary nucleic acids (43) for sequence-dependent trapping (6), should further increase the specificity of carbon nanotube membranes in separation devices.

We thank Drs. Attila Szabo, Alexander Berezhkovskii, and Christine Peter for many helpful discussions and suggestions. This study used the biowulf pc/linux cluster at the National Institutes of Health.