PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Faraday Discuss. Author manuscript; available in PMC 2010 July 20.
Published in final edited form as:
Faraday Discuss. 2009; 143: 47–93.
PMCID: PMC2907245
NIHMSID: NIHMS217418

Molecular control of ionic conduction in polymer nanopores

Abstract

Polymeric nanopores show unique transport properties and have attracted a great deal of scientific interest as a test system to study ionic and molecular transport at the nanoscale. By means of all-atom molecular dynamics, we simulated the ion dynamics inside polymeric polyethylene terephthalate nanopores. For this purpose, we established a protocol to assemble atomic models of polymeric material into which we sculpted a nanopore model with the key features of experimental devices, namely a conical geometry and a negative surface charge density. Molecular dynamics simulations of ion currents through the pore show that the protonation state of the carboxyl group of exposed polyethylene terephthalate residues have a considerable effect on ion selectivity, by affecting ionic densities and electrostatic potentials inside the nanopores. The role of high concentrations of Ca2+ ions on charge inversion effects, earlier reported in experiments, were investigated in detail.

Keywords: polymer nanopores, polyethylene terephthalate, ion current, charge inversion

1 Introduction

Solid-state nanopores are small pores with radii of a few nanometers that are built in synthetic materials, such as silicon nitride1,2, silicon oxide3, silicon4, and polymer films58. An atomic model of a polymer nanopore immersed in electrolyte solution is shown in Fig. 1. Due to their tunable geometry, resistance to harsh solvent conditions, and selectivity for surface modifications, solid-state nanopores have become promising tools in nanotechnology. To date, they have been used to unravel the fined-tuned interactions between DNA and proteins9,10, to build proteinaceous channels that mimic trafficking in cells11 and to evaluate the feasibility of a fast-sequencing DNA technique1214. Besides this wide range of applications, solid-state nanopores have also been used as models to test fundamental theories about ion dynamics in nanoscale confinements15,16. Studying ionic transport on the nanoscale is crucial not only for controlling the flow of ions across the nanopore1520, but also for understanding how biological channels function. In this regard, theoretical and computational modeling have been of critical importance, since experimental techniques are not yet able to resolve the dynamics of the few ions contained in the nanopore volume.

Fig. 1
Atomic model of a polymer nanopore. The figure shows a polyethylene terephthalate (PET) nanopore immersed in electrolyte solution. The nanopore bulk material is shown in gray color. PET residues with deprotonated carboxyl groups are colored in purple. ...

Modeling based on the Poisson-Nernst-Planck (PNP) equations, Monte Carlo (MC) as well as molecular dynamics (MD) simulations have accompanied experimental efforts. The continuum approach of PNP successfully described the effect of ion current rectification observed in assymetric pores2123 as well as ionic selectivity, i.e., a pore preference for transporting one type of ion, cations or anions24,25. However, continuum modeling cannot provide a molecular understanding of how ions compete for space and surface charges in the restricted volume of a nanopore. Monte Carlo simulations, which take into account the finite size of ions, offer more detail on the behavior of ions in nanopores. Applied to calcium-selective channels these simulations successfully described channel selectivity properties26. MD approaches reveal most details about ionic trajectories and their interactions with other ions, water and nanopore surface in equilibrium as well as non-equilibrium conditions. Increases in computational power have extended the duration of MD simulations, approaching today a time scale that is accessible in ion current measurements. Indeed, MD simulations successfully described transport properties of proteinaceous nanopores such as α-hemolysin2729 and solid-state silicon-based nanopores3033. The MD methodology had not been extended initially to polymer nanopores due to the lack of an accurate model for the chemical features of the pore surface, a key factor that determines the transport properties.

Here, we report all-atom simulations of ion transport through a polyethylene terephthalate (PET) nanopore. PET nanopores have emerged as a nano-scale model system of interest, because of the experimental observation of multiple effects allowing for the control of ion transport properties, i.e., ion current rectification8,34, selectivity22,24,25,35, as well as the realization of ionic diodes19,20 and ionic transistors36.

Recently, it has been shown that local charge inversion can be a powerful tool to control ionic transport properties in PET nanopores37. Charge inversion means that the total positive charge brought by the cations close to the surface becomes larger than the total negative surface charge on the walls3841. Charge inversion can be observed through reversal of the sign of the electrophoretic mobility of charged colloids42,43. In a nanopore with negative walls, charge inversion results in switching the selectivity properties of the nanopore. While a negatively charged nanopore in contact with e.g., KCl, would be cation selective, the same nanopore in contact with multivalent cations, such as Ca2+, will become selective for anions. Such a selectivity switch induced through charge inversion has also been observed in protein pore OmpF44 as well as in PET nanopores37. Therefore, studying the interactions of multivalent and monovalent ions with the pore structure is crucial for gaining a better understanding of induced charge inversion.

