The Motion of Ions Confined in a Molecular Channel

In order to understand the features governing the motion of ions in a molecular environment the migrational features of Na+ and Cl− ions in a molecular channel composed of stacked crown ether 6-CE-18 rings is followed using molecular dynamics, which shows that Na+ is subject to a much greater dynamic hindrance than the Cl− ion. The effects of the fluctuating electric fields of the atomic constituents in the channel on the motion of the migrants are investigated by clamping them so as to remove the fluctuations. The dynamic system is simulated both in vacuo and in water. For both it is found that the fluctuating electric fields of the channel and water atoms play a significantly greater role in the ion motions than do fluctuations in the ‘non-bonded’ interactions. The effect of temperature on the dynamics is investigated. Oscillatory trajectories are followed via the force and the potential energy profiles of the system over the timestep range of the molecular dynamics.


Introduction
Details of the behaviour of ions in small volumes is of importance not only to the understanding of functioning of the cell but also in investigations of particles which can migrate under the user's control through artificial channels in the design of novel nano-scale molecular devices [1,2]. Molecular dynamics of the latter systems, which include crown ether channels and molecular helices, have been performed in this laboratory [3]. In the latter investigations the migrations of the ions were 'active' by requiring an external action such as an applied electric field. Many natural migrations, however, are 'passive', i.e. their driving force is internal to the migration system, operating as a consequence of temperature, osmotic effects or dynamic conformational changes in the molecular structure of the channel. Because of its role in channel migration it would be valuable to ascertain the origins of the complete driving force. However as the complexities of such molecular systems entail considerable labour for their detailed elucidation, such studies may more profitably be conducted on models in which a migrant species enters a confined space such as exists in a molecular channel, a procedure which is followed in the current work. The motion of an ion in the channel, for example, is dependent on several 'active' and 'passive' factors such as ionic charge, electric field, differential concentrations of molecular species in the medium such as pH changes and osmosis, non-bonded interactions between migrant and channel species or transporter, and the environmental conditions of temperature and pressure. While most of these can be studied by molecular simulation and dynamics, another factor is the fluctuating motion of the (partially) charged atoms of the channel that are in the proximity of the migrant. Supposing that the inhibitions to the migrant's progress in the sterically confined space originates in its requirement to surmount the energy barriers that are associated with local energy wells, it would be of interest to examine the role played by the fluctuating motion and charges on atoms in the molecular environment that facilitate the ion's progress by a mechanism that is analogous to that producing the Brownian motion of pollen grains by the random buffeting of solvent molecules [4]. It would therefore be useful to pursue a quest leading to an understanding of the effects that govern the motion of ions through a solid. A dynamics simulation of actual biological ion channels would be excessively complex to provide the required fine details of the motions. An analogous molecular model of reduced complexity is therefore presented that however still retains sufficient features of the 'real' ones to permit the derivation of a useful analysis.

Molecular System
For this purpose a molecular confined space is considered in which 12 amino crown ether 6-CE-18 like rings (CE) depicted in Figure 1(a) which are stacked mutually parallel between two similar 'capping' CE rings which differ from the other rings by the substitution of a carboxylic acid COOH group as shown in Figure 1(b) (N atoms have been included in the rings since the NH groups cause lass steric crowding of H atoms than do CH 2 ). In the Figure C atoms are black, O is red and N is green.) Part (c) of the Figure shows a channel with a length of 31.26 Å between the centres of the terminal rings; the internal diameter is between 5.9 and 7.7 Å, depending on its definition. Figure 1 parts (c) and (d) also indicate the initial positions of the Na + (red) and Cl − ion migrants. The number of component CE rings (12) was chosen since the overall length which they confer on the channel corresponds to the thickness of several types of lipid membranes that hold ion channels in biological cells. The purpose of including the COOH groups in the terminal rings is twofold. Firstly, ion channels are commonly used to transfer ions from one aqueous medium to another, as in a biological cell, where the channel is embedded in a hydrophobic lipid membrane; hence for consistency with the aqueous medium the channel ends must be polar. The second function of the carboxylic groups is to generate the fluctuating electric field referred to in the Introduction with a view to gaining insight into its possible role in influencing ion migration.

