Electric-Field and Friction Effects on Carbon Nanotube-Assisted Water Self-Diffusion Across Lipid Membranes

1. Introduction

The accelerated appearance of applications and theoretical studies of carbon nanotubes (CNTs)¹ in the literature over the last two decades reveals their broad attraction in fields as diverse as nanoelectronics, bionanotechnology and solar energy conversion nanodevices^2-5. In particular, they provide the basis for simplified models of nanofluidic machines and membrane channels such as aquaporins^6-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 CNTs^10-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 Hummer^12,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 CNTs^14-16 and in model pores containing charge distributions which are known to be ubiquitous for biological channels^17,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 aquaporins¹⁹, 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 attention^20-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 flux²⁶. In this work it was also observed that the non-field (6,6) flux results were much higher than those reported by Zhu and Schulten¹⁴, 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 [27] 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 MD²⁸ to extend the simulation timescales to longer than was previously possible^26,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 simulations^23-25 were carried out using the NAMDv2.7^28,29 program suite, in conjunction with the CHARMM27 parameter set³⁰ and the TIP3P model for water³¹. 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 +z and -z-direction to solvate the system. Each atom of the SWCNT was modelled as a neutral sp² aromatic carbon with uncapped ends (in a similar way as shown in Fig. 1, although this is for arrays of four CNTs). For all of these systems, the plane of the POPC membrane was placed in the x-y plane, and the axis normal to the membrane was along the z-direction, with the CNT oriented along this z-axis (perpendicular to the plane of the membrane – cf. Fig. 1). The final size of the (three-dimensional-periodic) simulation box was 34.84 x 34 x 66.29 ų, with total atom numbers of 8,415 (1337 water molecules and 30 phospholipids).

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 z-axis of a POPC membrane, again solvated with layers of water molecules of around 10 Å in thickness in the +z and -z-directions. In each of these cases, the membrane was also placed in the x-y plane, with the principal axis of the SWCNT oriented perpendicular to the membrane. The final sizes for the simulation boxes in the array-systems were 45.75 x 60.32 x 66.68 ų, 49.36 x 49.21 x 71.27 ų and 63.10 x 59.99 x 66.27 ų for respective separations of 15, 20 and 25 Å, with total atom numbers of 18,979 (3045 water molecules and 62 phospholipids), 17,901 (3043 water molecules and 54 phospholipids) and 26,309 (4327 water molecules and 88 phospholipids), respectively. All molecular systems were constructed using the VMDv1.86³² program.

Figure 1.

a), (c) and (e) Cross-section of the periodic cells simulated for the (6,6) CNT-arrays separated by 15, 20 and 25 Å, respectively. (b), (d) and (f) Plan view (along the z-axis, with the x-y plane in that of the page) of the (6,6) CNT-array systems separated by 15, 20 and 25 Å, respectively.

The Particle Mesh Ewald³³ 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 method³⁴, 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 dynamics³⁵ 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 SHAKE³⁶ 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 applied³⁷ using anisotropic cell variation, with Langevin dynamics for piston fluctuation control³⁸; 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 – vide infra). Three replicates were performed for each system in each case (zero- and in-field).

Use of GPU-accelerated MD implemented in NAMD v2.7^28,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 +z-direction (i.e., along the CNT-axis). The EF exerts a force over partial atomic charges ia:

$${\text{f}}_{\text{ia}}={\text{ q}}_{\text{ia}}\text{E}$$ E1

Non-equilibrium MD simulations were also carried out in the presence of alternating electric fields³⁹ along the same z-axis (i.e., along the axes of the nanotubes). In this case, the RMS intensity was set to 0.0065 V/Å, the same as that of the peak frequency for the static EF, but the field does change direction every period between the +z and –z axis, ipso facto ^40,41. Three different frequencies were used: 2.45, 50 and 100 GHz. Here, the force exerted over the partial atomic charges is given by:

$${\text{f}}_{\text{ia}}{=}_{ }{\text{q}}_{\text{ia}}{\text{E}}_0\text{cos}\left(\omega \text{t}\right)$$ E2

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 Ų 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, p _d and p _f , respectively (see [10] for a review of details). These are defined through the flux expressions:

$$j_d=p_d\Delta c_{tr}$$ E3

and

$$j_f=\frac{1}{k_{B}T}{p}_f\Delta P$$ E4

where $$\Delta c_{tr}$$ is the tracer (or colour) concentration difference across the membrane, P is the trans-membrane liquid pressure difference (due, under normal circumstances, to an osmotic bulk solute concentration difference), T is absolute temperature and k _B is Boltzmann’s constant.