In this article, we present a procedure to build all-atom models of PET nanopores and use it to accurately replicates experimental devices. After validation of the model, we study the dynamics of K+, Cl and Ca2+ ions inside the nanopore using the MD approach. The following questions are discussed: How does the PET nanopore surface affect the ratio of the current carried by K+, Cl and Ca2+ ions? How does Ca2+ interact with the charged groups of the pore walls?.

2 Computational methods

In this section, the protocol to build PET nanopores is presented, along with the needed procedures for solvation/ionization, and setup of ionic conduction simulations.

The MD simulations were performed using the program NAMD45, the post-processing analysis of the MD trajectories was performed with custom VMD46 and MatLab47 scripts. The protocols for the MD simulations have been described in detail before33,45 and are briefly stated here. For all simulations, an integration time step of 1 fs was used. For periodic boundary conditions, full electrostatic forces were computed using the particle-mesh-Ewald method with a grid density of 1/Å3 or finer. For van der Waals and electrostatic interactions in non-periodic boundary conditions, a cut-off of 12 Å with a switching function starting at 10 Å was assumed. Langevin dynamics was utilized to control system temperature. For MD simulations in the NpT ensemble, the pressure (p) was kept at 1 atm using a hybrid Nosé-Hoover Langevin piston. Force field parameters for K+, Cl, and Ca+2 ions were taken from the CHARMM force field48; the TIP3P49 model was used for water.

2.1 PET model

A PET monomer can exist as a trans or a cis isomer, depicted in Fig. 2a and 2b, respectively. Monomers connect to form polymers and terminate in carboxyl groups that are protonated at pH 4 and deprotonated at pH 7 (cf. Fig. 2c and 2d). Several MD models of PET polymers had already been developed5052 and successfully used to study the atomic origin of PET mechanical properties5355. However, these models had been designed to describe properties of the bulk material and used force fields incompatible with explicit water models. To develop a PET model suitable for simulation of ion currents in a nanopore system, we built the PET moiety with the widely used CHARMM force field48. This force field was emplyed extensively to study a variety of biomolecular channels29,5658, the limitations of the model in representing interactions with water and ions being well-known. The force field parameters for PET were obtained in analogy to previously parameterized model compounds. Topology and parameter files as well as atomic models of the cis and trans isomers are provided in Supplementary material.

Fig. 2
PET chemical structure. Shown are the chemical structure of a PET monomer in (a) trans and (b) cis isomeric forms. Also shown are the terminal PET monomers with (c) protonated and (d) deprotonated carboxyl groups, respectively.

In previous studies55,59, PET bulk models were built using a stepwise procedure60, where monomers are added to a growing chain confined within a periodic cell. A disadvantage of such a method is that the accessible space is reduced as the polymer grows, which can lead to residues crossing through the center of aromatic rings. To avoid this problem, we built a periodic PET bulk using a collapsing/annealing MD procedure, presented schematically in Fig. 3.

Fig. 3
Building bulk PET. The figure illustrates the collapsing/annealing procedure used to build amorphous PET. The first stage of the procedure involved collapsing a grid of PET 9-mers. (a) Model of a PET 9-mer. (b) Initial conformation of the PET grid. Ten ...

First, ten PET polymers were generated. The number of subunits for each polymer was set to 9, corresponding to the most abundant polymer length observed in amorphous PET materials55. The initial conformation for each 9-mer was linear, the isomeric form (cis or trans) for each residue was assigned randomly, and carboxyl groups at both ends were protonated. To randomize the structure, each 9-mer was minimized for 200 steps and equilibrated for 1 ns at 1000 K. A single PET 9-mer obtained by this method is shown in Fig. 3a. Subsequently, one of the ten PET 9-mers was randomly selected, reoriented in space, and assigned to a point on a 7×7×5 grid with 75 Å spacing between grid points. After that, the grid was collapsed into a 90 Å ×90 Å ×60 Å rectangular box by applying forces of 5 pN towards the center of the box to each aromatic carbon atom located outside of the box. This simulation was performed at 1000 K to avoid folding of the 9-mers when they were moving toward the rectangular region. The collapsing grid at 0, 0.4, 0.7 and 2 ns is shown in Fig. 3b, 3c, 3d and 3e, respectively. The collapsed structure was placed within a periodic cell of dimension 120 Å ×120 Å ×90 Å, i.e., sufficiently large to avoid contacts among periodic images, as shown in Fig. 3g. To obtain homogeneous density, the PET structure was annealed in the NpT ensemble. First, the PET cell was simulated for 5 ns at 1000 K. At this temperature, the PET structure is a fluid and fills the entire cell volume. Snapshots of that simulation after 1 and 5 ns are shown in Fig. 3h and 3i, respectively. After that, the temperature was gradually reduced to 400 K, with a cooling step of 0.2 ns per −100 K. The annealing cycle finished with 5 ns of equilibration at 300 K. The resulting amorphous PET is shown in Fig. 3j, the periodic cell having a dimension of 107 Å ×84 Å ×64 Å.

A nanopore was molded into the PET material by defining a conical pore geometry8, according to

x2+y2=[(z+h2)tanθ+Rmin]2,h2zh2
(1)