Computations
The system energy is calculated by summing contributions for all the constituent atoms that are derived from specified force fields. These are: (i) intramolecular two-centre contributions such as those describing the stretching of a bond of length x AB , the threecentre deformation of bond angle θ ABC and an n-fold 4-centre torsion φ in the group AB-CD around bond B-C: The force field also contains (ii) intermolecular ('non-bonded') energy terms between atoms AB with separation r which are often written as a combination of a repulsive and an attractive term and often referred to as a 'van der Waals' energy contribution, and (iii) coulombic terms between the atoms A and B which bear electrostatic partial charges q A and q B: The force constants k AB , K ABC and torsional barriers V 0 and the potential parameters P and Q for the C, H and N atoms are taken from the Universal Force Field (UFF) [5]. The water solvent was described by Berendsen's SPC model which has been shown to simulate successfully the solvation of the ions by the water molecules so that no additional electronic polarization effects were allowed for [6]. Åquist potentials were used for Na + , and those of Lee and Rasaiah for Cl − [7,8]. Potential functions for heteroatomic pairs were calculated by the geometric-mean rule.
The interatomic Coulomb contributions to the lattice electrostatic energy were calculated from the partial charges on the atoms of the system and their separations. The Gaussian package [9] was employed to calculate atomic charges which the MD code (DL_POLY version 2.2) used for the coulombic eneZZrgy components; these quantities were summed according to the Ewald method [10,11] which is incorporated in DL_POLY. The definition used to calculate the sets of atomic charges {q A } was the natural population analysis (NPA) as their prediction of acid pK a values was shown to be superior to those from other definitions [12]. The molecular geometries (bond angles and lengths) required for the MD calculation were obtained from diffraction studies in the literature (where these were available; where they were not Gaussian was used to derive geometries from Hartree-Fock calculations on geometry-optimised partially 'hydrogenated' crown ether molecular components of the channe). In order to restore the original structures, due allowance was made for the removal of those H atoms which had been added to create the components upon which to perform the quantum calculations.
MD structural relaxation over a million timesteps (10 -9 s) might raise a doubt whether the original channel structue is sufficiently preserved for the intended investigation to be performed. However despite some dropping of the original translational symmetry of the channel that is visible in Figure 1(c) and in the potential profile in Figure 2 (see below) a comparison of the atomic structure of the CE molecule after such an MD run assures that although a small curvature develops along the z axis the final relaxed structure is nevertheless sufficiently similar to the intial one that the main features of the original condition are sufficiently preserved that the model may be used in the computations. The ionic migrations were followed after firstly allowing the complete MD system to relax around the initial positions of the Na + and Cl − ions, fixed at the axial coordinates z = 0.0 and z =25.0 Å shown in Figure 1(c).

Axial Electrostatic Potential
The electrostatic potential V(z) at points between the centres of successive rings along the z axis of the crown ether channel that is generated by the charges q i on the n atoms of the 12 CE rings was calculated from = ∑ ∑ (6) V(z) was averaged over the million timesteps of the MD run and the result is shown in Figure 2. The overall negative value of the potential in the resulting profile is due mainly to the partial charges of the oxygen atoms, which are closer to the axis than the other atoms. The effect of the positive axial potential, and its asymmetry might therefore be expected to distinguish between the Na + and Cl − ions and possibly their directions of migration in the channel.

Electric Fields and the Motions of Atoms in the Channel
As ions, the migrants should respond to the axial potential generated by the channel atoms but also be sensitive to an externally applied electric field.

Figures 3 and 4
show the response of the Na + and Cl − ions to electric fields ε applied initially along the positive z axial direction (though this was changed in the investigation of migrant directionality that will be described below). The units of the fields in the Figures are VÅ −1 and the response of the Na + and Cl − is shown when fields of the indicated strengths are applied along the positive axial direction (z). While these fields are larger than would be applied in a laboratory experiment, nevertheless in a simple ionic system a calculation would show that over short distances in the vicinity of a bare ion electric fields may indeed be of the order that is considered here. The behaviour of the two ionic species is clearly different, the Cl − ion showing significantly more mobility than Na + , whose progress is hindered by the attractive forces of the channel's negatively charged O atoms. The results of the calculations indicate that both the Na + and Cl − ions require a minimum axial field to of about 0.5 VÅ −1 to show any displacement but that when such a field is applied the Na + ion moves about 6 Å along the axis but atomic plots show that it then 'adheres' to a CE oxygen atom from which it is not dislodged by greater electric fields. However the same fields enable the Cl − ion to migrate almost freely through the negative potential of the channel (and beyond).