In [10] 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 Δμ or an external force, the osmotic permeability was computed by equilibrium MD using a collective diffusion model. In [26, 27], a more general model was developed for the diffusive permeability (which is not restricted in application to single-file pores) and this more general approach is described below. In addition, the collective transport permeability is analysed in this chapter in greater detail using the generalized linear transport theory for slip flow in pores developed in earlier work^42,43.

3.1. The diffusive permeability

To compute the diffusive permeability, p _d , of water in a strongly inhomogeneous system constituting two bulk liquid reservoirs separated by a lipid membrane with embedded carbon nanotubes, the simplest and statistically most accurate approach is to employ a colour diffusion technique. In this approach the flux of ‘black’ water molecules diffusing from left to right in Fig. 2 is monitored during the course of a long molecular dynamics simulation (a similar analysis applies to ‘white’ water molecules diffusing in the opposite direction). If the cross-sectional areas in the bulk and in the pore are A and A _p , respectively, and L is the length of the mixing zone between the black reservoir on the left of Fig. 2 and the pore mouth, then at steady state the net flow (particles per unit time) between the black reservoir and the outer perimeter of the CNT is given by

$$j_d=\frac{DA}{L}(c_B(1)-c_B(2))$$ E5

where c _B (1) is the concentration of black water molecules in the black reservoir and c _B (2) is the concentration of black water molecules just outside the CNT pore mouth.

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

$$r_{\text{2}-\text{2p}}=k_f{c}_B(2)-k_r{c}_{Bp}(2)$$ E6

where the forward and reverse rate constants k _f and k _r are in units of sec^-1 and c _Bp (2) is the concentration of black water molecules just inside the CNT pore adjacent to the pore mouth. The net rate in Eqn. (6) is a volumetric measure (these rate constants will be determined to a great extent by the hydrogen bond structure and dynamics). Multiplying by a volume V _i (an interfacial volume corresponding to A _p δ with δ = length scale of the order of the molecular dimensions representing the scale of the concentration inhomogenities) then for steady state conditions we have

$$\begin{array}{l}{V}_i{r}_{\text{2}-\text{2p}}=j_d\\ =V_i(k_f{c}_B(2)-k_r{c}_{Bp}(2))\\ =V_i{k}_f(c_B(2)-\frac{1}{K_p}{c}_{Bp}(2))\\ =\frac{DA}{L}(c_B(1)-c_B(2))\end{array}$$ E7

Figure 2.

Schematic of the ‘black’ and ‘white’ water reservoirs on either side of the CNT, with the CNT oriented along the z-axis, normal to the membrane patch. The concentration of the black water molecules is specified at the edge of the black reservoir (c _B (1)), at the pore mouth (c _B (2)), just inside the pore (c _Bp (2)), with analogous concentrations on the other side of the simulation cell (the positive z-coordinate half), i.e., c _Bp (3), c _B (3) and c _B (4). L _c was 3.5 nm in all of the systems simulated.

Eliminating c _B (2) then

$$j_d=\frac{p_1}{K_p}(K_p{c}_B(1)-c_{Bp}(2))$$ E8

where K _p is defined as k _f /k _r and p ₁ is the partial permeability for diffusive transport from the black reservoir to just inside the CNT:

$$\frac{\text{1}}{p_{\text{1}}}=\frac{L}{DA}+\frac{\text{1}}{k_f{V}_i}$$ E9

The steady-state flow along the pore may also be defined using an intra-pore diffusion coefficient, D _p (which is considered to exist in the same sense as proposed in [10, 12-14] but which has a deeper meaning than indicated by the simple CTRWM)

$$j_d=\frac{D_p{A}_p}{L_p}(c_{Bp}(2)-c_{Bp}(3))$$ E10

where L _p is the length of the pore and c _Bp (3) is the concentration of black water molecules just inside the CNT pore adjacent to the pore mouth near the white particle reservoir.

Also, for the right hand side of the system it follows from Eqn. 7 that

$$j_d=\frac{p_1}{K_p}(c_{Bp}(3)-K_p{c}_{Bp}(4))$$ E11

where c _Bp (4) is the concentration of black water molecules in the white particle reservoir.

Now combining Eqns. 8, 9 and 10 eliminating the intermediate concentration terms gives

$$j_d=p_d(c_B(1)-c_B(4))$$ E12

where