Here, θ is the angle between nanopore wall and pore axis, h is the height of the pore and Rmin is the minimum radius of the cone. Through deletion of residues, a pore was created with Rmin=16 Å, h=105 Å, θ=3°. The height of the pore was achieved through replication of the bulk structure unit cell, followed by deletion of residues above z = h/2 and below z = −h/2. Protonated carboxyl groups were added to all polymer terminals.

Deletion of residues cut some of the 9-mers into smaller segments. To prevent that small segments would be released from the nanopore surface, polymers with a length of three subunits or less were removed. To relax the surface while preserving the optimal geometry, a force of 1 pN guided non-terminal residues entering the pore volume back into the PET bulk during 2000 steps of minimization and 1 ns of simulation. After that, the system was equilibrated in a constant volume simulation without guiding forces for 4 ns. Figure 4b shows the radius of the pore along the Z-axis after the initial cut (circles) and after surface rearrangement (squares). As it can be observed, the conical geometry is preserved during the entire procedure. The narrowest pore region has a radius of 1 nm, resembling the smallest nanopores that can be manufactured by track-etching62.

Fig. 4
Nanopore radius. Panel (a) shows a snapshot of a solvated PET nanopore. The PET structure is pictured as gray beads and water as a white surface. Panel (b) shows the radius of the PET nanopore along the Z-axis at two different stages of creation and after ...

2.2 Ionic current simulations

Simulations of ionic currents closely follow protocols used in earlier MD studies 29,33,45,5658. The PET nanopore was solvated in a rectangular water box using the Solvate plugin of VMD46. To allow the solvent molecules fill all small voids at the PET surface, the system was simulated for 5 ns at 300 K in the NpT ensemble, with the PET structure restrained using harmonic forces with a spring constant of 5 kcal/mol/Å2. Subsequently, the restraints were removed, and the system simulated further for 7 ns at 300 K in the NpT ensemble. The radius of the pore along the Z-axis after solvation is shown in Fig. 4b (solid line).

The solvated nanopore was used to study ionic conduction at pH 4 and pH 7. To mimic pH 7 conditions, terminal PET residues in contact with water were deprotonated, resulting in a surface charge of −1.03 e nm−2, which agrees well with experimental values63. The system was neutralized by adding counter ions, either K+ or Ca2+, using the program cionize64 that located the ions according to the electrostatic potential of the PET structure. To mimic pH 4 conditions, all carboxyl groups remained protonated. Subsequently, K+ and Cl ions were randomly placed in solvent compartments, corresponding to 1 M KCl concentration. After adding ions, the systems were equilibrated at 300 K for 1 ns in the NpT ensemble, followed by 1 ns equilibration in the NV T ensemble. The final systems have periodic lattice vectors of 107 Å ×84 Å×265 Å.

To simulate ionic currents, a uniform electrostatic field EZ was applied to all atoms of the system along the Z-axis. The voltage difference [var phi] across the simulated cell is

φ=LZ·EZ,
(2)

where Lz is the dimension of the system in the Z-direction. Table 1 presents a summary of the ion conduction simulations performed. The voltage signs are defined relative to the lower solvent compartment which is in contact with the narrowest pore opening. A +1 V bias induces cations to move from the lower to the upper solvent compartment. To prevent the PET membranes from drifting, the benzene carbons located in a belt region of 40 Å height and 30 Å away from the pore axis were restrained by harmonic forces with a spring constant value of 1 kcal/mol/Å2.

Table 1
Ionic current simulations of PET nanopores. The table shows the conditions for the various simulations. The number of K+ ions was 994 for all simulations. The number of Ca2+ ions in simulations 5–8 is sufficient to saturate half of the carboxyl ...

3 Results and Discussion

In this section we present and discuss the results of our simulations. First, we compare our PET bulk model with previously reported PET simulations and the expected values from experimental PET materials. Then, we report the ionic current for different pH conditions and electrolyte solution and relate the observed ionic currents to the individual ionic behavior and surface properties of the pore.

3.1 PET structure

Experimentally, single nanopores in PET are prepared by the track-etching technique described previously7. Briefly, PET films are irradiated with single swift heavy ions, which were accelerated to a total kinetic energy of 2.2 GeV (UNILAC, GSI Darmstadt, Germany). This irradiation process causes the formation of a single damage track through the films. It is expected that the impact of the heavy ion modifies the biaxially-oriented PET into an isotropic amorphous structure65. Conically shaped nanopores are obtained by asymmetric etching of the irradiated foils in 9 M NaOH while the other side of the membrane is in contact with an acidic stopping medium6. The heavy ion irradiation and chemical etching cause formation of 1 carboxyl group per nm2 63.

An atomic-model of a PET nanopore should be built from amorphous bulk material, rather than from crystalline material. The procedure for building an amorphous bulk cube has been described in Computational Methods (see Fig. 3). The final step of this procedure involves an annealing cycle, shown in Fig. 5a, with 5 ns equilibration at 1000 K, cooling in 100 K steps with 0.2 ns equilibration per step, and 5 ns equilibration at room temperature (300 K).