The Channel Environment
The role of the fluctuating charges on the channel was examined by controlling the motions of the CE atoms by clamping atoms in four conditions A, B, C and D, which are defined as follows: A No restraints were imposed on any channel atoms. B The six O atoms of each CE ring were fixed. C The atoms of the six COOH atoms on each terminal CE ring were clamped. D The CE oxygen ring atoms and the terminal-ring COOH were clamped As well as revealing the greater mobility of the Cl − ion over Na + the trajectories in Figures 3 and 4 show other distinguishing features. A consistent behaviour of the A to D conditions is that for small applied axial electric fields ε << 1.0 the migrations of both migrant ion species cannot occur without the fluctuations of local fields that are generated by the vibrations of the O atoms. At higher axial fields although (non-intuitively) of these, case D, where all the O atoms are stationary, results in faster Na + migration than cases B and C where less O atoms are clamped the Cl − ion behaves similarly: case A, in which none of the O atoms are clamped, results in a mobile ion, which passes through the channel with little inhibition (save a kink in the trace in Figure 4 when the Cl − ion reaches the end of the channel and enters into the vacuum). The other traces in Figure 4 reveal that the Cl − migration is inhibited by the O clamping, although at higher fields the inhibition does not always occur. Interestingly, at ε = 0.2 Figure 3(g) shows the trajectory of the Na + ion conducting an oscillatory motion which is reminiscent of ion behaviour in the 5-CE-15 channel [3].

Solution Environment
In order that the investigations more closely simulate commonly-studied molecular systems the channel was next considered to be an aqueous solute in which it might be expected that the coulombic atomic charges on the fluctuating water molecules might play a role in the dynamics of the ion migrants. Accordingly, water molecules were placed to surround the channel and occupied sites in its interior with a density approximating that outside.
In the following investigations when water is present the electric field ε field is greater than 0.75 and the migration of the ions is shown in Figure 5. Three cases were considered:  Figure 3 for Na + .
A There are no restrictions on the vibrations of the surrounding atoms. B The water atoms (only) are clamped. C Water atoms were clamped and all the charges on the channel O atoms were assigned as zero. Figure 5 shows that the general results of the investigation were broadly similar to those in the 'dry' environment. In case A where the dynamics of the atomic environment was unrestricted, ε fields of minimum strengths ca. 0.7 (as in the 'dry' case) were needed to initiate the ion migration.
As the effect was commented on in Section 4.2 we note that in condition B, where the water molecules are clamped, here too there is a range of fields, roughly from 0.9 to 1.1, causing oscillatory ion motion (Figure 6 below). In Section 4.2 above this had been noted to have occurred at higher fields in restricted dynamics C and D in the dry system (see Figure 3).

Temperature
The horizontal traces in Figures 2 to 5 denote halting sites of the migrants along the channel axis, albeit where the ions may undergo oscillations. The sites suggest the presence of energy wells in which the ions are (perhaps temporarily) entrapped during their migration. If this is correct it would imply that the ions require the acquisition of thermal energy in order to be liberated from the wells and continue their motion along the channel. The results of performing MD simulations at different temperatures respectively are given in Figure 7 (a) and (b), which show the channel trajectories of the Na + and Cl − migrants respectively in water at the temperatures indicated on the traces. Because of the randomness underlying the dynamics the results were too disparate to make an Arrhenius plot involving the rate and temperature to determine an activation energy for the motion. Nevertheless the conclusion is evident that the high-temperature traces are to the left of the Figures, indicating that at low temperatures (on the right) the ions require longer times to pursue or continue their motions along the channel. The inherent randomness that affects the dynamics of both ion species (and particularly because Cl − travels almost unimpeded through the channel, showing little additional features of interest) decides for us that the subsequent description of ion migration in this work will dwell more on the behaviour of the Na + than the Cl − ion. In addition, all investigations below will be in the water-inserted channel.    Figure  respectively show the axial displacement and the force on the ion over this interval. The force plot shows that after reaching the point A 1 the ion is subjected to a negative force, causing it to rebound in the reverse direction to B 1 . Continuing its forward progress to A 2 , Na + then responds to a second slight negative force that brings it back to B 2 .