$$\frac{1}{p_d}=\frac{2L}{DA}+\frac{2}{k_f{V}_i}+\frac{L_p}{D_p{K}_p{A}_p}$$ E13

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 c _B (4) = 0 (in the simulations the ‘black’ water molecules entering the ‘white’ reservoir are simply switched in colour to ‘white’) and c _B (1) = molecular/molar density of water in the bulk i.e. c _H2O = 33.5 nm^-3 or 55.6 mol.litre^-1. Therefore

$$j_d=p_d{c}_{H_{2}O}$$ E14

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 D/L terms.

In the development of the CTRW model reported by Berezhkovskii and Hummer¹², there is the implicit assumption that L is sufficiently small to minimise bulk-phase diffusional resistance while retaining a sub-nanometre range close to the pore mouth (δ ~ a) to ensure that entry/exit kinetics are fully taken into consideration. This condition with the additional assumption that the single file is considered to be densely packed leads to

$$L_p=Na;\delta =a$$ E15

where a is a length scale characteristic of the size of the water molecules. Furthermore, the diffusion coefficient within the pore is given by the CTRW model as

$$D_p=\frac{a^{\text{2}}}{\text{2τ}}=\frac{a^{\text{2}}}{\text{2}}{k}_r$$ E16

In our model, we employ the physically reasonable assumption that the intra-pore hopping rate k is equal to k _r , i.e., the reverse rate constant in Eqn. 5. Within the pore, the forward and reverse hopping rate constants are equal. At the pore mouth, however, we distinguish between the forward (pore entry) and the reverse (pore exit) rate constants. For the CNT’s under consideration here, one may anticipate k _r > k _f or K _p < 1 and if we assume the molecular volume $$v_w=a{A}_p$$ , then we find

$$\text{Permeation rate}=\frac{k_{f}a{A}_p}{\text{2}}\frac{\text{1}}{N+\text{1}}{c}_{H2O}=\frac{v_w{k}_o}{N+\text{1}}{c}_{H2O}$$ E17

which is the result provided by the CTRW model in [12]. 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 in turn originates from an osmotic pressure gradient, may be expressed as

$$j_f=\frac{c_{H_{2}O}}{k_{B}T}{p}_f\Delta \mu$$ E18

which is equivalent to Eqn. (4). Here, $$c_{H_{2}O}$$ is the molecular density of water while p _f is osmotic permeability as defined, inter alia, in [10, 14-16]. Using results from linear response theory^42,43 the permeability of membranes investigated in this work may be expressed as

$$\begin{array}{l}{c}_{H_{2}O}{p}_f=\frac{c_{H_{2}O,M}{D}_M^{(z)}A}{L_p}={\left(\frac{N_T}{L_p}\right)}^2{\displaystyle \underset{0}{\overset{\infty }{\int }}{u}_z(t)u_z(0)}dt\\ ={\left(\frac{N_T}{L_p}\right)}^2\underset{t\to \infty }{\lim }\frac{d{(z(t)-z(0))}^2}{dt}\end{array}$$ E19

where $$c_{H_{2}O,M}$$ is the water concentration within the membrane (which is related to $$c_{H_{2}O}$$ via the partition coefficient⁴⁴ K between the bulk liquid phase on both sides of the membrane and the membrane pores), N _T is the total number of water molecules within the membrane at any given time and $$D_M^{(z)}$$ is the axial slip-diffusion (or viscous slip) coefficient of the water molecules within the CNT-pores of the membrane. The function in the integral of Eqn. 12 is the autocorrelation function (ACF) of the axial, collective centre-of-mass velocity of all of the water molecules present in the pores of the membrane, with the collective (number averaged) velocity given by

$$u_z(t)=\frac{1}{N_T}{\displaystyle \sum _{i=1}^{N_T}{v}_{i,z}(t)}$$ E20

where $$v_{i,z}\left(t\right)$$ is the velocity of water molecule i along the z-axis of the CNT-pores. In Eqn. 12, z(t) is the corresponding collective, axial centre-of-mass position of the water molecules in the pore at any point in time:

$$z(t)=\frac{1}{N_T}{\displaystyle \sum _{i=1}^{N_T}{z}_i(t)}$$ E21

For a membrane containing an array of N _p pores, each of which is accessible to water to the same extent, as studied in this work, we note that Eqn. 12 may be expressed as:

$$\begin{array}{l}{c}_{H_{2}O}{p}_f={\left(\frac{1}{L_p}\right)}^2{\displaystyle \underset{0}{\overset{\infty }{\int }}{\displaystyle \sum _{j=1}^{N_p}{\displaystyle \sum _{i=1}^{N(j)}{v}_{i,z}(t){\displaystyle \sum _{k=1}^{N_p}{\displaystyle \sum _{i=1}^{N(k)}{v}_{i,z}(0)}}}}}dt\\ ={\left(\frac{1}{L_p}\right)}^2{\displaystyle \underset{0}{\overset{\infty }{\int }}{\displaystyle \sum _{j=1}^{N_p}N(j)u_z^{(j)}(t){\displaystyle \sum _{k=1}^{N_p}N(k)u_z^{(k)}(0)}}}dt\end{array}$$ E22

where N(j) is the number of water molecules in pore j. The correlation function in Eqn. 16 may be considered a total one (TCF), in that it applies to all water molecules in the membrane (i.e. present in all of the CNT-pores). Considering Eqn. 16, it is clear that this may be written in terms of the ACFs of the collective velocities in individual pores as well as cross-correlation functions (CCFs) of collective velocities between particular pores. Eqn. 16 demonstrates the possibility of significant cross-interactions between the individual pores in the array, unless they are well separated. Should this latter condition apply, then Eqn. 16 simplifies to:

$$\begin{array}{l}{c}_{H_{2}O}{p}_f|{}_{\infty }={\left(\frac{1}{L_p}\right)}^2{\displaystyle \underset{0}{\overset{\infty }{\int }}{\displaystyle \sum _{j=1}^{N_p}{N}^2}\left(j\right)u_z^{\left(j\right)}\left(t\right)u_z^{\left(j\right)}\left(0\right)}dt\\ =N_p{c}_{H_{2}O}{p}_f^{\ast }\end{array}$$ E23

where the subscript $$\infty$$ denotes the limit of infinite separation of all pores, and $$p_f^{*}$$ 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 work²⁷ it is appropriate to consider a simplified model of friction. For a given collective velocity within the pores, $$u_z$$ , the friction may be approximated by a term due to isolated pores (the ‘ $$\infty$$ ’ 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:

$$\frac{\text{d}{u}_z}{\text{d}t}=-\left(\text{1}/{\tau }_{\infty }+1/{\tau }_f\right)u_z$$ E24

the solution of which is

$$u_z\left(t\right)=u_z\left(0\right)\exp \left(-\left(\text{1}/{\tau }_{\infty }+1/{\tau }_f\right)t\right)$$ E25

Here, the isolated-case $$\text{1}/{\tau }_{\infty }$$ 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 $$\text{1}/{\tau }_f$$ 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, $$\text{1}/{\tau }_f$$ , to be inversely proportional to separation d, i.e. $$\text{1}/{\tau }_f$$ ~ 1/d. Exploring this collective-velocity correlation further for the CNT-array:

$$\begin{array}{l}{\displaystyle {\int }_0^{\infty }{\left({\displaystyle \sum _{j=1}^{N_p}{u}_{j,z}(0)}\right)}^2\exp \left(-\left(1/{\tau }_{\infty }+1/{\tau }_f\right)t\right)dt}={\left({\displaystyle \sum _{j=1}^{N_p}{u}_{j,z}(0)}\right)}^2\frac{1}{1/{\tau }_{\infty }+1/{\tau }_f}\\ =\frac{1}{N_T}\frac{k_{B}T}{m}\frac{1}{1/{\tau }_{\infty }+1/{\tau }_f}where1/{\tau }_f~\text{ 1}/d\end{array}$$ E26

so that the permeability may be expressed in terms of the frictional terms:

$$\begin{array}{l}\Rightarrow c_{H_{2}O}{p}_f=\frac{N_T}{L^2}\frac{k_{B}T}{m}\frac{1}{1/{\tau }_{\infty }+1/{\tau }_f}\\ =c_{H_{2}O,M}\frac{A}{L}\frac{k_{B}T}{m}\frac{1}{1/{\tau }_{\infty }+1/{\tau }_f}\end{array}$$ E27

leading to an effective slip-diffusion coefficient, D _eff :

$$p_f=\frac{K{D}_{eff}A}{L_p}\text{ where }\frac{1}{D_{\text{eff}}}=\frac{m}{k_{B}T}\left(\text{1}/{\tau }_{\infty }+1/{\tau }_f\right)$$ E28

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 study²⁷ although longer MD timescales from GPU-acceleration has led to reduced error bars in this work. In our previous work^26,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 field²⁶. Although it is difficult to reach a conclusion for the (more statistically accurate vis-à-vis [27]) 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) data^26,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.