Fig. 5
Analysis of the annealing cycle for periodic bulk PET. Panel (a) shows the variation of temperature during the annealing cycle and the respective variation in density (b), characteristic ratio (c) and percentage of glycol dihedral in trans conformation ...

To compare our bulk material with other PET models and experimental PET material, we analyzed three measures commonly used to characterize PET bulk properties, namely, volume density, the so-called characteristic ratio for a 9-mer, C9, and the glycol dihedral angle between two joint PET residues.

The density of our PET bulk model during the annealing cycle is shown in Fig. 5b. At the initiation of the annealing cycle, densities in the center (x marks) and periphery (solid line) of the bulk cube differed significantly, but both converged to a constant value, demonstrating homogeneity of the bulk material. The equilibrium density of 1.20 g cm−3 is in good agreement with previously reported PET models51,52,55 with densities of 1.28 g cm−3 and 1.13 g cm−3, but still slightly lower than the experimental value53 of 1.34 g cm−3.

The characteristic ratio is defined as the ratio between the mean-square end-to-end length R of a polymer and the length li of each rigid bond within the polymer51. We calculate the characteristic ratio C9 according to

C9=R2ili2.
(3)

At 1000 K, the PET bulk behaves as a fluid and C9 fluctuates, reaching values up to 5. As the system becomes more dense, C9 decreases, obtaining a final value of 3.96, which is close to the experimental value of 4 and in good agreement with previously reported simulations51.

Panel 5d shows the percentage of glycol dihedral in the trans conformation, the trans angle being defined as 180° ± 20°. At the end of the annealing cycle, we obtained 39.4% of trans content. Here we found our largest discrepancy with previous atomistic models and experimental values that report 25% trans55. The percentage of trans content is strongly influenced by the values of the O-C-C-O dihedral parameters. Cail et al.54 have reported that for dihedrals with barrier heights of 1.40, 1.00 and 0.24 kcal mol−1, the trans content at 541 K is 19%, 24% and 39%, respectively. The value we used was taken from the propane-1,2-diyl diacetate compound model present in CHARMM, which has a barrier height of 0.2 kcal mol−1 (dihedral parameters available in Supplementary material). In order to obtain 25% trans content, a further refinement of this dihedral parameter would be needed.

Finally, the accuracy of the bulk model was evaluated using the pair distribution function g(r). Figure 6 shows g(r) for all carbon-carbon pairs. The overall shape of the curve compares well with previous studies51. The peaks observed between 1 and 4 Å describe the carbon- carbon distances defined by the connectivity within each polymer60. The inset in Fig. 6 shows the g(r) computed for carbon-carbon pairs, where each carbon belongs to a different polymer chain. The shoulder around 5 Å is related to the van der Waals contact between two carbons that are not covalently bonded.

Fig. 6
Pair distribution function, g(r), of carbon-carbon pairs separated by a distance r. g(r) was computed over the last 2 ns of a 5 ns MD equilibration at 300 K. The inset shows the g(r) for carbon-carbon pairs, each carbon belonging to different PET 9-mers. ...

From the PET bulk structure we built a conical PET nanopore with a minimum radius of 10 Å and 100 Å height and with a surface charge density as observed in experiments 63. The detailed procedure is described in Computational Methods.

3.2 Ionic currents

Previous sections presented the results from building an amorphous form of the PET polymer material and from modeling a nanopore with conical shape and surface charge density closely matching the experimentally studied systems. Following the validation of material characteristics, we used the nanopore system to study its transport properties. Table 1 lists the ionic conduction simulations that were performed under various electrolyte conditions.

As mentioned in the introduction, we employed MD simulations to obtain information about how the surface charge of the pore walls affects the ion current passing through the pore, and about how the presence of divalent cations such as Ca2+ changes the ionic selectivity of the pore.

We performed ten simulations, varying pH values between pH 4, where the carboxyl groups are protonated, and pH 7, where the carboxyl groups are deprotonated, also varying voltage biases between +1 V and −1 V, and choosing different Ca2+ concentrations. For each simulation, the instantaneous ionic currents were computed using

I(t)=1ΔtLzi=1Nqi[zi(t+Δt)zi(t)],
(4)

where qi is the charge of ion i, zi its Z-coordinate as a function of time t, Lz is the length of the system along the Z-axis, and Δt is 1 ps. The averaged ionic currents over the last 5 ns of each MD simulation are presented in Table 1. Individual ionic currents for each MD simulation are presented in Fig. 7.

Fig. 7
Instantaneous ionic currents through PET nanopores. The figure shows the ionic currents computed using Equation 4, each point representing an average over 0.5 ns. The ionic currents were calculated for each ionic species, K+ (blue), Cl (red) ...