Energy Changes
Intermolecular energy profiles over the same interval were calculated using eqn. (5) for the electrostatic component and eqn. (4) for the van der Waals or 'non-bonding' energy and the results are displayed in the two parts of Figure 9. The results indicate in accessing its first halting site the Na + ion encounters an activation energy of about 7 ev. A comparison of the energy scales of (a) and (b) shows that the dominating contribution is from the coulombic energy in part (a).
Recall that the process of clamping atoms in Sections 4.3 and 4.4 showed that the motions of charged atoms in the environments of the ion migrants influenced migration by providing fluctuating fields. These impose forces that boost dr liberation of the migrants from energy wells created by the atomic configurations that have entrapped the ions at sites along the channel. A comparison of the magnitudes of the energy scales of parts (a) and (b) of Figure 9 shows that the energies imparted to the ions arise overwhelmingly in the fields produced by the rapid changes in the positions of the atomic charges rather than of those accruing from the motions of the atoms that generate non-bonded forces.

Directionality of Channel Migration
The translational asymmetry of the channel and of its axial electrostatic potential were referred to in Section 4.1, and there it was suggested that features of the ion migrations might consequently depend on their direction of their motion in the channel. The effect of the channel's asymmetry on the ion dynamics was explored as follows. At various strengths of axial electric field ε the ion migrants conducted their trajectories in the same MD system as before, which consisted of the aqueous-immersed channel.
For this exploration two cases were considered. In the case 1, as described in Section 4, the Na + ion began its migration in a relaxed lattice in which its axial coordinate at low z and the Cl − ion from high z. In this case E was aligned along positive z so as to drive Na + in this direction and the Cl − ion towards the point z = 0.0. In case 2 the starting positions of the ions were interchanged and the direction of the electric field was reversed. Figure 10 displays the results of this investigation where four magnitudes of ε from 0.6 to 0.9 were considered, the blue and red colours of the traces distinguishing the field directions. The most notable observation in the Figure is that at the lower magnitudes of field (ε = 0.6 and 0.7) the Na + does not depart from its home site in case 1 until ε ≥ 0.8 but in case 2 a reversal of the initial sites of the ions and the of field direction has the effect of initiating the migration of Na + but preferentially enhancing the mobility of the Cl − ion in comparison to that of Na + .
The traces in Figure 10 are of two kinds. Some show little or no change to the applied electric field; others form sharp 'steps' indicating a more rapid response to the field. The Figure shows that the 'stepped' traces are more common for Cl − than for Na + at low ε. Thus the greater field-response of the former ion is consistent with the observation in Section 4.2 of less inhibited trajectories shown by the Cl − ion. It suggests a property of the aqueous channel-field system to act as a partial 'rectifier' of the (Na + , Cl − ) ion current at low values of the axial electric field.
Because of the non-uniform 'step' responses to the field it is not possible to ascribe dynamic variables to the ion dynamics as could be done if they were in a homogeneous medium. An example will illustrate this point. The slow progress shown by the blue trace for the Cl − ion as it migrates in the case 1 condition at ε = 0.7 leads to a diffusion coefficient of 8.1× 10~9 m 2 s~1. But when the condition is changed to that of case 2 the ion's diffusivity (which now responds sharply to the new field direction and channel conditions) jumps to 53.1× 10~9 m 2 s~1. This seems to obviate the use of diffusion parameters to characterise the directionality of the migrations