Figure 3.

Steady-state permeation rates under the influence of static and time-varying electrical fields. The dashed black line represents zero-field conditions with the corresponding error bars shown at either end thereof.

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 z-axis, α. In agreement with our earlier results, the statistically more accurate results reported here illustrate that at high frequencies the distributions become more symmetric with peaks at -1 and 1, in very good agreement to the zero-field case. In the static field case the distribution is shifted almost entirely to the left side of the diagram.

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 (i.e., longer period) relative to permeation and loading timescales allows the field to influence the dynamics of these events. On the other hand, high-frequency electric fields do not readily permit enhancements in permeation or alterations in dipole alignment vis-à-vis the zero-field case, due to their faster oscillation and shorter period relative to permeation timescales in (6,6) CNTs.

Figure 4.

Normalised probability distributions of the cosine of the angle which the dipole vector of the water molecules makes with the positive z-axis, α.

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 x-y and z-directions, and a translational ‘bumping’ mode vis-à-vis the CNT walls is clearly visible in the x-y plane at around 35 cm^-1 (circa 0.8-0.9 ps) in the oxygen-atom spectra, along with corresponding librational modes of 35 and 160 cm^-1 (around 0.18 and 0.8-0.9 ps). The axial-spectra, although somewhat more bulk-like, also exhibit important differences vis-à-vis the bulk state: notably a greater translational ‘rattling’ at 260 cm^-1 (around 0.12 ps), and a marked reduction in higher-frequency librational motion, in accord with retarded rotational motion (evident in the preferred dipole distributions of Fig. 4).

Figure 5.

Power spectra of single-water molecule VACFs for oxygen atoms in the isolated CNT under zero-field conditions, showing (a) the total spectrum along with corresponding results for bulk water, and (b) the x-y and z-direction spectra in the CNT.

Figure 6.

Power spectra of single-water molecule VACFs for hydrogen atoms in the isolated CNT under zero-field conditions, showing (a) the total spectrum along with corresponding results for bulk water, and (b) the x-y and z-direction spectra in the CNT.

4.2. Water permeation through nanotube-arrays

Steady-state water diffusion rates, j _d , and hopping rates, k, were computed for the (6,6) SWCNT array simulations and are reported in Table I. The intra-pore water self-diffusion coefficients, D _p, provided in this Table were computed independently as described in [26] using the Einstein relation for self-diffusion (a special case of Eqn 12)

$$D_p=\underset{t\to \infty }{\lim }\frac{1}{2}\frac{d{(\Delta z)}^2}{dt}$$ E29

The results summarized here are essentially in agreement with previously reported values^25,26 albeit with reduced error bars due to the longer timescales rendered possible from GPU-acceleration. There is an incremental increase in j _d when the separation between the pores is augmented, which is a direct consequence of the incremental increase in D _p with the separation of the nanotubes for separations greater than 15 Å. The lower than anticipated diffusive flux for the 15 Å close packed structure when comparing the magnitude of j _d relative to D _p arises in part due to the distance-dependent friction effect²⁷ that retards water movement within the pores as demonstrated by the osmotic permeability results (the right-hand column of Table I) but also reflects a greater degree of interference between the pores at the pore mouths, namely a lowering of the forward (pore entry) rate constant k _f in the second term in Eqn. 11.

Estimates of the hopping rate constant, k, suggested by the CTRW model were also computed to compare our results with those reported in refs. [12-14], since, as previously shown^26,27 our results for the zero-field (6,6) CNT were in quite good agreement with those reported by Berezhkovskii and Hummer^12,13 but differed with the ones reported by Zhu and Schulten¹⁴, even though both our results and the ones previously mentioned are based on the TIP3P/CHARMM potential. The hopping rate result, k = 57.4 ns^-1, reported here for simulations involving the 15 Å separation is more than a factor of 2 higher than the corresponding value reported by Zhu and Shulten (k = 26.9 ns^-1). As discussed in [27] only four tubes were employed in our simulations (embedded in a non-diffusing phospholipid membrane) while a close packed periodically imaged hexagonal array of twelve tubes was employed in the work conducted by Zhu and Schulten.¹⁴ The much lower k observed in ref. 14 is explained by an increase in the friction effect generated by the twelve single files of water interacting via long range electrostatics, rather than the four nanotubes employed in our work.