Our simulations suggest a strong dependence of ion selectivity properties of PET nanopores on the nanopore surface charge. At pH 4, when the nanopore surface is not charged, the anionic and cationic currents are of about the same magnitude (simulations 1 and 2), as shown in Fig. 7a and Table 1. As expected, there is no ion current selectivity in a neutral nanopore. At pH 7 however, where the carboxyl groups are negatively charged, the ion current carried by K+ is about three times larger than the current carried by Cl (simulations 3 and 4), as shown in Fig. 7b and Table 1. Ion selectivity is typically expressed as the ratio of the current carried by cations and the total current carried by cations and anions. The ion selectivity value of 0.79, as predicted by our simulations, agrees well with the values measured experimentally at elevated KCl concentrations37,66. Due to the comparatively small volume accessible in all atom simulations, the ion concentration of 1 M KCl is chosen larger than the concentrations typically used in experimental measurements. However, experimental values of ion currents and ion selectivity have been determined at a very wide range of KCl concentrations including 1 M KCl 37,66. One should note that due to the required electroneutrality of the system, the number of Cl ions changes from 995 in simulations 1 and 2 (pH 4) to 703 in simulations 3 and 4 (pH 7). This change can account for some of the reduction in anionic current. However, the reduction in anionic current by a factor of two is larger than would be expected solely due to the changed concentration.

The pH dependence of the ion selectivity properties of the PET pore can also be seen by analyzing the concentration profiles of different ion species as a function of distance from the center of the nanopore (Fig. 8a and 8b). At pH 4, the concentrations of K+ and Cl ions are practically identical, while at pH 7 there is a huge increase of the K+ concentration, which agrees well with the observed increase of the cationic current (Table 1).

Fig. 8
Ionic concentrations as a function of radial distance from the nanopore center. Concentrations are averaged over an 8 nm [−40; 40 Å] compartment along the Z-axis. The figures show the ion concentration inside the pores for each MD simulation ...

We also examined ion currents in the presence of both KCl and CaCl2. In a control condition at pH 4 (simulation 5 and 6), with an uncharged surface, addition of CaCl2 had little effect on the total current and the ion selectivity when compared with a pure KCl solution (simulation 1 and 2). The cationic current was simply split between K+ and Ca2+ ions with a ratio roughly corresponding to their relative concentrations. The increase of the anionic current compared to the results of simulations 1 and 2 is due to the increase of the Cl concentration (Fig. 8c).

At pH 7, with negatively charged surfaces due to the deprotonated carboxyl groups in the system, the addition of CaCl2 resulted in a more complex and interesting behavior. In simulations 7 and 8, 145 Ca2+ ions were added to neutralize the surface charge. In simulations 9 and 10, 291 Ca2+ were added, which could, in principle, saturate every carboxyl group, thereby possibly resulting in overcharging of the system walls. In all of the simulations 7–10, Ca2+ competed with K+ ions for locations close to the negative carboxyl groups (see below)67. During the simulations, most Ca2+ ions remained practically immobile, in contact with carboxyl groups. For the smaller number of 145 Ca2+ ions, indeed only a very small fraction of Ca2+ ions was mobile and their contribution to the total ion current was negligible, as shown in Fig. 7c and Table 1. With Ca2+ ions linked to carboxyl groups, the system is effectively neutralized (simulations 7 and 8). Differences between cationic and anionic currents are within the standard error of the simulations, indicating that there is no selectivity of ion current in this case.

Addition of 291 Ca2+ ions (simulations 9 and 10) resulted in a further change of the ionic conduction behavior. The currents carried by Cl are now larger than the combined K+ and Ca2+ currents (Table 1 and Fig. 7d). This suggests that the nanopore in the latter system has changed selectivity, possibly due to charge inversion in the nanopores. To analyze the reason underlying the changed ionic currents, we determined the ionic concentration profile in the pore as well as the electrostatic potential along the Z axis of the nanopore.

An analysis of the ion densities in the nanopore shows that Ca2+ outcompeted K+ in the region close to the pore walls, the K+ concentration being significantly less than in the situation without Ca2+ (compare Fig. 8b with 8d and 8e). However, the Ca2+ concentration at the pore walls increased only slightly when the overall concentration was changed from 145 to 291 ions. This indicates that the fraction of the carboxyl groups with an adsorbed Ca2+ did not change much for the higher Ca2+ concentrations, i.e., the pore walls remained roughly electroneutral. But there arose an increase of the total Ca2+ concentration in the pore, accompanied by a corresponding increase in the amount of chloride ions. The distribution of the Cl ions, however, does not show any maximum in the vicinity of the position with the increased Ca2+ concentration. Such maximum of Cl concentration near the pore wall is typically indicative of charge inversion68. Since one does not see such a maximum here, one can rule out overcharging of the nanopore walls as the cause for increased anionic currents (see Table 1). We suggest then that the slight increase of the adsorbed Ca2+ causes changes in the electrostatic potentials. These changes, together with the increased concentration of Cl ions, can explain the increased anionic current. Our results are in agreement with experimental data40,67 and previous modeling69, which showed that the charge inversion and ion selectivity switch do not occur if a system contains a high concentration of monovalent cations. This effect, called salting-out, is explained by the competition of K+ and Ca2+ ions for the negative charges at the walls. The K+ ions concentration at the wall remains high enough so that surface overcharging does not occur.