Trajectory Oscillations
We should like to find the cause of the oscillations in the Na + ion's axial migration that appear in the migration trace in Figure 8(a) and in the potential energy profile in Figure 9(a). In Figure 11(a) the axial coordinate z(Na + ) is reproduced from Figure 8 while Figure 11(b) shows how its radial coordinate r(Na + ) (i.e. its distance from the z axis) varies over the same time interval. In parts (a) and (b) of the figure the numbers in red indicate the serial number of the CE ring which the migrant ion has reached at the timestep indicated. When the oscillations start at ca. 1000 timesteps Na + succumbs to the attractive forces of the negatively charged O atoms of the CE rings which cause it to veer off its original axial path (which until this point had been less than 2 Å from the z axis). In Section 2 the channel diameter was stated to be 5.9 to 7.7 Å, so in the subsequent part of the ion trajectory with radial coordinate r < 3.0 Å describes how the ion departs for the channel walls by a presumably complex motion after which molecular plots illustrate an association with the CE's oxygen atoms. The radial coordinate of the Cl − ion in Figure 11(c), which is plotted with the same scale as (a) and (b) of the Figure, shows that the anion, in contrast, travels on a near-axial path for 4000 timesteps which was discussed above by reference to Figure 4(b) as an almost unimpeded trajectory over the 31 Å length of the channel.