Separation (Å)	jd (ns^-1)	pd (nm^3.ns^-1)	Dp (nm^2 ns^-1)	^ak (ns^-1)	^bpf (nm^3.ns^-1)	^c p f F M(nm^3.ns^-1)
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

Table 1.

The hopping rates defined by the CTRW model, k (≡ k _r ) = 2D _p /a ² , where a = 0.26 nm is the projected length for a given water molecule along the pore axis (for more details see [12]).^b The osmotic permeability predicted by integration of the viscous slip total velocity-CF (TCF) (Eqn. 16).^c The osmotic permeability predicted by the simplified friction model (Eqn. 19).Steady-state tracer permeation rates, tracer intrapore diffusion coefficients, hopping rates and osmotic permeabilities for the CNT arrays.

To demonstrate the existence of significant frictional effects between the aligned water files in [27] 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 relation²³ 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:

$$\eta =\frac{V}{2k_{B}T}\underset{t\to \infty }{\lim }\frac{d}{dt}{\left({\text{P}}_{zz}(t)-{\text{P}}_{zz}(0)\right)}^2{\text{ where P}}_{zz}=\frac{1}{V}{\displaystyle \sum _{i=1}^{N_T}{r}_{iz}{p}_{iz}}$$ E30

The long-time slopes of the mean square displacement of the dynamical variable $${\text{P}}_{zz}$$ 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 study²⁷ 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 η on the separation distance between the nanotubes is more clearly demonstrated in Fig. 8 where the fit to a simple 1/d dependence is shown to be excellent.

Figure 7.

Mean square displacement of the dynamical variable $$P_{\alpha \beta }$$ (with α = z and β = z) for the four hexamer water files as functions of time and for three different nanotube-array separations.

Figure 8.

Plots of the CNT array effective viscosity coefficient (from Eqn. 20), filled circles, and the 1/τ_f distance-dependent friction parameter (inferred from the fitted peak-decay lines in Fig.9b), filled squares, versus the inverse tube-separations (i.e., 15, 20, 25 Å).

To investigate friction further, collective, axial VACFs of all of the water molecules in each individual CNT along the z-axis are shown in Fig. 9a, including the isolated tube case. A log-plot of this in Fig. 9b shows an essentially negative-exponential decay, along the general trend of an $$A\exp (-t/\tau )\cos \omega t$$ motif. Fitting to the successive peaks in the log-plot yields an isolated tube frictional relaxation time of 0.16 ps. From the other slopes, the distance-dependent frictional contribution may be determined for each of the four-CNT arrays (inferred from the fitted peak lines), and these are detailed in Table II. The 1/τ_f distance-dependent friction parameter declines with increasing separation, and a plot of this versus the inverse separation as shown in Fig. 8, coupled with a regression line, shows an approximate 1/d relationship. Direct evaluation of the z-axis collective-velocity cross correlation function (VCCF) for each tube with the others in the three CNT-array cases confirms that the magnitude of the VCCF declines with increasing separation, and this is shown in Fig. 10.

Figure 9.

a) Normalised collective TCFs of all of the water molecules in each CNT along the z-axis; this is averaged over each of the four CNTs in the case of the arrays. (b) Log-plot of the normalised collective TCFs, showing negative-exponential fits to the peaks.

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 p _f /p _d in all cases (with p _f given by the exact results in the second last column of Table I) is approximately 18 in magnitude. This corresponds reasonably with the estimate of 15 obtained from an approximation suggested in reference [10].

Figure 10.

Collective z-axis VCCFs of one pore with the other three in the CNT-array systems, averaged over all such cross-interactions and normalised vis-à-vis the four-averaged TCFs (shown in Fig. 9a).

Separation (Å)	τeff (ps)	τf (ps)^*
15	0.128	0.64
20	0.134	0.825
25	0.142	1.26
^∞	0.16	Infinite

Table 2.

Inferred from 1/τ_eff = 1/τ_eff, $$\infty$$ + 1/τ_f Table II. Fitted values for effective frictional relaxation times, τ _eff , and values of the resultant separation-friction relaxation time, τ _f . See plot of 1/τ _f versus 1/d (Fig. 8).

5. Conclusions

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) CNTs^26,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 fields³⁹ and also weak relative to fields generated by charge distributions¹⁴, 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.¹⁴ 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 approach^42,43 to the concept of the osmotic permeability based on linear response theory leads to the recognition that p_f 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 p_f ~ (N + 1)p_d suggested by Tajkhorshid et al ¹⁰ where N is the average number of water molecules within a single carbon nanotube.