To study how the adsorbed ions affect the electrostatic properties of the nanopore channel, we computed the electrostatic potential along the pore Z-axis, Vsystem(Z). The procedure used has been described in detail before33. Briefly, the averaged three-dimensional electrostatic potential was computed for the entire system using the PMEpot plugin29 available in VMD46 for the last 5 ns of each trajectory. From the three-dimensional electrostatic potentials, we extracted the electrostatic potential across the pore Z-axis at +1 V and −1 V biases, denoted as V+1V (Z) and V−1V (Z), respectively. Then, Vsystem(Z) was calculated using:

Vsystem(Z)=V+1V(Z)+V1V(Z)2.
(5)

Adding V+1V (Z) and V−1V (Z) cancels out the influence of the external applied voltage biases +1 V and −1 V, leaving only the electrostatic contributions from all charged species in the system. Therefore, Vsystem describes the electrostatic potential of the system itself acting on a line coincident with the pore Z-axis.

Figure 9 shows Vsystem for all systems. There is a significant change in Vsystem at different pH conditions. At pH 4, the PET surface is not charged, there is no ion adsorption, and the resulting Vsystem is small and positive (Fig. 9a and 9c). However, at pH 7, the PET surface becomes negatively charged and in the presence of 1 M KCl, Vsystem becomes negative (Fig. 9b), which implies cation selectivity of the pore. After adding Ca2+, the negatively charged surface binds the Ca2+ ions, Vsystem becomes positive again as shown in Fig. 9d and 9e. The reversal of the sign of Vsystem after Ca2+ binding suggests that charge inversion has occurred; however, we do not see selectivity for Cl ions, a characteristic feature of charge inversion. The absence of Cl selectivity can not be observed in our simulations due to the high ionic concentration assumed; only at a sufficiently low ionic concentration would the effect of charge inversion and Cl selectivity arise.

Fig. 9
Electrostatic potential of the system Vsystem(Z). Vsystem(Z) is defined in Equation 5. Figures show Vsystem(Z) for the systems at: (a) pH 4 and 1 M KCl, (b) pH 7 and 1 M KCl, (c) pH 4 and 1 M KCl and 145 Ca2+ ions, (d) pH 7 and 1 M KCl and 145 Ca2+ ions, ...

4 Conclusions

In this article, we have reported the first all-atom MD simulations of PET polymer nanopores. The study presents several methodological and scientific accomplishments and opens an avenue for future research. First, we have presented a collapsing-annealing procedure to produce polymeric PET bulk structures. This procedure overcomes the difficulties identified for the growing-chain methodology60 and can be easily used to build other polymeric bulk structures. Secondly, we have built a PET nanopore model that accurately reproduces geometry and surface charge density of the experimental systems. Third, we have built an MD model based on the CHARMM force field, which permits extension for future simulations, for instance, to study functionalized PET nanopores20,70. Our all-atom MD simulations give insight into transport properties of these pores. We have described the ion dynamics inside PET nanopores immersed in KCl solution at different pH conditions, which correspond to different surface charge properties of the pore walls. Our results reproduced the cation selectivity observed at neutral pH conditions, at which the PET surface is negatively charged. Furthermore, we have calculated the ionic current carried by each ion species separately as well as the ionic distribution inside the pore. We have shown a very high affinity of Ca2+ ions to the surface deprotonated carboxyl groups, which is important for understanding the effect of charge inversion. Our simulations showed that mobility of Ca2+ ions in the pore is very low and that Ca2+ ions close to the carboxyl groups are practically immobile. Since our simulations were performed in the presence of a high concentration of KCl, charge inversion could not be observed.

Supplementary Material

ESI

Acknowledgments

We thank the members of the Theoretical and Computational Biophysics Group for helpful discussions. Z. S. is an Alfred P. Sloan Fellow. This work was supported by grants from NIH (P41-RR05969) and NSF (CCR-02-10843, CHE-0747237). We acknowledge supercomputer time provided at the National Center for Supercomputing Applications through Large Resources Allocation Committe grant MCA93S028 and provided on the Turing Xserve Cluster at the University of Illinois at Urbana-Champaign.

References