The Dynamics of an Ion migrant Responding to Charge Fluctuations in the Channel
The effect of the motions of the COOH groups on the terminal CE6 rings was investigated by considering the changing electrostatic interactions of the atoms of each COOH group as it undergoes rotations about the C-C bond linking the carboxylic C atom to the carbon of the CE6 ring (see Figure 1).We shall denote this 'pivotal' carbon C*. The carboxylic rotations on the end rings were facilitated in the simulation by initially rotating each COOH group into a configuration with its C-C* bond along the channel's z axis. This COOH configuration is referred to as perpendicular since this torsional bond is 90˚ to the CE6 ring, while that in which the carboxylic group makes normal (tetrahedral) angles with the ring is termed the natural configuration. Imposing a torsion on the carboxylic group is then the computationally simple exercise of rotating it around the z axis. After conducting the torsion in the perpendicular conformation the group is restored to its natural configuration.
The coordinates of the COOH group which is to be rotated by an angle φ about the torsional C-C bond are (x*, y*, z*). Those of COOH in the 'perpendicular' configuration in which the C*--COOH bond is parallel to z are (x', y', z'). The resulting coordinates (x, y, z) of each of the four COOH atoms in the natural configuration are then given by x = (x' -x*)cos φ + (y' -y*)sin φ + x* y = [(x' -x*)sin φ + (y' -y*)cos φ]cos θ -(z' -z*)sin θ + y* (7) z = [(x' -x*)sin φ + (y' -y*)cos φ]sin θ -(z' -z*)cos θ + z* where θ = ½(τπ) is the angle through which C*-COOH is rotated to produce the perpendicular configuration, and τ is the tetrahedral bond angle. With the Na+ migrant ion in Figure 1 placed on the channel axis at the origin of the coordinate system, the equations (6) are used to calculate the distance r = (x 2 + y 2 + z 2 ) ½ and direction of each of the four carboxylic atoms from the migrant and hence the coulombic forces that they exert on the migrant ion.
In order to calculate the effect of a torsion φ on the energetics of the system we require the derivatives of the coordinates in eqn. (6) with respect to φ which can then be used in the expression of the φ-dependent force f i exerted by atom i of COOH. Here atom i (i = 1 to 4) with partial charge q i is at a distance r i from the migrant ion at the origin of the coordinates. The force fi from the four atom i atoms of the six COOH groups is where Q is the charge on the migrant ion. (There would of course be another, much smaller, contribution from the carboxylic groups on the terminal ring at the opposite end of the channel, which is not considered). When the φ-dependent distances r i from the carboxylic atoms to the Na + migrant are calculated using eqns. (7) the force f(φ) on the ion is provided by eqn. (8). This is plotted in Figure 12 for three revolutions around the C* ̶ COOH bond providing a series of impulse spikes which alternate between positive and negative values as the carboxylic groups rotate. The magnitudes of these impulses, which are ̶ 100 nanonewtons in one direction and +50 nN in the other, occur once in a 360˚ revolution of the carboxylic group. They are of the order of 10 nN, which are in turn consistent with the values calculated from an analysis of the oscillatory behaviour of ions in a CE channel, which were also found to have magnitude 10 nN [3]. Let us consider how consistent this ~10 nN force found from Figure 12 is with the ion dynamics in the channel system that is encountered in the current investigation. In Figure 8(a) the Na+ migrant, after an induction time that keeps it in the site where it was placed, responds to an impetus that abruptly confers on it a velocity that, after 250 fs, displaces it to a point 20 Å along the channel axis, indicating a mean speed of 0.40 Å fs -1. In order to get to the midpoint of its trajectory, the Na+ ion with mass 3.82×10-26 kg, has acquired an acceleration of 2.2×10-4 Å fs -2 and hence a force of 0.76 nN. Naturally we cannot claim that this force, which is of the order of 1 nN, is consistent with the profile in Figure 12, which shows force impulses that are at least an order of magnitude greater. The discrepancy between the two dets of results may originate in the fact that as the six COOH groups in the model described in this section are necessarily rotating mutually in phase, requiring the introduction of the factor 6 in eqn. (7). As a result, the overall effect would be larger than that produced by the actual rotations or those generated by the molecular dynamics simulation. There is also the possibility that the impetus spikes of negative force, occurring briefly in a complete revolution of C* ̶ COOH produce indirect effects on the lattice which accumulate to result in the passage of the migrant ion but which are too fine to be resolved by the MD. But importantly, we have found that the torsional rotations of the carboxylic groups generate substantial positive and negative impulse fields that can influence the passage of a migrating ion, a fact which is consistent with the results described in Section 4 of this work.

Discussion
The investigation has shown that the migratory motion of Na + and Cl − ions through a molecular channel is sensitive to fluctuating electric fields that are generated by the thermal motions of the atoms of the solution environment. We have seen that ion motion in a channel is strongly dependent on coulombic forces but that even if the migration is 'passive' it may not be simple, being subject to structural, dynamic and coulombic forces associated with the channel. Thus the axial electrostatic potential is not the sole factor governing the ion migration. The relatively light resistance to the migration of the Cl − ion in comparison to that on Na + that was noted in Sections 4.2, 4.7 and 4.8 may perhaps be ascribed to the coulombic 'sticking' of the cation to the channel's electronegative atoms (here the O atoms of the CE rings) that is found elsewhere [3].
The importance of the channel structure is clear when it is considered that (a) when the applied axial electric fields are small compared to those produced by the disparate fields of the surrounding atom the atomic forces play a greater role (it was shown that on increasing the axial field beyond ε ~ 1.0 the migration rates do not always behave as expected by simply clamping a selection of charged atoms in the ions' environment) and (b) the ion dynamics is reflected in the channel's translational asymmetry that is shown in Figure  1(c) and electrostatic potential profile in Figure 2 as it is responsible for a difference in the dynamics between the forward and backward ion motion in the channel. The result (b) may be of importance in the design of nanomolecular systems that involve the transportation of ions.
The force calculations in Section 4.9 demonstrate that the torsional motion of the chemical groups on the end CE rings may indeed be instrumental in providing the fluctuating fields that are necessary to displace a migrant ion in the channel. This result may be allied with the biological fact that the presence of a phosphate group (provided by an energy-rich ATP) which becomes affixed to the end of an ion channel in a neuron cell, subsequently playing a role as a pump that delivers ions by passive transport through the channel; this occurrence may have a similar effect to that of the COOH groups, the effect of whose motions have been investigated in this work Because of the near-random nature of the fluctuating electric fields that are generated by the vibrations of all the atoms in the environment of the ion migrants it is impractical to attempt to quantify the fluctuating forces ab initio; a firstprinciples treatment of the problem would use an equivalent of the Langevin equation which expresses the net force on the ion of mass M tavelling with velocity v where the external force F ext is here Eq that is experienced by the ion of charge q in applied electric field ε [13]. The term F st is a rapidly varying stochastic force originating in the atomic fluctuations and ᴦ is a parameter associated with the viscosity or friction of the system. In the present work the 2 nd and 3 rd term in eq. (7) are those that are implied in the molecular dynamics. Without postulating a specific form for F st the MD uses the motions of the atoms to supply the pseudo-random force that in the form of eq. (7) that was used by Einstein to rationalise the Brownian motion of pollen grains in a liquid [4].

Conclusions
1) The channel progress of ions is sensitive to the motions of the channel atoms.
2) The progress of the ions is affected by the motions of the solvent atoms.
3) The interaction of the migrant ions with their surroundings is strongly coulombic. 4) Cl − ions migrate almost freely through the channel. 5) The Na + migrants encounter an activation energy of ~7 ev in the channel.