The hopping rates defined by the CTRW model,
The accelerated appearance of applications and theoretical studies of carbon nanotubes (CNTs)1 in the literature over the last two decades reveals their broad attraction in fields as diverse as nanoelectronics, bionanotechnology and solar energy conversion nanodevices2-5. In particular, they provide the basis for simplified models of nanofluidic machines and membrane channels such as aquaporins6-9 and have been employed in a range of simulation studies at the nanoscale to investigate the transport of water undergoing single-file or highly confined diffusion and permeation in a single CNT or in arrays of CNTs10-22.
The simulation method of choice in the study of water transport through nanostructures is molecular dynamics (MD) (see, for example, [23-25]) which provides real-time trajectories of all of the atoms within the system under investigation and, through the statistical framework of linear response theory, permits a direct link with macroscopic observables. One of the first studies of water confined within CNTs employing MD, supplemented with an analysis using a continuous-time random walk model (CTRWM), was reported by Berezhkovskii and Hummer12,13 and demonstrated that water fills a (6,6) CNT spontaneously in a single-file fashion. In later work, MD simulations focused on pressure-driven osmotic water flow in single CNTs and arrays of CNTs14-16 and in model pores containing charge distributions which are known to be ubiquitous for biological channels17,18. Water channels made of double-walled carbon nanotubes embedded in a lipid membrane, where the inner tube is modeled as a (6,6) CNT with a charge distribution emulating the NPA motif of aquaporins19, have also been investigated, and anomalous observations relating to the dynamics, configurational states and apparent pore-length dependence of the properties of water confined within CNTs have received attention20-22.
In [26,27], molecular simulations were reported for unidirectional single-file diffusion of water within (5,5) and (6,6) single-walled carbon nanotubes (SWCNT) and normal diffusion in (8,8) and (11,11) SWCNT’s across a lipid membrane. These studies employed both equilibrium MD simulation techniques and an adaptation of MD to investigate conditions involving the influence of both static and alternating electric fields. While the fields employed in this work were relatively weak vis-à-vis fields generated by charge distributions, it was observed that the permeation of water for the (5,5) SWCNT embedded in a lipid membrane is enhanced by static electric fields due to a decrease in the fluctuations of the number of water molecules inside the SWCNT. In the larger single-file (6,6) CNT pore, dipole alignment with the imposed external field is more easily achieved due to a greater degree of rotational freedom of the water molecules and this was observed to result in a significant reduction in the water self-diffusion flux26. In this work it was also observed that the non-field (6,6) flux results were much higher than those reported by Zhu and Schulten14, even though the same potential/force field was used in both studies. We suggested that this difference may be due to the array of parallel nanotubes employed by Zhu and Schulten which could cause friction (via long-range electrostatics) between the single files of water. This was confirmed in  where it was demonstrated that in these CNT-arrays there is a significant friction effect that retards movement of the single files of water.
In this chapter, we report further studies of non-equilibrium MD simulations of water permeation within (6,6) SWCNT’s across a lipid bilayer under the influence of both static and alternating electric fields, taking advantage of GPU-accelerated MD28 to extend the simulation timescales to longer than was previously possible26,27, while also investigating single-molecule translational and rotational frequency modes, collective properties and friction in tube-arrays in greater detail, allied with straightforward, mechanistic models to interpret our results.
2. Simulation methodology
All of the GPU-accelerated MD simulations23-25 were carried out using the NAMDv2.728,29 program suite, in conjunction with the CHARMM27 parameter set30 and the TIP3P model for water31. Initially, a periodically replicated single-nanotube system was studied, comprising a (6, 6) armchair SWCNT, with a diameter of 8.20 Å and approximate length of 36.9 Å, embedded in a 1-palmytoil-2-oleoyl-sn-glycero-3-phosphatidylcholine (POPC) membrane patch. Layers of water of around 10Å were added in the +
To model the arrays of CNTs, three-dimensional-periodic cells consisting of four (6, 6) SWCNT arrays were constructed with separations of 15, 20 and 25 Å between the geometric central axes of each (see Fig. 1), and inserted along the
The Particle Mesh Ewald33 method was utilised to treat long-range electrostatics to within a relative tolerance of 1x10-6 with a cut-off distance of 12 Å applied to real-space Ewald interactions; it was found that a real-space screening coefficient of 0.258 Å-1 and (three-dimensional Fourier transform) grid spacings in the 0.9 to 1 Å range in each Cartesian direction led to this desired accuracy level, along with close-to-optimal execution speeds. A cut-off distance of 12 Å was also applied to van der Waals interactions, using a smooth switching function between 10 and 12 Å. Using the r-RESPA method34, multiple time steps were applied to MD, with 1 fs for bonded interactions, 1 fs for short-range non-bonded interactions (including real-space Ewald interactions) and 2 fs for the full (real- and reciprocal-space) electrostatics evaluation. All production runs were carried out using Langevin dynamics35 coupled to an NVT reservoir set at 300 K with a damping coefficient of 1 ps-1, representing rather mild coupling which affects diffusion little. The SHAKE36 algorithm was used to constrain bond lengths to all hydrogen atoms.
Prior to production, an NPT relaxation run was carried out for 1 ns; a Nosé-Hoover approach was applied37 using anisotropic cell variation, with Langevin dynamics for piston fluctuation control38; here, the set points were 1 atm and 300 K. During this NPT relaxation, the nanotube atoms were fixed, while the lipid bilayer tails and water molecules were allowed to move freely. Following this, the CNT and lipid-head atoms were released slowly, keeping a harmonic constraint thereon with a strength of 1 kcal/mol Å². This strength coefficient was reduced progressively to 1.5, 1.0 and 0.5 kcal/mol Å² every 100 ps for 0.4 ns. After this initial relaxation, during which the pressure and volume had stabilised, MD simulations of each of the (6, 6) CNT-systems (single-tube and array) were carried out for 50 ns in the NVT ensemble (in the absence and presence of external fields –
Use of GPU-accelerated MD implemented in NAMD v2.728,29 on an NVIDIA double-precision C2050 platform led to a substantial acceleration (approximately 2.5-fold) vis-à-vis the equivalent number of CPU-cores on quad-core Intel Xeon Harpertown-Stoakley nodes connected via sub-microsecond-latency Infiniband. The acceleration was primarily, though not only, in the evaluation of intermediate-range non-bonded van der Waals and electrostatics interactions, where the computational overhead was reduced approximately six-fold relative to that required with the equivalent number of CPU-cores in a CPU-only approach (on the abovementioned Intel Xeon chips).
Static electric fields (EFs) were applied in the non-equilibrium MD simulations with an intensity of 0.0065 V/Å in the +
Non-equilibrium MD simulations were also carried out in the presence of alternating electric fields39 along the same
No constraints were applied to the membrane during MD in either the absence or presence of static and alternating electric fields, while a light 0.5 kcal/mol Å2 harmonic constraint was applied to the SWCNT atoms.
3. Modelling diffusion and permeation of water in single-file channels
The transport of water through single file channels embedded in lipid bilayers may be analysed from two perspectives (a) self-diffusion of water (e.g. labeled isotopically for direct physical measurement or, for simulation purposes, simply assigned a ‘colour’) and (b) water permeation under the influence of a pressure gradient (osmotic or otherwise). The two parameters which characterise these transport processes are the diffusive and osmotic permeabilities,
where is the tracer (or colour) concentration difference across the membrane,
In  and references cited therein (see also [12-15]) the diffusive permeability has been investigated using a continuous time random walk model (CTRWM) and, in the absence of a chemical potential difference
3.1. The diffusive permeability
To compute the diffusive permeability,
At the pore mouth a simple reversible rate process takes place for the molecules to either enter the pore or return into the reservoir on the left and the net rate is given by
where the forward and reverse rate constants
The steady-state flow along the pore may also be defined using an intra-pore diffusion coefficient,
Also, for the right hand side of the system it follows from Eqn. 7 that
Eqn. 8 simply states that the permeability is governed by a sum of mass transfer resistance terms due to bulk diffusion, incorporation into the pore mouth and intra-pore diffusion. Normally, one sets
The overall permeability can be simply determined from the MD simulations and by selecting different positions for the ‘black’ and ‘white’ reservoirs, one can eliminate the bulk diffusion
In the development of the CTRW model reported by Berezhkovskii and Hummer12, there is the implicit assumption that
In our model, we employ the physically reasonable assumption that the intra-pore hopping rate
which is the result provided by the CTRW model in . This analysis demonstrates the level of the assumptions involved in the CTRWM approach. As noted above, in a more general setting, the in-pore forward hopping frequency is not the same as at the pore mouth. Furthermore, the more general model of Eqn. 11 applies to wide (i.e. multi-pass) pore conditions of arbitrary cross-sectional area and it is also simple to extend Eqn. 11 to include rate processes at other junctions (e.g. local charged sites) within the pore and to relate the rate constants to the energy barriers associated with the particle movement along the channel.
3.2. The osmotic permeability, collective properties and friction in pore-arrays
In general, the flux of water molecules across a membrane subject to a chemical potential driving force, Δ
which is equivalent to Eqn. (4). Here, is the molecular density of water while
whereis the water concentration within the membrane (which is related to via the partition coefficient44
where is the velocity of water molecule
For a membrane containing an array of
where the subscript denotes the limit of infinite separation of all pores, and is the permeability for a single, isolated pore.
In order to explore further the underlying nature of the collective pore-velocity correlations of Eqn. 16 which govern the distance-dependent friction observed in our previous work27 it is appropriate to consider a simplified model of friction. For a given collective velocity within the pores, , the friction may be approximated by a term due to isolated pores (the ‘’ case), as well as an additional, distance-dependent term (a finite separation ‘f’ case), which will be zero in the case of the limit of infinite separation:
the solution of which is
Here, the isolated-case friction term may be considered as a relaxation time which might possibly be expected in the decay of the ACF of the isolated-case collective velocity, while that of would be additional friction that might be manifest in the decay of the TCF of Eqn. 16, vis-à-vis that of the isolated case. At its simplest, one might model the distance-dependent friction term, , to be inversely proportional to separation
so that the permeability may be expressed in terms of the frictional terms:
leading to an effective slip-diffusion coefficient,
4. Results and discussion
4.1. Static and alternating electric field effects and water self-diffusion
The permeation rates for water molecules through the SWCNT (in molecules per ns) embedded within the lipid layer at steady-state, in the absence and presence of axially-applied static and alternating electric fields, are provided in Fig. 3. Although the rate in the static field is marginally higher (notwithstanding overlapping error bars vis-à-vis the zero-field case), the alternating field results show no discernible trend when the electric field frequencies are incremented. It is important to note that each point represents the average fluxes of three independent simulations. Statistically, with the possible exception of the static-field case, we observe that there is no correlation between the field frequencies and the computed water permeation rates. This is in essential agreement with our previous study27 although longer MD timescales from GPU-acceleration has led to reduced error bars in this work. In our previous work26,27 monotonic decreases of water permeation rates with field frequencies were found for axial field application to (5,5) CNT-systems towards the zero-field rate at higher frequency, but not to (6,6) CNT-cases (or larger CNTs). For (5,5) CNT-systems, it was found that axial static electric fields enhanced water self-diffusion across them due to a decrease in the fluctuations of the number of water molecules within the pores, arising from the water molecules’ dipole orientation with the field26. Although it is difficult to reach a conclusion for the (more statistically accurate vis-à-vis ) slightly larger static-field rate in the case of this work for (6,6) CNTs, due to overlapping error bars with respect to the zero-field result, the single-file nature of water transport in (6,6) CNTs, similar to the (5,5) case, suggests that a more minor decrease in the fluctuations of the number of water molecules inside is at play in the static field. This was observed in this study with results slightly lower than previous (6,6) data26,27 although error bars for the size of the fluctuations were again overlapping vis-à-vis the zero-field case, so no firm statistical conclusion can be drawn in that regard. It would be necessary to carry out a greater number of longer runs to investigate this further, somewhat prohibitive even with ample GPU-accelerated computing resources.
The water dipole orientations within the pore (see Fig. 4), illustrate more clearly the influence of the external field on the probability distributions of the cosine of the angle of the dipole with the positive
From consideration of the fluxes and dipolar orientations, it is clearly evident that the results of the low-frequency electric field simulations are more similar to the static-field case, whereas high-frequency fields resemble zero field conditions. These results support the finding that the effects of low-frequency electric fields should be closer to those of static electric fields (if any), because the slower time variation of the field (
Considering single-molecule translational and librational modes in the isolated CNT system, the power spectra of the oxygen- and hydrogen-atom VACFs are shown in Figs. 5 and 6, respectively, under zero-field conditions, along with corresponding results for bulk water. The spectra are decomposed in terms of
4.2. Water permeation through nanotube-arrays
Steady-state water diffusion rates,
The results summarized here are essentially in agreement with previously reported values25,26 albeit with reduced error bars due to the longer timescales rendered possible from GPU-acceleration. There is an incremental increase in
Estimates of the hopping rate constant,
|15||1.69 (0.18)||0.050||1.94 (0.07)||57.40||0.896||1.74|
|20||2.12 (0.28)||0.063||2.00 (0.06)||59.17||1.156||1.85|
|25||2.32 (0.13)||0.069||2.14 (0.11)||63.31||1.268||1.93|
|2.47 (0.29)||0.074||2.38 (0.12)||70.41||1.368||2.18|
To demonstrate the existence of significant frictional effects between the aligned water files in  we estimated an effective viscosity along the axis of the water confined inside the CNTs. In summary, this was done by use of the appropriate Einstein relation23 applied to tracking chains of six water molecules initially nearest the centre of each of the four nanotubes. The four hexamers of water so defined were treated as a fluid of four ‘particles’ and their motion was monitored until one or more of the water molecules in any given hexamer left its respective nanotube. The hexamers were then redefined and sampling continued until the next disruption in one of the hexamer chains. Averaging over each such occurrence in all four tubes in the zero-field simulations then enabled the determination of the effective viscosity of the files using:
The long-time slopes of the mean square displacement of the dynamical variable as functions of time (see Fig. 7) for each nanotube separation provided the effective viscosities 7.5 ± 0.45 (15 Å), 5.5 ± 0.4 (20 Å) and 4.4 ± 0.32 (25 Å) × 10-3 Nsm-2. These are agreement with our previous study27 although the standard deviations are reduced due to longer simulations. The significantly larger effective viscosity the closer the tube-tube separation confirms the strong interaction induced by the surroundings on the water molecules confined within each individual CNT pore. The dependence of
To investigate friction further, collective, axial VACFs of all of the water molecules in each individual CNT along the
The osmotic permeabilities were evaluated in two independent ways, and the results are provided in the last two columns of Table II. Integration of the Green-Kubo-type representation of the viscous slip total velocity correlation function (TCF) (Eqn. 16) (see [42,43] and also [23, 45,46] for additional details on Green-Kubo analysis) leads to convergence within 10-15 ps and the results for the osmotic permeability are provided in the second-last column of Table II. These results are the exact values for the osmotic permeabilities. The larger values obtained from the approximate friction analysis (Eqn. 19) are to be expected given that the friction model does not take into account negative backscattering parts of the velocity correlation function clearly observed in Fig. 9a. Nonetheless, the general trend with nanotube separation is clearly evident. Also note that the ratio of the osmotic and diffusion permeabilities
Single-file diffusion results through (6,6) SWCNTs, embedded in a lipid bilayer, under the influence of static and time-varying electric fields, have been reported. No substantial influence of alternating electric fields on water flux through SWCNTs of larger diameters was observed, although dipole orientations were affected markedly, with higher-frequency distributions resembling those of the zero-field case. In the case of static fields, however, the water flux appears somewhat higher than that of the zero-field case, owing possibly to a similar reduction in water-number fluctuations therein, as found in previous work for (5,5) CNTs26,27; an overlapping of error bars vis-à-vis the zero-field case, however, renders any firm statistical conclusion elusive in this respect. The (RMS) field intensities of 0.0065 V/Å are weak compared to previous molecular simulation studies employing static or time-varying electric fields39 and also weak relative to fields generated by charge distributions14, and are substantially more realistic with respect to experiment.
When arrays of CNTs are employed in the simulations the interpretation on the transport properties of the fluid within the individual nanotubes must take cognisance of inter-tube interactions. The (6,6) CNT-array simulations reported in this work have demonstrated the presence of a friction effect due to long-range electrostatic interactions between the moving water single-files within the pores, with better statistics than was possible in previous work.14 It was shown that the generated friction is inversely proportional to tube separation, reaching values closer to the single-pore simulations (which should be considered as an infinitely-separated system) at the largest separation between the tubes, and estimation of an effective viscosity of the water molecules confined within the nanotubes confirmed this, along with consideration of the collective total velocity correlation functions (TCFs) for the water files within the arrays of CNTs.
A more general approach42,43 to the concept of the osmotic permeability based on linear response theory leads to the recognition that pf is simply related to the slip-diffusion (or viscous slip) coefficient of a fluid undergoing pressure-driven transport in a nanopore structure. The results reported in this work also support the approximate relationship pf ~ (N + 1)pd suggested by Tajkhorshid et al 10 where N is the average number of water molecules within a single carbon nanotube.