1. Li J, Stein D, McMullan C, Branton D, Aziz MJ, Golovchenko JA. Nature. 2001;412:166–169. [PubMed]
2. Heng JB, Ho C, Kim T, Timp R, Aksimentiev A, Grinkova YV, Sligar S, Schulten K, Timp G. Biophys J. 2004;87:2905–2911. [PubMed]
3. Storm AJ, Chen JH, Ling XS, Zandergen HW, Dekker C. Nat Mater. 2003;2:537–540. [PubMed]
4. Park SR, Peng HB, Ling XSS. Small. 2007;3:116–119. [PubMed]
5. Fleischer RL, Price PB, Walker RM. Nuclear tracks in solids. Principles and applications. University of California Press; Berkeley: 1975.
6. Apel PY, Korchev YE, Siwy Z, Spohr R, Yoshida M. Nucl Instr Meth Phys Res B. 2001;184:337–346.
7. Spohr R. German Patent DE. U S Patent. 2951376 C2. 4369370. 1983.
8. Siwy Z. Adv Funct Mater. 2006;16:735–746.
9. Zhao Q, Sigalov G, Dimitrov V, Dorvel B, Mirsaidov U, Sligar S, Aksimentiev A, Timp G. Nano Lett. 2007;7:1680–1685. [PMC free article] [PubMed]
10. Smeets RMM, Kowalczyk SW, Hall AR, Dekker NH, Dekker C. Nano Lett. 2008 doi: 10.1021/nl803189k. [PubMed] [Cross Ref]
11. Jovanovic-Talisman T, Tetenbaum-Novatt J, McKenney AS, Zilman A, Peters R, Rout MP, Chait BT. Nature. 2009;457:1023–1027. [PMC free article] [PubMed]
12. Healy K. Nanomedicine. 2007;2:459–4811. [PubMed]
13. Branton D, Deamer DW, Marziali A, Bayley H, Benner SA, Butler T, Ventra MD, Garaj S, Hibbs A, Huang X, Jovanovich SB, Krstic PS, Lindsay S, Ling XS, Mastrangelo CH, Meller A, Oliver JS, Pershin YV, Ramsey JM, Riehn R, Soni GV, Tabard-Cossa V, Wanunu M, Wiggin M, Schloss JA. Nature Nanotech. 2008;26:1146–1153. [PMC free article] [PubMed]
14. Sigalov G, Comer J, Timp G, Aksimentiev A. Nano Lett. 2008;8:56–63. [PMC free article] [PubMed]
15. Lev A, Korchev Y, Rostovtseva T, Bashford C, Edmonds D, Pasternak C. Proc R Soc Lond B. 1993;252:187–192. [PubMed]
16. Gillespie D, Boda D, He Y, Apel P, Siwy ZS. Biophys J. 2008;95:609–619. [PubMed]
17. Dekker C. Nature Nanotech. 2007;2:209–215. [PubMed]
18. Daiguji H, Oka Y, Shirono K. Nano Lett. 2005;5:2274–2280. [PubMed]
19. Karnik R, Duan C, Castelino K, Daiguji H, Majumdar A. Nano Lett. 2007;7:547–551. [PubMed]
20. Vlassiouk I, Siwy Z. Nano Lett. 2007;7:552–556. [PubMed]
21. Cervera J, Schiedt B, Ramírez P. Europhys Lett. 2005;71:35–41.
22. Cervera J, Schiedt B, Neumann R, Mafé S, Ramírez P. J Chem Phys. 2006;124:104706-1–104706-9. [PubMed]
23. White SH, Bund A. Langmuir. 2008;24:2212–2218. [PubMed]
24. Plecis A, Schoch RB, Renaud P. Nano Lett. 2005;5:1147–1155. [PubMed]
25. Vlassiouk I, Smirnov S, Siwy Z. Nano Lett. 2008;8:1978–1985. [PubMed]
26. Boda MVD, Eisenberg B, Nonner W, Henderson DJ, Gillespie D. Phys Rev Lett. 2007;98:168102-1–168102-4. [PubMed]
27. Wells D, Abramkina V, Aksimentiev A. J Chem Phys. 2007;127:125101–125110. [PMC free article] [PubMed]
28. Mathé J, Aksimentiev A, Nelson DR, Schulten K, Meller A. Proc Natl Acad Sci USA. 2005;102:12377–12382. [PubMed]
29. Aksimentiev A, Schulten K. Biophys J. 2005;88:3745–3761. [PubMed]
30. Heng JB, Aksimentiev A, Ho C, Marks P, Grinkova YV, Sligar S, Schulten K, Timp G. Nano Lett. 2005;5:1883–1888. [PMC free article] [PubMed]
31. Heng JB, Aksimentiev A, Ho C, Marks P, Grinkova YV, Sligar S, Schulten K, Timp G. Biophys J. 2006;90:1098–1106. [PubMed]
32. Zhao Q, Comer J, Yemenicioglu S, Aksimentiev A, Timp G. Nucl Acids Res. 2008;36:1532–1541. [PMC free article] [PubMed]
33. Cruz-Chu ER, Aksimentiev A, Schulten K. J Phys Chem C. 2009;113:1850–1862. [PMC free article] [PubMed]
34. Siwy Z, Fuliński A. Phys Rev Lett. 2002;89:198103-1–198103-4. [PubMed]
35. Siwy Z, Gu Y, Spohr HA, Baur D, Wolf-Reber A, Spohr R, Apel P, Korchev YE. Europhys Lett. 2002;60:349–355.
36. Fan R, Yue M, Karnik R, Majumdar A, Yang P. Phys Rev Lett. 2005;95:086607-1–086607-4. [PubMed]
37. He Y, Gillespie D, Boda D, Vlassiouk I, Eisenberg RS, Siwy ZS. J Am Chem Soc. 2009;131:5194–5202. [PMC free article] [PubMed]
38. Besteman K, Zevenberger MAG, Heering HA, Lemay SG. Phys Rev Lett. 2004;93:170802-1–170802-4. [PubMed]
39. Besteman K, Zevenberger MAG, Lemay SG. Phys Rev E. 2005;72:061501-1–061501-5. [PubMed]
40. van der Heyden FHJ, Stein D, Besteman K, Lemay SG, Dekker C. Phys Rev Lett. 2006;96:224502-1–224502-4. [PubMed]
41. van der Heyden FHJ, Stein D, Besteman K, Lemay SG, Dekker C. J Phys Chem C. 2007;111:15575–15585.
42. James R, Healy T. J Coll Interf Sci. 1972;40:53–64.
43. Martín-Molina A, Quesada-Pérez M, Galisteo-González F, Hidalgo-Álvarez R. J Chem Phys. 2003;118:1–7.
44. Alcaraz A, Nestorovich EM, López ML, García-Giménez E, Bezrukov SM, Aguilella VM. Biophys J. 2009;96:56–66. [PubMed]
45. Phillips JC, Braun R, Wang W, Gumbart J, Tajkhorshid E, Villa E, Chipot C, Skeel RD, Kale L, Schulten K. J Comp Chem. 2005;26:1781–1802. [PMC free article] [PubMed]
46. Humphrey W, Dalke A, Schulten K. J Mol Graphics. 1996;14:33–38. [PubMed]
47. Inc. The MathWorks. Matlab v.6. 2002.
48. MacKerell A, Jr, Bashford D, Bellott M, Dunbrack RL, Jr, Evanseck J, Field MJ, Fischer S, Gao J, Guo H, Ha S, Joseph D, Kuchnir L, Kuczera K, Lau FTK, Mattos C, Michnick S, Ngo T, Nguyen DT, Prodhom B, Reiher IWE, Roux B, Schlenkrich M, Smith J, Stote R, Straub J, Watanabe M, Wiorkiewicz-Kuczera J, Yin D, Karplus M. J Phys Chem B. 1998;102:3586–3616. [PubMed]
49. Jorgensen WL, Chandrasekhar J, Madura JD, Impey RW, Klein ML. J Chem Phys. 1983;79:926–935.
50. Sun H. J Comp Chem. 1994;15:752–768.
51. Hedenqvist MS, Bharadwaj R, Boyd RH. Macromolecules. 1998;31:1556–1564.
52. Boyd SU, Boyd RH. Macromolecules. 1998;31:1556–1564.
53. Zhou J, Nicholson TM, Davies GR, Ward IM. Comput Theor Polym Sci. 2000;10:43–51.
54. Cail JI, Stepto RFT. Polymer. 2003;44:6077–6087.
55. Roberge M, Prud’homme RE, Brisson J. Polymer. 2004;45:1401–1411.
56. Khalili-Araghi F, Tajkhorshid E, Schulten K. Biophys J. 2006;91:L72–L74. [PubMed]
57. Sotomayor M, Vasquez V, Perozo E, Schulten K. Biophys J. 2007;92:886–902. [PubMed]
58. Gumbart J, Schulten K. J Gen Physiol. 2008;132:709–719. [PMC free article] [PubMed]
59. Gestoso P, Brisson J. Polymer. 2003;44:7765–7776.
60. Theodorou DN, Suter UW. Macromolecules. 1985;18:1467–1478.
61. Smart OS, Neduvelil JG, Wang X, Wallace BA, Sansom MSP. J Mol Graphics. 1996;14:354–360. [PubMed]
62. Siwy Z, Apel P, Baur D, Dobrev DD, Korchev YE, Neumann R, Spohr R, Trautmann C, Voss K. Surf Sci. 2003;532–535:1061–1066.
63. Wolf-Reber A. Ph D Thesis. University of Frankfurt; Frankfurt, Germany: 2002.
64. Stone JE, Phillips JC, Freddolino PL, Hardy DJ, Trabuco LG, Schulten K. J Comp Chem. 2007;28:2618–2640. [PubMed]
65. Singh V, Singh T, Chandra A, Bandyopadhyay SK, Sen P, Witte K, Scherer UW, Srivastava A. Nucl Instr Meth Phys Res B. 2006;244:243–247.
66. Ramírez P, Gómez V, Cervera J. J Chem Phys. 2007;126:194703-1–194703-9. [PubMed]
67. van der Heyden FHJ, Stein D, Dekker C. Phys Rev Lett. 2005;95:116101-1–116101-4. [PubMed]
68. Torrie GM, Valleau JP. J Phys Chem. 1982;86:3251–3257.
69. Chen Y, Ni Z, Wang G, Xu D, Li D. Nano Lett. 2008;8:42–48. [PubMed]
70. Xia F, Guo W, Mao Y, Hou X, Xue J, Xia H, Wang L, Song Y, Ji H, Ouyang Q, Wang Y, Jiang L. J Am Chem Soc. 2008;130:8345–8350. [PubMed]