Non Local Dielectric Effects on the Dielectric Barrier of an Ion Channel by Willem Bruce Krayenhoff B.Sc., The University of British Columbia, 2004 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF Master of Science in The Faculty of Graduate Studies (Physics) The University Of British Columbia (Vancouver) March, 2009 c© Willem Bruce Krayenhoff 2009 Abstract In this thesis we study non-local dielectric models on the sub-nanometer scale. In particular, we focus on the effect of non-local electrostatics on the potential barrier of a water-filled ion channel through a lipid cell membrane, and on the interaction between ions within such an ion channel. A polar- ization energy functional is used with its parameters calibrated to roughly reproduce the wave-vector-dependent dielectric function ǫ(q) of water. The lipid membrane is still modelled as a local dielectric. This energy functional is discretized onto a lattice and minimized using local moves to find the energy and the electric and polarization field config- urations for a given charge distribution and dielectric profile. This method is used to successfully reproduce known theoretical results with and with- out an ion channel. Necessary, but not necessarily sufficient conditions for obtaining good results are also derived. The dielectric barrier of an ion channel is studied with these non-local dielectric properties of water, and it is found that the results are very sen- sitive to how the water-membrane interface is modeled. We conclude that more molecular dynamics simulations are needed to provide guidance as to how to implement the polarization energy density functional at this inter- face, as well as whether these energy density functionals, which match the properties of bulk water, need to be modified in order to apply inside narrow ion channels. However, in all instances that we explored the self-energy of the ion traversing the ion channel was substantially modified from its self-energy under the local model, suggesting that non-local dielectric effects may be one significant factor determining the dielectric barrier of an ion channel. ii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . viii 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Review of Local Electrostatics . . . . . . . . . . . . . . . . . 1 1.2 The Dielectric Barrier of an Ion Channel . . . . . . . . . . . 3 1.2.1 Literature Review of Classic Solutions for the Dielec- tric Barrier . . . . . . . . . . . . . . . . . . . . . . . . 5 1.2.2 Short Channels . . . . . . . . . . . . . . . . . . . . . 5 1.2.3 Long Channels . . . . . . . . . . . . . . . . . . . . . . 7 1.2.4 The Pair Potential of Two Ions in an Ion Channel . . 8 1.3 Nonlocal Dielectric Effects . . . . . . . . . . . . . . . . . . . 9 1.3.1 Linear Dielectrics . . . . . . . . . . . . . . . . . . . . 9 1.3.2 Non-Local Dielectrics . . . . . . . . . . . . . . . . . . 10 1.3.3 Calculating ǫ(q) for the Water Model, and Fitting it to Molecular Dynamics Results . . . . . . . . . . . . 12 2 Model and Numerical Methods . . . . . . . . . . . . . . . . . 14 2.1 Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.2 Minimization . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.3 Minimization by Local Moves . . . . . . . . . . . . . . . . . . 18 2.3.1 A Sufficient Minimization Algorithm . . . . . . . . . 19 2.3.2 Additional Minimization Techniques . . . . . . . . . . 19 2.3.3 Criteria for Determining when the Minimization is Complete . . . . . . . . . . . . . . . . . . . . . . . . . 24 iii Table of Contents 3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.1 Units and Default Values . . . . . . . . . . . . . . . . . . . . 27 3.2 Local Electrostatics . . . . . . . . . . . . . . . . . . . . . . . 27 3.2.1 Pair Potentials . . . . . . . . . . . . . . . . . . . . . . 27 3.2.2 Pair Potentials in an Ion Channel . . . . . . . . . . . 29 3.2.3 Parameter Constraints and Numerical Artifacts . . . 33 3.3 Non-Local Electrostatics in Bulk Water . . . . . . . . . . . . 35 3.4 Conditions for the Successful Simulation of a Water-Filled Ion Channel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.4.1 The Water-Lipid Transition . . . . . . . . . . . . . . 36 3.4.2 Testing the Scaling Relation (Eqn. 2.6) for Varying Water-Lipid Transition Lengths . . . . . . . . . . . . 39 3.4.3 Testing the Finite Size Effects . . . . . . . . . . . . . 40 3.4.4 The Infinite Channel . . . . . . . . . . . . . . . . . . 45 3.5 The Dielectric Barrier of an Ion Channel . . . . . . . . . . . 52 3.5.1 The Local Dielectric Approximation . . . . . . . . . . 52 3.5.2 The Lorentzian Dielectric Approximation . . . . . . . 56 3.5.3 Water Model 2 . . . . . . . . . . . . . . . . . . . . . . 57 4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 Appendices A Verifying Scaling Results and Determining Acceptable Mesh Sizes for the Water Model . . . . . . . . . . . . . . . . . . . . 69 iv List of Figures 1.1 Diagram of a water-filled ion channel with an ion at its center. 4 1.2 A simplified model of the field of a single ion in an ion channel. 6 1.3 A simplified model of the field when two ions of equal and opposite charge are of intermediate separation R < d < ξ inside an ion channel. . . . . . . . . . . . . . . . . . . . . . . 8 2.1 A diagram showing links and nodes from a portion of the grid in a z = constant plane. . . . . . . . . . . . . . . . . . . . . . 15 2.2 The X and Y coordinates throughout the minimization of U(x, y) = x2 + 10(x − y)2 when performed by individually minimizing U with respect to each of x and y in turn. . . . . 20 2.3 The effect of the over-relaxation constant on the time needed for minimization. . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.4 The decrease in energy during each sweep of the grid. The X- axis shows the time, normalized so that the entire simulation takes a time of one. . . . . . . . . . . . . . . . . . . . . . . . . 24 2.5 The first graph shows the decrease in energy during each set of twenty consecutive sweeps... . . . . . . . . . . . . . . . . . 26 3.1 A diagram showing the grid and ion configuration we are using. 28 3.2 A plot of ǫU VS d for several values of ǫ and grid size L. The mesh spacing a = 1Å, q = ±1 and we are using the local polarization model (κ2 = α = 0). . . . . . . . . . . . . . . . . 30 3.3 A plot of the numerical discrepancy from theoretical results, (ǫU + 14pid − C), VS ion separation d. . . . . . . . . . . . . . . 30 3.4 A plot of Ex/Ex,theoretical VS ion separation d for several grid sizes L. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.5 A diagram of the channel configuration showing how the ions can see each other in both directions due to the circular boundary conditions. . . . . . . . . . . . . . . . . . . . . . . . 32 3.6 Theoretical (eq. 3.8) and calculated values of ǫwEx VS ion separation d1 for various values of ǫw . . . . . . . . . . . . . . 34 v List of Figures 3.7 Energy VS ion separation for the local model (eq. 1.19), the two Lorentzian models (eq. 1.25), and the two water models (eq. 1.27). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.8 A diagram showing the dielectric profile of the grid, as well as the position of the ions inside the water-filled ion channel. 38 3.9 Energy VS ion separation for various mesh sizes a for the second water model. . . . . . . . . . . . . . . . . . . . . . . . 41 3.10 Same as Fig. 3.9, but with l = 0.75Å in (c), and l = 1Å in (d). 42 3.11 Ex due to the positive ion at the location of the negative ion VS ion separation d for various grid sizes L, using the second water model and the local model. . . . . . . . . . . . . . . . . 43 3.12 Corrected electric field Ex+Ec VS ion separation for various grid sizes L using the second water model and the local model. 44 3.13 Corrected pair potential U +Uc VS ion separation for various channel radii R using the first water model and the local model. 46 3.14 Corrected electric field E+Ec VS ion separation d for various channel radii R, using the first water model and the local model. 47 3.15 Ex,withchan − Ex,bulkwater VS the ion separation d for various channel radii R using the first water model and the local model. 48 3.16 Ewithchan − Ebulkwater VS the ion separation d for various channel radii R, using the second water model and the local model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.17 Corrected Pair Potential U+Uc VS ion separation for various channel radii R for the second water model and the local model. 50 3.18 Ex,withchan − Ex,bulkwater VS the ion separation for various channel radii R using the second Lorentzian model and the local model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.19 Diagram of the dielectric profile of our grid used to model an ion channel of finite length, and of the trajectory of the ion we are moving. . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.20 U(X)−Ubulkwater VS the X-position of the single ion present (see diagram 3.19). The local, linear model for polarization is used and the channel radius R is varied. . . . . . . . . . . . 54 3.21 Ubarrier = Umax − Uout VS channel radius R, using the local dielectric model. . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.22 Ubarrier = Umax − Ubulkwater VS the X-position of the single ion present (see diagram 3.19). The Lorentzian models for polarization are used, and the channel radius R is varied. . . 56 vi List of Figures 3.23 U(X)−Ubulkwater VS the x-position of the single ion present (see diagram 3.19). The second water model is used, and the channel radius R is varied. . . . . . . . . . . . . . . . . . . . . 57 3.24 A diagram illustrating the trajectory of the ion we are moving, now with a Y offset. . . . . . . . . . . . . . . . . . . . . . . . 58 3.25 Uwithchan−Ubulkwater VS the x-position of the single ion present (see diagram 3.24). The second water model is used, and the y-offset is varied. . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.26 A diagram with an ion near the edge of the ion channel that illustrates how the polarization may be able to take advantage of what is energetically cheap in each medium to lower the overall energy. . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 3.27 U(x) − U(0) VS the x-position of the single ion present (see diagram 3.24). The second water model is used, amembrane = αwater, and the y-offset is varied . . . . . . . . . . . . . . . . 62 3.28 Uwithchan−Ubulkwater VS the x-position of the single ion present (see diagram 3.19) for various channel radii. The local model is used above and the second water model withαmembrane = αwater (not 0 in the nonlocal case) is used below. . . . . . . 63 3.29 Uwithchan−Ubulkwater VS the x-position of the single ion present (see Fig. 3.19) for various channel radii and usingκ1,membrane = 10 so that ǫm = 1.1, not the default value of ǫm = 2. The local model is used above and second water model is used below. 65 A.1 Energy VS ion separation for various mesh spacings a, using the water model (eq. 1.27) with A = 3.5 and L = 12Å. In the first graph q0 = 2.5, whereas in the second graph q0 = 4. . 70 A.2 Energy VS ion separation for various mesh spacings a. Uses the water model (eq. 1.27) with q0 = 2.5, A = 3.5, and L = 30Å. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 vii Acknowledgements I would like to thank my supervisor, Joerg Rottler, for helping give me direction when I had no idea how to proceed, and for assistance in making this thesis as clear as possible. viii Chapter 1 Introduction Most of physics involves approximations, and dielectrics are no different. The typical linear dielectric is in fact an approximation that is valid for many materials, but only over long length and time scales, and for weak electric fields. In local dielectric models, the polarization at a point is assumed to de- pend only on the electric field at that exact point, and not on the polar- ization or electric field of surrounding points, and the relationship between these two fields is mediated by electric susceptibility χe(r) such that P(r) = χe(r)E(r). However, when electric fields vary rapidly on the nanoscale, this model is not a valid description of highly polar fluids such as water, in which the orientation of each water molecule is highly dependent on the orientation of nearby water molecules. In more realistic models, the polarization at a point also depends on the polarization at nearby points, and is thus ’non-local’. These non-local effects are only significant when electric fields vary rapidly over short length-scales on the order of Angstroms (Å). In the case of ion channels, which are responsible for ion-transport in and out of cells, the relevant lengths are measureable in Angstroms. Thus the non-local dielectric properties of water may play an important role in deter- mining the dielectric barrier of an ion channel, provided they are significant over length scales greater than the channel radius. Thus far, research into the dielectric barriers of ion channels has largely ignored these non-local polarization effects, so this thesis is primarily an exploration of whether the non-locality of the dielectric response of water is indeed an important effect that must be taken into account when considering ion-transport through ion channels. 1.1 Review of Local Electrostatics When subject to an electric field, matter tends to polarize. When the matter is non-conducting, the field is not too strong, and the field also doesn’t 1 Chapter 1. Introduction change too quickly in time or space, the dipole moment per unit volume P is roughly proportional to the applied electric field E at each point in space [3]: P(r) = χe(r)E(r), (1.1) where the electric susceptibility χe is the constant of proportionality (Note that we are assuming ǫ0 = 1). This polarization may be due to the orientation of already polarized molecules along the electric field lines, or due to the induced polarization of atoms or molecules that didn’t have a pre-existing dipole moment. In any event, changing polarization creates accumulations of bound charge, which in turn modify the electric field, which in turn modifies the polarization, ad infinitum. This circularity is solved by the introduction of the electric displacement field D, defined by D = E+P. By Gauss’ law we have ∇ · E = ρ = ρf + ρb, (1.2) where ρf and ρb are the free and bound charges respectively. Now ρb = −∇ ·P, (1.3) and so we get ∇ ·D = ρf , (1.4) which is Gauss’s law in terms of D. Note that, while ∇×E = 0 in electro- statics, D is not subject to this constraint because P is not subject to this constraint. Assuming we are considering a linear dielectric, D, like P, is proportional to E D(r) = ǫ(r)E(r), (1.5) where the constant of proportionality ǫ(r) = χe(r) + 1 (1.6) is the permittivity of the linear dielectric. It is convenient to work in terms of D because it depends only on the free charge, which is generally what we control or know. Thus it is convenient to also express P in terms of D as 2 Chapter 1. Introduction P(r) = χ(r)D(r), (1.7) Calculating Energies Additional results that we shall be using from un- dergraduate electrostatics include that the vacuum energy due to a charge configuration can be calculated from its electric field using: UV = ∫ 1 2 E2dV, (1.8) and that the full energy due to a charge configuration, including the potential energy stored in linear dielectrics, can be calculated from E or D using: U = ∫ ǫ 2 E2dV = ∫ 1 2ǫ D2dV, (1.9) If we imagine that there are many dipoles in a dielectric composed of positive and negative charges which are held together by springs, the extra energy in eq. (1.9) over equation (1.8) is the potential energy stored in these springs. 1.2 The Dielectric Barrier of an Ion Channel Fig. 1.1 shows a toy model of the configuration of an ion channel in a biological cell. The smallest ion channels have length M ≈ 2.5nm and radius R ≈ 3Å. The cell membrane has a low dielectric constant of ǫm ≈ 2, whereas the water in the channel has a high dielectric constant of ǫw ≈ 80, for a dielectric contrast of 40. Due to the low dielectric constant of the cell membrane, the self energy of an ion in an ion channel is much higher than that of an ion in open water (i.e. either inside or outside of a cell), meaning that there is a potential barrier that must be overcome for an ion to traverse the ion channel, which has implications for ion transport in and out of cells. Thus far we are neglecting many things, such as the presence of salt in the water, of multiple ions, of other transport mechanisms, and of non-local polarization effects over short length scales. In this thesis our goal is to explore whether this last effect must be taken into account, as the length scales involved may be short enough that non-local polarization effects play an important role. 3 Chapter 1. Introduction Figure 1.1: Diagram of a water-filled ion channel with an ion at its center. M is the thickness of the cell membrane, and thus is the length of the channel. R is the radius of the channel, which is assumed to be cylindrical. ǫw and ǫm are the permittivities of water and the cell membrane, respectively. 4 Chapter 1. Introduction 1.2.1 Literature Review of Classic Solutions for the Dielectric Barrier The Intuition Because of the large dielectric contrast between the water inside an ion channel (ǫw ≈ 80), and the surrounding cell membrane (ǫm ≈ 2), over intermediate distances the electric field of an ion in an ion channel is largely confined within the ion channel. This results in a 1D (linear) potential since, by Gauss’ law, this confined electric field is constant, and thus the potential changes approximately linearly as one moves lengthwise (in the x direction) within the channel, and we obtain: V = k + Dconfined ǫw x, (1.10) where k is an arbitrary constant. For long enough channels the field eventually ’leaks out’ of the channel and into the cell well, and we revert to the 3D potential corresponding to ǫ = ǫm: V = k + q 4πǫmr , (1.11) where r is the distance from the charge q and k is an arbitrary constant. To determine the dielectric barrier of the ion channel, we must consider the self-energy of an ion in the channel. One way to do this is to integrate over the electric displacement field, using eq. (1.9). 1.2.2 Short Channels We first assume that the channel is short enough and the dielectric contrast large enough that the electric field is largely confined within the channel until it exits from one of the channel’s two ends. If the ion is in the center, then the field in both directions will be equal and thus we get D = q 2A (1.12) by Gauss’ law (eq. 1.4), where A is the cross-sectional area of the channel and q is the ion’s charge. If we shift the ion somewhat to the right, the electric flux will tend to shift somewhat more to the right direction as well, since the distance it must travel before exiting the channel is shorter in this direction, and so this will lower the overall energy. The same thing occurs if we shift the ion to the left. 5 Chapter 1. Introduction Figure 1.2: A simplified model of the field of a single ion in an ion channel. We assume that the electric field is entirely excluded from region D, and that the fields of regions A and C together have the same energy as the ion in open water, so that the potential barrier is determined by the extra energy of the field in region B. Thus the self energy of an ion is highest when it is at the center of the channel, and decreases as it moves towards either end The Dielectric Barrier Given a path that an ion takes through an ion channel, the potential barrier of this path is the maximum potential attained on that path. The potential barrier of the ion channel is then the minimum over all paths through the channel of the potential barrier of each path. Since the self-potential of an ion increases as the ion approaches the edges of a channel, the minimal path is clearly along the center of the channel, and as discussed earlier is maximal along this path when the ion is in the middle of the ion channel lengthwise as well. Thus determining the potential barrier of the ion channel amounts to determining the potential of an ion at the exact center of the ion channel and subtracting from this the potential of an ion in open water. 6 Chapter 1. Introduction A Simplified Model The exact barrier cannot be calculated analytically, however we shall calcu- late it for a simplified model where the electric field configuration is shown in Fig. 1.2. In this model the electric field is entirely excluded from region D, and the fields of regions A and C together have the same energy as the ion in open water. Thus the potential barrier is determined by the extra energy of the electric field in regions B, whose volume is the same as the ion channel’s. Using equation (1.12) we get: Ubarrier = V D2 2ǫw = Lq2 8Aǫw , (1.13) where V = AL is the volume of the channel. There are certainly more sophisticated analytical models [6], however this captures the basic physics and the limiting behaviour of this configuration. 1.2.3 Long Channels For very long ion channels the electric field does eventually ’leak out’ of the channel over a characteristic length scale of ξ, where ξ increases with the dielectric contrast and with the diameter of the channel. This case may be realistically modeled as an infinitely long ion channel, which, assuming a cylindrical channel, can then be tackled by using a Fourier Transform on the differential equation for the potential φ. It is then possible to see three regimes: On length scales much less than the channel radius R, the potential is approximately that of an ion in open water. On length scales between R and ξ the potential is approximately 1D with dielectric constant ǫw ≈ 80, and on length scales greater than ξ the potential is 3D in nature with dielectric constant ǫm ≈ 2. The exact value of ξ is given by the equation: ξ = R √ ǫw 2ǫm √ log(2ξ/R)− γ, (1.14) where R is the channel radius, and γ ≈ .577 is Euler’s constant [6]. Thus ξ is proportional to the channel radius R, and increases non-linearly with the dielectric contrast. For a water channel (ǫw ≈ 80) and a lipid membrane (ǫm ≈ 2), we find ξ ≈ 6.35R. (1.15) 7 Chapter 1. Introduction Figure 1.3: A simplified model of the field when two ions of equal and opposite charge are of intermediate separation R < d < ξ inside an ion channel. The Dielectric Barrier Having identified the behaviour of the potential in these three regimes, we can find the approximate electric field by taking a gradient. We can then find the self-energy of an ion in the center of an ion channel using equation (1.9), and subtract from this the self-energy of an ion in open water, thus finding the channel’s dielectric barrier. This is done in reference [6], where they found that Ubarrier = 1 2 e2 ǫmξ (log( 2ξ R )− γ) + .26 e2 2ǫmξ , (1.16) where γ ≈ 0.577 as before, and ξ is give by eq. (1.14). This is just an analytical approximation, but they found that it agrees fairly well with numerical results [6]. 1.2.4 The Pair Potential of Two Ions in an Ion Channel One may also want to consider the pair potential of two ions of equal and opposite charges inside an ion channel. For simplicity we shall assume that this channel length can be approximated as infinite. If the ions were widely separated so there was little interaction between them, the potential would simply be the sum of the two charges’ self-energies individually, so we shall focus on the case where the ion separation is of intermediate distance d such that Rchannel ≪ d ≪ ξ. Thus the channel is significant but d is also small enough that there is negligeable leakage of the electric field into the surrounding cell membrane (Fig. 1.3). Thus the field in between the ions is D = q A , and the potential is Upair = V D2 ǫw = dq2 Aǫw , (1.17) 8 Chapter 1. Introduction where V is the volume inside the channel that is between the two ions. Thus we expect the potential U to increase linearly with ion separation d while d is in this intermediate region. 1.3 Nonlocal Dielectric Effects The linear dielectric theory discussed earlier describes the average polariza- tion in a region based on the average electric field in this region. Over very short length scales such as the molecular length scale, such an averaging procedure is no longer valid, and short-ranged correlations in the polariza- tion field become significant. For example, in highly polar fluids such as water, neighboring molecules orient themselves such that the polarization is strongly correlated over short lengthscales. Thus when the electric field changes rapidly over short length scales, the linear dielectric model of matter is no longer adequate. Here we shall discuss generalizations of the linear model of dielectric response that are more accurate over short length scales. 1.3.1 Linear Dielectrics First we ’derive’ the linear dielectric properties of matter using a simple model, to which we will later add more terms to account for non-local di- electric properties. One basic model used to visualize linear dielectrics and gain some in- tuition for them is that there are many many tiny dipoles imbedded in the material with the positive and negative charges connected by a spring. Thus polarizing the material comes at an energy cost related to the ’spring constant’ of these tiny springs, giving an energy density of: UP = 1 2 ∫ V κ1P 2dV. (1.18) Adding in the energy due to the electric field (using ǫo = 1 for conve- nience) we get U = 1 2 ∫ V (D−P)2 + κ1P 2dV. (1.19) Minimizing with respect to P, we find that P = D/(1 + κ1) (1.20) 9 Chapter 1. Introduction for each point r. Since ǫ is defined by D = ǫE = ǫ(D−P), we find that ǫ = 1 + 1 κ1 (1.21) or κ1 = (ǫ− 1) −1 = χ−1e . (1.22) We fit κ1to the macroscopic properties of water (ǫ ≈ 80) using eq. (1.22), and call this the ’local model’. So far the polarization at a point depends only on the electric field and κ1at that exact point, and not on the polariza- tion or the electric field of surrounding points. Thus eq. (1.19) is entirely ’local’. 1.3.2 Non-Local Dielectrics To make this model more accurate when the electric field changes rapidly over short length scales, we can add ’nonlocal’ terms which involve deriva- tives of P. This causes the polarization at each point to depend on the po- larization in its neighborhood, as well as on the electric field at that point, and causes rapid changes of the polarization over short distances to come at an energy cost. Thus in Fourier space ǫ(q) becomes wavevector-dependent. In isotropic materials, only the longitudinal components (D||P||E||q) are relevant, and their relationship is described by D(q) = ǫ(q)E(q). (1.23) We shall be working with several non-local energy functionals. In Fourier space, ǫ(q → ∞) = 1 for all materials, because molecules are incapable of responding to such a rapidly varying field and see only the average field of 0. For water, ǫ(q → 0+) = ǫmacroscopic ≈ 80. A popular starting point for introducing non-local effects has been to assume exponentially decaying polarization correlations, which leads to a Lorentzian form of the dielectric permitivity: ǫ(q) = 1 + 1 κ+ ( q q0 )2 , (1.24) which interpolates between the q → 0+ and q →∞ limits over a character- istic lengthscale of 1 q0 . This dielectric permitivity can be reproduced by minimizing the follow- ing energy density functional suggested by reference [4]: 10 Chapter 1. Introduction U = 1 2 ∫ V [ (D−P)2 + κ1P 2 + κ2(∇ ·P) 2 ] dV, (1.25) where κ1 determines ǫmac according to eqn. 1.21, and κ2 = 1 q20 (1.26) is chosen to reproduce the length-scale 1 q0 . We shall call this the ‘Lorentzian model’. We shall be using q0 = 3Å −1 and q0 = 2Å −1, which we will call Lorentzian model 1 and Lorentzian model 2, respectively (so the second Lorentzian model with q0 = 2Å −1 will have much more significant non-local effects). We choose these values of q0 so that the length scale 2π/q0 is of the order of two to three Å, which is comparable to the diameter of a water molecule. However, for water, neutron diffraction data and computer simulations [2] have shown that as q increases ǫ(q) goes through two poles before con- verging to 1 as q → ∞, and that between these poles ǫ is negative due to ’overscreening’. That is, P > D, meaning E is in the opposite direction to what you would expect. To more accurately fit this behaviour, we use the following energy functional also suggested by ref. [4]: U = 1 2 ∫ V [(D−P)2 + κ1P 2 + κ2(∇ ·P) 2 + α[∇(∇ ·P)]2]dV. (1.27) Again, κ1 determines ǫmacroscopic, while κ2 < 0 and α > 0 are chosen to match the main features of water obtained through molecular dynamics computer simulations and experimental data. We shall call this the ‘Water Model’. Notice that the new non-local terms in the energy functional both involve derivatives, and as such are not significant over long length scales with slowly varying electric fields and thus slowly varying polarization, which leaves only the local polarization term we are familiar with. In fact, transverse correlations of the polarization field are also present, and can be modeled by the inclusion of an additional energy density term: 1 2κ3(∇ × P) 2. However this term only contributes to the energy in inho- mogeneous media or when there is an interface, and we as yet lack clear guidance from experiments and molecular dynamics simulations as to how to model interfaces, so we shall be neglecting this term in this thesis. 11 Chapter 1. Introduction 1.3.3 Calculating ǫ(q) for the Water Model, and Fitting it to Molecular Dynamics Results We move into Fourier space to calculate ǫ(q), which we will later match to values determined through molecular dynamics simulations. First we substitute D = D0e iq·r and P = P0e iq·r into eqn. (1.27), assuming that D0, P0 and q all point in the same direction (So that we are only considering longitudinal correlations, as discussed earlier). We drop the integral as well, so that we are looking at the energy density functional. This gives U(q) = 1 2 (D0 −P0) 2 + 1 2 κ1P0 2 + 1 2 κ2(q ·P0) 2 + 1 2 αq2(q ·P0) 2. (1.28) Taking the gradient with respect to P (to find the minimum), we get: U(q) = P0 −D0 + κ1P0 + κ2(q ·P0)q+ αq 2(q ·P0)q = 0. (1.29) Solving for P0 we get (assuming all vectors point in the same direction): P0 = 1 1 + κ+ κ2q2 + αq4 D0. (1.30) Thus χ(q) = 1 1 + κ+ κ2q2 + αq4 , (1.31) and ǫ(q) = 1 + 1 κ+ κ2q2 + αq4 . (1.32) The maximum susceptibility (equation 1.31) occurs at q0 = √ κ2/(2α), where it has a value of A ≡ χ(q0) = 1/(1 + κ + 1 2κ2q 2 0). We have values for the wavevector q0 and amplitude A = χ(q0) of peak susceptibility from molecular dynamics simulations, and so we choose κ2 and α to match them as follows: κ2 = 2 q20 ( 1 A − 1− κ) (1.33) α = − κ2 2q20 (1.34) 12 Chapter 1. Introduction Computer simulations by Bopp et al [2] suggest q0 = 3Å −1 and A = 2, which we shall refer to as ’Water Model 1’ or simply ’Water 1’. However, subsequent work has found that this result is highly sensitive to the local charge distribution of water molecules, and suggests instead q0 = 2Å −1 and A = 5 [1], which we shall refer to as ’Water Model 2’ or ’Water 2’. Given the variability in the literature, we shall study both sets of values, and in some cases their average. We shall also study the Lorentzian model discussed earlier with q0 = 3Å −1 (’Lorentzian 1’) and q0 = 2Å −1 (’Lorentzian 2’). 13 Chapter 2 Model and Numerical Methods 2.1 Discretization To find D and P we shall be minimizing the energy functional eq. (1.27) subject to the Gaussian constraint: ∫ S D · dA = Qenclosed. (2.1) Note that, at the minimum of the energy functional (eq. 1.27), ∇×(E) = ∇×(D−P) = 0, so we do not need to include this as an additional constraint for electrostatics. We use an on-lattice scheme with a mesh spacing of aÅ and with an auxiliary curl component to the electric field, which is reduced to zero by the minimization process. The ions reside on the nodes, whereas the dis- placement (D) and polarization (P) fields reside on the links connecting the nodes (see Fig. 2.1). We denote the node with cartesian coordinates (x,y,z) by Nx,y,z, where we use the mesh spacing a as the unit of length. Let ei denote the (i + 1)th unit vector: e0 = (1, 0, 0), e1 = (0, 1, 0), e2 = (0, 0, 1). The link that connects the Nx,y,z and N(x,y,z)+ei nodes is then denoted by Lix,y,z. We use periodic boundary conditions, so that if the dimensions of the grid are n× n × n (n nodes in each direction), the L0n,y,z link connects the Nn,y,z andN1,y,z nodes. For simplicity, in our description below the subscript n+1 shall be interpreted as 1, and the subscript 0 shall be interpreted as n, where appropriate, in accordance with these periodic boundary conditions. Each node Nx,y,z has a free charge, which may be zero, which we shall denote by qx,y,z. Each link has a displacement field, which we denote by D(Lix,y,z) or D i x,y,z, and a polarization field which we denote similarly. Since we shall be working with systems of inhomogeneous dielectric pro- file, such as water-filled ion channels through lipid membranes, each node 14 Chapter 2. Model and Numerical Methods Figure 2.1: A diagram showing links and nodes from a portion of the grid in a z = constant plane. The nodes are denoted by their coordinates in 3D cartesian space. The links have an additional superscript that denotes the direction the link is pointing in (x = 0, y = 1, z = 2). Thus, while it is not shown in the diagram, the L2x,y,z link connects the Nx,y,z and Nx,y,z+1 nodes. 15 Chapter 2. Model and Numerical Methods also has values for κ1, κ2, and α, which we shall denote by κ1,x,y,z = κ1(Nx,y,z), κ2,x,y,z = κ2(Nx,y,z), and αx,y,z = α(Nx,y,z). Since P2 = P 2x + P 2 y + P 2 z , in discretizing it is natural to express κ1P 2 as κ1,xP 2 x + κ1,yP 2 y + κ1,zP 2 z , which can then be easily discretized onto our grid. We then require a value of κ1 for each link, which we obtain from the node values using simple arithmetic averaging, and denote by κ1(L i x,y,z) or κi1,x,y,z. We discretize the divergence term in equation (1.27) around each node. We denote the flux of polarization out of node Nx,y,z by a 2divP (Nx,y,z) ≡ a2divPx,y,z, where divPx,y,z is: divPx,y,z = P 0 x,y,z + P 1 x,y,z + P 2 x,y,z − P 0 x−1,y,z − P 1 x,y−1,z − P 2 x,y,z−1. (2.2) Using this notation, Gauss’ law (eq. 2.1) translates into a linear equation relating the values of D on the six links connected to each node: a2divDx,y,z = qx,y,z (2.3) We also discretize the grad-div term (α[∇(∇ · P)]2) in equation 1.27 around the links. Let us denote the component of the grad-div term centered around Lix,y,z by gradDivP i x,y,z. Choosing i = 0, for example, we get: gradDivP 0x,y,z = divPx+1,y,z − divPx,y,z (2.4) Since the grad-div term is calculated around each link, we also interpolate α onto the links using arithmetic averaging. The Discretized Energy Functional Let Links denote the set of all Links, and Nodes the set of all nodes in our grid. Using the above definitions and discretizations, the energy functional (eq. 1.27) becomes: U = .5a3[ ∑ N∈Nodes κ1( divP (N) a )2, (2.5) + ∑ L∈Links (D(L)− P (L))2 + κ1P (L) 2 + α( gradDivP (L) a2 )2] subject to constraint eq. (2.3). All instances of the mesh spacing ‘a’ are explicit in equation 1.27 and in its constraint, so we can see that we would get the same answer if we 16 Chapter 2. Model and Numerical Methods substituted κ2,x,y,z → κ2,x,y,z a2 , αix,y,z → αix,y,z a4 , and qx,y,z → qx,y,z a2 , ran our minimization procedure with mesh spacing 1Å and the same number of nodes, then multiplied the resulting energy by a3. We can furthermore undo the last substitution, qx,y,z a2 → qx,y,z, for each charge and multiply the resulting energy by an addition factor of a−4. To sum up, if we express the energy as a function of the mesh spacing a, of the number N of nodes spanning the grid, and of the non-local dielectric constants on each node or link, U = U(a,N, κ2,x,y,z, α i x,y,z, qx,y,z), we have the relation: U(a,N, κ2,x,y,z , α i x,y,z, qx,y,z) = 1 a U(1, N, κ2,x,y,z a2 , αix,y,z a4 , qx,y,z). (2.6) This is, in fact, how we ‘change’ the mesh size. 2.2 Minimization In this section we discuss how the grid is initialized as well as the various minimization techniques we use to find the equilibrium displacement and polarization fields and the equilibrium energy efficiently. Initializing D Our strategy is to first initialize D on the lattice such that the Gaussian constraint eq. (2.3) is satisfied on all nodes for the given charge distribution, and then to make no changes during minimization that would change divDx,y,z on any node. Since the location of each charge is fixed, this ensures the Gaussian constraint is always zero. Note that satisfying this constraint (eq. 2.3) for all nodes requires that the net charge on the lattice be zero. When there are only two equal and opposite charges, this is accomplished by charting a path of links that connects them, and then initializing all of those links with field D = ±q, as needed to satisfy eq. (2.3). When there are more than two charges, as for example when we have a single charge q and a smeared out countercharge, we initialize D as follows: We first sweep the x-oriented links, setting D on each link such that excess divergence accumulates on all the nodes N1,y,z, and Gauss’ law is satisfied elsewhere. We next sweep the y-oriented links, then the z-oriented links until Gauss’ law is satisfied on all links (assuming charge neutrality). Our strategy can be described using pseudocode as follows: 17 Chapter 2. Model and Numerical Methods for i = 1 to n for all y and z, set D0i,y,z such that Gauss’ law is satisfied around node Ni,y,z end for # Comment - Gauss’ law is now satisfied on all nodes other than those # with an x-coordinate of one: N1,y,z. for i = 1 to n for all z, set D11,i,z such that Gauss’ law is satisfied around node N1,i,z end for # Comment - Gauss’ law is now satisfied on all nodes other than those # with x and y coordinates of one: N1,1,z. for i = 1 to n set D21,1,i such that Gauss’ law is satisfied around node N1,1,i end for # Comment - Gauss’ law is now satisfied everywhere, # provided the grid is charge neutral overall. 2.3 Minimization by Local Moves It is clear from the form of eq. (2.6) that the energy functional is a parabolic function of the polarization P on each link, and of the transverse degrees of freedom of the electric displacement field. Thus there are no local minima other than the global minimum. In theory, we should be able to calculate the first derivatives vector and the second derivatives matrix of this energy functional, and solve analytically for the fields and energy at the minimum. Unfortunately, this would be very difficult to implement, and is likely very innefficient computationally. We instead use a simple relaxation method with respect to local mini- mization moves on the D and P fields that satisfy the Gaussian constraint 2.3. A Local Energy Functional It is clear that changing the polarization or electric field on a single link only affects a small number of terms in equation (2.6) corresponding to links or nodes connected to the link in question. This formulation is thus entirely local, meaning we can calculate the change in 18 Chapter 2. Model and Numerical Methods energy due to a change in P or D on a single link with a few calculations instead of having to re-calculate the entire energy functional eq.(2.6). 2.3.1 A Sufficient Minimization Algorithm A fixed repertoire of three types of local moves suffices to find the minimum energy and the equilibrium polarization and electric fields. The polarization is unconstrained, and so the only move needed for the polarization is P → P + δP on each single link. The most obvious choice for δP is that which minimize U with respect to this one variable. However, as we shall discuss later, other choices are often preferable. Since the electric field is constrained, we require local and global moves on the transverse (curl) degrees of freedom of D which ensure that Gauss’ Law is still satisfied. We define a local curl move on each plaquette consisting of four links that form a loop. Thus, an example of a local curl move would be D0x,y,z → D 0 x,y,z + δD, D 1 x+1,y,z → D 1 x+1,y,z + δD, D 0 x,y+1,z → D 0 x,y+1,z − δD, and D1x,y,z → D 1 x,y,z − δD (The reader may find it useful to consult Fig. 2.1). Note that this move will ensure that ∇ × E = ∇ × (D−P) = 0 at the minimum of energy functional eq. (2.6), which is the expected result for electrostatics. A global move which wraps all the way around the grid with its periodic boundary conditions is also needed. An example is D0x,y,z → D 0 x,y,z+ δD for all x and for fixed values of y and z. Again, for these two moves on D, δD can be chosen to minimize U with respect to this variable. Any allowable change in the D and P fields can be produced through a combination of these three moves. In addition, since the potential is parabolic, by simply choosing δD to decrease U as much as possible for each move and repeating all the moves over and over, we will eventually approach the minimum. 2.3.2 Additional Minimization Techniques One finds that the minimization technique described above, while sufficient, is fairly slow. Here we give a two-dimensional example of the underlying reason. Suppose we are minimizing the energy functional U(x, y) = x2 + 10(x− y)2 (2.7) by individually minimizing U with respect to each of x and y in turn. What happens is that during the minimization process the x-y coordinates hug the line x− y = 0 because the 10(x− y)2 term is so much greater in magnitude 19 Chapter 2. Model and Numerical Methods 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 Y X Figure 2.2: The X and Y coordinates throughout the minimization of U(x, y) = x2 + 10(x − y)2 when performed by individually minimizing U with respect to each of x and y in turn. 20 Chapter 2. Model and Numerical Methods than the x2 term (Fig. 2.2). This is very innefficient and uses many many steps to approach the equilibrium configuration (x = y = 0). There are a couple of possible approaches to ameliorating this problem. In minimizing equation (2.7) it would clearly be very useful to have a move that changed both x and y at once, holding x − y constant. That is, we could couple moves so that they leave the most energetically costly terms unchanged. These coupled moves would then be in a new diagonal direction. Clearly, the above problem occurs because we don’t have moves that go in the directions needed to make minimization efficient in the space of possible minimization moves. We shall call this the direction problem. Coupled Moves The first approach to making our code more efficient is to somehow couple moves to obtain moves in new directions. Going one step further, we can instead solve the direction problem by performing combined moves that minimize U with respect to two or more moves/variables at once. We use two such combined moves: Firstly, on each plaquette we do simultaneous curl moves of δD onD and δP on P. We can then choose δD and δP to go straight to the minimum of U(δD, δP ). These moves affect only the (D − P )2 and κ1P 2 terms in the energy functional eq. (2.6), and this avoids the problem just described which would otherwise be caused by the mismatch κ1 ≪ 1 for water. Secondly, we couple the global move of magnitude δD onD with a global move of magnitude δP on P, and minimize U with respect to δD and δP simultaneously. This is again helpful because of the mismatch κ1 ≪ 1 for water. We also tried a third combined move. The nonlocal terms in the energy functional eq. (2.6) makes combined moves on the polarization of adjoining links very attractive as a method of avoiding the direction problem. If we did a combined move on the polarization of all the adjoining links in a row around the grid, the math should be tractable. For example, we could change P 0x,y,z → P 0 x,y,z + δPx for each x and fixed values of y and z. Let δP denote a vector of all the δPx for x = 1, 2, ...n and some fixed y and z. Then the gradient UδP = ∂U ∂(δP) depends only on the values of P and D on nearby links, and the second derivative matrix UδP×δP is a constant function of the dielectric profile on these links. More importantly, UδP×δP is a cyclical five-diagonal matrix, and so it is not difficult to solve for the location of the minimum of U with respect to δP: UδP×δP(δP) +UδP = 0. (2.8) 21 Chapter 2. Model and Numerical Methods We wrote our own cyclical five-diagonal matrix solver. Unfortunately, when applied it produced an occasional error whose source we were not able to identify. This could potentially have made our algorithm far more efficient when the nonlocal terms κ2 and α are large, as occurs for finer mesh sizes eq. (2.6). Local curl moves on D and P could be similarly combined, and would require us to solve cyclical tri-diagonal matrices. Unfortunately, we must leave this to future research. Over-Relaxation The second approach for solving the direction problem for the high-dimensional function we are working with is an over-relaxation technique. Simply put, when doing a minimization move on the polarization of a single link for example, choose δP to be m times the value that would minimize U(δP ), where 1 ≤ m < 2. It is necessary to sweep the links in an organized fashion in both directions, so that the changes can propagate quickly. As when applied to the Laplace equation, this technique aids faster equilibration over long length scales. We also apply this technique to the curl moves as well. When applying over-relaxation to the Laplace equation, it is possible to determine the ideal m theoretically from the dimensions of the grid. Unfor- tunately, the functional we are working with is much more complex. There are likely two ideal values for m, one for curl moves and one for polarization moves, and these ’over-relaxation constants’ likely depend on κ2 and α in addition to the dimensions of the lattice. We use the same value of m for both types of moves, and test the mini- mization time for bulk water with two ions of equal and opposite charge. We choose the directions in which we sweep the grid randomly for each sweep, which means there is some randomness to the speed of minimization. As Fig. 2.3 shows, the optimal value of the over-relaxation constant is in the neighborhood of m = 1.75 for these values of a, L, κ2, and α. Comparing with the time used form = 1, it appears that this over-relaxation reduces the time required by half. The range m = 1.4–1.9 appears to provide significant time reductions, so we will henceforth use values of m in the upper end of this range (m = 1.8–1.9), since we expect the optimal m to increase as L, κ2, and α increase. 22 Chapter 2. Model and Numerical Methods 300 400 500 600 700 800 900 1000 1100 1200 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 M in im iz at io n Ti m e Over-relaxation Constant Figure 2.3: The effect of the over-relaxation constant on the time needed for minimization. a = 1Å, L = 40Å, κ2 = −1, α = 2 23 Chapter 2. Model and Numerical Methods 1e-12 1e-10 1e-08 1e-06 0.0001 0.01 1 100 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 En er gy C ha ng e D ue to S w ee p Time Figure 2.4: The decrease in energy during each sweep of the grid. The X- axis shows the time, normalized so that the entire simulation takes a time of one. 2.3.3 Criteria for Determining when the Minimization is Complete During the minimization process we simultaneously sweep the grid with all the moves discussed above in both directions (except for global moves where sweeping in both directions provides no extra benefit), and we call this ’one sweep’. For each type of move we randomize which orientation of faces/links we sweep first, second, and third, and consequently there is a definite randomness to the amount by which the energy is reduced in each sweep, in addition to a clear long-term trend that appears to be roughly exponentially decreasing with respect to the number of sweeps (Fig. 2.4). While the energy decrease during one sweep has a significant random component to it, if we look at the energy change during groups of a number (i.e. 20) of sweeps at a time, the randomness tends to average out, as shown by Fig. 2.5 (left). The ratio between subsequent energy decreases also becomes fairly stable when we look at groups of twenty sweeps (Fig. 2.5), 24 Chapter 2. Model and Numerical Methods and so we use the roughly exponentially decreasing nature of the energy decreases with respect to the number of sweeps to estimate the amount that the energy will yet decrease by, and thus the error, using the formula 1 + r + r2 + r3 + ... = 11−r . We terminate the simulation when this error is estimated to be less than 10−8. Variations of the above graphs with the same essential features we are relying on are produced independent of the number of charges, their separa- tion, or any of the other parameters passed to the program, thus this appears to be a very reliable method of determining when we are sufficiently close to the minimum and can stop the minimization. 25 Chapter 2. Model and Numerical Methods 1e-08 1e-07 1e-06 1e-05 0.0001 0.001 0.01 0.1 2 3 4 5 6 7 8 9 10 11 δ E Number ’n’ of Sets of Sweeps 0.01 0.1 1 2 3 4 5 6 7 8 9 10 11 δE n /δ E n -1 Number ’n’ of Sets of Sweeps Figure 2.5: The first graph shows the decrease in energy during each set of twenty consecutive sweeps of the grid, while the second graph shows the ratio of the decrease in energy during each set of 20 sweeps to the decrease during the previous set of 20 sweeps. 26 Chapter 3 Results 3.1 Units and Default Values Unless otherwise specified, throughout this thesis all lengths are in angstroms (Å), all charges are in units of the charge of a proton (e), all fields are in units such that ǫ0 = 1 so that∫ E · dA = qenclosed, (3.1) and the energy U is measured in [e]2/[Å], so that the pair potential of two ions a distance R apart in vacuum is U = q1q2 4πR . (3.2) All charges q = ±1, unless otherwise specified. 3.2 Local Electrostatics We first check that we get expected results and relations in certain known situations and limits as a test of whether our code is working properly. 3.2.1 Pair Potentials Using the configuration depicted in Fig. 3.1, we first reproduce some ex- pected results for the pair potential of two ions using the purely local model of dielectric effects (κ2 = 0 and α = 0 in eq. (2.6)). The charge q = 1, so we expect that, when immersed in a medium with purely local dielectric properties, the energy should be U = −1 4πǫd + C ǫ , (3.3) where d is the ion separation and C is a constant. Thus ǫU = − 14pid + C should be independent of ǫ. As shown by Fig. 3.2, our computational results 27 Chapter 3. Results Figure 3.1: A diagram showing the grid and ion configuration we are using. L is the length of the width, height and length of the grid, a is the mesh spacing, and d is the distance separating the ions, all measured in Å. 28 Chapter 3. Results show that ǫU is indeed independent of ǫ to a very high degree of precision, and also agrees quite well with our theoretical prediction provided the ion separation d is at least several lattice spacings apart and d≪ L/2. Since ǫU = 14pid + C, in theory ǫU − 1 4pid − C = 0, so we graph ǫU − 1 4pid − C VS d to see what errors and corrections discretization and the circular boundary conditions have introduced (Fig. 3.3). The discrepancies from theoretical results are about an order of magnitude smaller than the variation in the original function. It appears that there are two primary origins: An error due to the discretization that is significant on short length- scales of under four lattice spacings, and a correction due to the finite volume and the circular boundary conditions, which appears to be significant for ion separations d > L/4 and varies roughly quadratically with d. We shall call these the discretization error and the finite grid-size correction, respectively. While the finite grid-size correction has only a limited absolute impact on the potential, we may sometimes be interested in the relative change in the electric field E due to the circular boundary conditions. According to Fig. 3.4, the relative correction in the x-component of the electric field Ex is less than 20% for ion separations d < L3 . Since the absolute correction is quite small L ≥ 3d or even L ≥ 2.5d should suffice for our purposes, and since the effects of the non-local terms in the energy functional eq. (2.6) decay far more quickly with distance than the local terms, we expect this condition to suffice in general. As for the short ranged discretization error, we have established that it exists and is significant, however we cannot make any prescriptions in general, as this error will also depend on the non-local parameters κ2 and α. Lacking clear guidance, we shall use d ≥ 3a. 3.2.2 Pair Potentials in an Ion Channel Here we expect that, given dielectric contrast sufficient to contain the field within the channel and ion separation significantly greater than the channel radius (d≫ R), the energy should depend linearly on d in accordance with eq. (1.17), and (ǫwU) should have a constant derivative with respect to ion separation d of: ∂(ǫwU) ∂d = 1 2 ǫ2wE 2A = 1 2 ( 1 A )2A = 1 2A , (3.4) where A is the cross-sectional area of the ion channel, ǫw is the permitivity of the water (inside the ion channel), and we assume q = ±1, as before. 29 Chapter 3. Results 0.16 0.17 0.18 0.19 0.2 0.21 0.22 0.23 0.24 0.25 0 2 4 6 8 10 12 14 16 εU d ε = 10, L = 30 ε = 80, L = 30 ε = 200, L = 30 ε = 10, L = 45 ε = 80, L = 45 ε = 200, L = 45 ε = 10, L = 60 ε = 80, L = 60 ε = 200, L = 60 Theoretical Value Figure 3.2: A plot of ǫU VS d for several values of ǫ and grid size L. The mesh spacing a = 1Å, q = ±1 and we are using the local polarization model (κ2 = α = 0). The results appear to be independent of ǫ. The expected theoretical value ǫU = − 14pid +C is also shown, with C chosen to match the numerical results at large ion separations. -0.007 -0.006 -0.005 -0.004 -0.003 -0.002 -0.001 0 0 2 4 6 8 10 12 14 16 εU + 1 /(4 pi d) - C d ε = 10, L = 30 ε = 80, L = 30 ε = 200, L = 30 ε = 10, L = 45 ε = 80, L = 45 ε = 200, L = 45 ε = 10, L = 60 ε = 80, L = 60 ε = 200, L = 60 Figure 3.3: A plot of the numerical discrepancy from theoretical results, (ǫU + 14pid − C), VS ion separation d. 30 Chapter 3. Results 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 0 2 4 6 8 10 12 14 16 E x /E x ,th d L = 30 L = 45 L = 60 theoretical Figure 3.4: A plot of Ex/Ex,theoretical VS ion separation d for several grid sizes L. Ex refers to the field due to the positive charge at the location of the negative charge in the direction of their separation (the x direction). The graph is independent of ǫ, thus the value of ǫ is not included. a = 1Å, κ2 = 0, and α = 0 (local polarization). Ex is calculated in two ways: Directly from the lattice, and via E ≈ ∆U∆d . For small values of d and relative to theoretical predictions, the E field values calculated directly from the lattice are larger than those calculated by finite difference, accounting for the sawtooth pattern observed. 31 Chapter 3. Results Figure 3.5: A diagram of the channel configuration showing how the ions can see each other in both directions due to the circular boundary conditions. The distance between the ions is then d1 or d2, depending on the direction, where d2 = L− d1. 32 Chapter 3. Results However, we are using circular boundary conditions, so the ions can see each other around the torus in both directions (Fig. 3.5), which should result in a significant finite grid-size correction. As Fig. 3.5 shows, the distance between the ions in one direction is d1 and the distance in the other direction is d2 = L− d1. There are corresponding displacement fields D1 and D2, which obey A(|D1|+ |D2|) = q (3.5) by Gauss’ Law. Integrating over the channel volume, we find that U = A 2ǫw (d1D 2 1 + d2D 2 2). (3.6) The field configuration and energy are found by minimizing this energy function with respect to D1 and D2 subject to constraint eq. (3.5), which gives U = 1 2ǫwAL (d1L− d 2 1) (3.7) and ǫwEx = ∂(ǫwU) ∂d1 = 1 2AL (L− 2d1), (3.8) which agrees with eq. (3.4) in the d1 ≪ L limit, as required. As Fig. 3.6 shows, for large values of ǫw and for d≫ R these theoretical calculations are in excellent agreement with the results from our model. 3.2.3 Parameter Constraints and Numerical Artifacts In this section and in Appendix A we have demonstrated that we get the expected behaviour in certain numerical limits, and have established the following rough numerical constraints on the parameters we use in our min- imizations that must be obeyed if we wish to get valid results. In this section we came up with the following constraints on d: 3a ≤ d ≤ 2 5 L (3.9) We do expect the electric field to be significantly altered due to local dielectric effects when L = 52d, however non-local effects decay very quickly, and so we do not expect finite size effects to substantially alter the non-local features of the potential, provided this is satisfied. 33 Chapter 3. Results 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0 5 10 15 20 25 30 ε w E x d1 εw = ∞ (theoretical) εw = 2 εw = 10 εw = 80 εw = 200 εw = 1000 Figure 3.6: Theoretical (eq. 3.8) and calculated values of ǫwEx VS ion separation d1 for various values of ǫw and with the ions situated in the center of the channel as shown in Fig. 3.5. We used a = 1Å, L = 75Å, ǫm = 2 for the cell membrane, and R ≈ 3.5, A = 37 ≈ π3.5 2. 34 Chapter 3. Results We have also established for local dielectric models that when we have an ion channel with a very large dielectric contrast betweem ǫw and ǫm, the potential and electric field are significantly modified by the circular boundary conditions, even when L is several multiples of d, as shown by the difference between eqs. (3.4) and (3.8). We explore whether non-local terms modify this effect later. In Appendix A, we find the additional constraint a ≤ 1 q0 (3.10) for the water models, for which the mesh spacing a must be small enough to resolve the oscillations in the potential. Keeping these constraints in mind, we now proceed to explore these non-local dielectric effects. 3.3 Non-Local Electrostatics in Bulk Water Fig. 3.7 shows the pair potential between two ions in water whose polariza- tion is modelled by each of the five energy functionals of section 1.3. For the Lorentzian models, the only effect of the non-local terms in the polarization energy functional appears to be an increased electric field at short ion separations. This is to be expected, since in the Lorentzian model ǫw(q) < ǫw(0) ≈ 80 ∀q > 0, and ǫw(q) → 1 as q → ∞, meaning that the effective dielectric permitivity is always less that 80 and really close to the ion where the field is changing very rapidly it approaches that of free space (ǫ = 1). Convergence to the local dielectric limit appears to occur within 3 and 4 Å for the two models, respectively. The water models are more interesting, exhibiting oscillations whose am- plitude decreases with increasing ion separation. Predictably, water model 1, which has the relatively smaller amplitude and higher frequency of maximum dielectric response, exhibits more frequent oscillations of smaller amplitude than model 2 does as the ions are drawn apart. For model 1, convergence to the local limit appears to have occurred by about 5Å separation, whereas for model 2 it appears that it will take significantly more than 8Å separation for convergence to occurr. Clearly non-local features dominate over short length scales. If one of the water models is indeed approximately correct and they still apply to wa- ter solutions of dissolved ions, this has implications for inter-ion distances: only certain inter-ion distances are energetically favourable, and ions in so- lution will tend to populate these minima. The presence of such potential 35 Chapter 3. Results oscillations in polar solvents has been pointed out many years ago [5]. Our results underline how sensitive the effects of the non-local terms are to the values of q0 and A used. 3.4 Conditions for the Successful Simulation of a Water-Filled Ion Channel First we would like to mention that, unless otherwise specified, the cell membrane shall be modelled as a purely local dielectric with ǫm = 2. 3.4.1 The Water-Lipid Transition Ion channels are water-filled passages, generally 3Å or more in radius, through lipid cell membranes. Non-local dielectric effects are due to short-ranged interactions between molecules, and we do not understand what kinds of interactions occur at the water-cell membrane interface. We are also un- aware of any molecular dynamics simulations that are capable of guiding us as to what implementation we should use for our implicite modeling of the polarization at this interface using an energy density functional. Let’s assume there are no such interactions between the lipid and water molecules. It is the ice-like structure of water over short lengthscales that causes the non-local polarization effects we are investigating, however at the interface the water molecules only have forces from water molecules on one side. This may substantially weaken all polarization effects (including local) of water, or may even have other effects (since forces from only one side may cause the edge water molecules to do something other than conform with the pattern established within the bulk water medium). Another effect that may be significant is that if you looked very closely at the water-cell membrane interface you might find that the boundary is wavy because it consists of individual water and lipid molecules, so to accurately model this it might be necessary to use a finer grid capable of resolving these ripples. Lacking clear guidance from molecular dynamics simulations, we shall ignore the above considerations and assume we are going to use an interpo- lation method for κ1, κ2, and α at the water-lipid interface over a transition region of length l which starts a distance of R − l2 from the center of the channel and ends a distance of R+ l2 from the center (Fig. 3.8). Recall that we first assign values to κ1, κ2, and α at each node, then use simple arithmetic averaging to determine κ1 and α on the links, as 36 Chapter 3. Results -0.03 -0.02 -0.01 0 0.01 0.02 0.03 0 1 2 3 4 5 6 7 8 U d Local Loretzian #1 Loretzian #2 Water #1 Water #2 -0.003 -0.002 -0.001 0 0.001 0.002 0.003 0 1 2 3 4 5 6 7 8 U d Local Loretzian #1 Loretzian #2 Water #1 Water #2 -0.0003 -0.0002 -0.0001 0 0.0001 0.0002 0.0003 0 1 2 3 4 5 6 7 8 U d Local Loretzian #1 Loretzian #2 Water #1 Water #2 Figure 3.7: Energy VS ion separation for the local model (eq. 1.19), the two Lorentzian models (eq. 1.25), and the two water models (eq. 1.27). The only difference between the graphs is the scale of the y-axis. a = 0.3Å, ǫw = 80, and L = 20Å. All energy graphs are shifted by a constant so that each graph’s right-most point has energy zero. 37 Chapter 3. Results Figure 3.8: A diagram showing the dielectric profile of the grid, as well as the position of the ions inside the water-filled ion channel. The dark regions represent the lipids in the cell membrane, the white region depicts the water- filled channel, and the gradiated transition regions show the region in which we interpolate between the dielectric properties of the lipid membrane and the water. The channel is cylindrical, the transition region has a thickness of l, and the radius R represents the distance from the center of the ion channel to the middle of the transition region. L and d are as before. 38 Chapter 3. Results described in section 2.1. However, this may introduce numerical artifacts in regions where these parameters are changing: If we look at the discrete energy functional eq. (2.6) we find that the value of κ1 on each link is only ’relevant’ to the polarization on that link, in that the terms in the sum that contain κ1 contain only the polarization on that same link. The value of κ2 on each node is relevant to the polarization on the six surrounding links, and the value of α on each link is relevant to the polarization on the eleven surrounding links. This assymetry between κ1, κ2 and α as well as the arithmetic averaging method used to obtain the link values of these parameters may introduce significant numerical artifacts if the mesh size a is too large compared to the transition length l. Thus we need a sufficiently long transition region, which raises the ques- tion of what to interpolate and how. One could for example use geometric or arithmetic interpolation of κ1. Geometric may be more appropriate, given the vastly different values of κ1 between the membrane and the water. An- other possibility is to interpolate χ (linearly?) and to determine κ1 from χ using eq. (1.31) with q = 0. We again lack any clear guidance from molecular dynamics simulations or experiment, and so we use linear interpolation of κ1, κ2, and α, just as we use arithmetic averaging of node values to determine κ1 and α on the links. However we raise all these possibilities to point out that we are forced to make many unjustified assumptions, whose validity future molecular dy- namics simulations will hopefully be able to shed light on. The following results are therefore merely exploratory. 3.4.2 Testing the Scaling Relation (Eqn. 2.6) for Varying Water-Lipid Transition Lengths We now expect the scaling relation (Eqn. 2.6) to hold so long as l≫ a. As Figs. 3.9 and 3.10 show, a significant transition length is necessary for this relation to hold. In the first graph l is so short that the actual transition length is determined by the mesh spacing a, and so as a decreases towards 0 the transition becomes sharper and sharper and the energies diverge (the minimization with a = 0.3Å didn’t finish). For the l = 0.5Å case it is not clear that convergence has occurred by the time a = 0.3Å, whereas in the l = 0.75Å case it appears that convergence has occurred for a ≤ 0.4Å. Since there is no apparent improvement in the convergence for the l = 1Å case, we conclude that the remaining discrepancy in the l = 0.75Å case is due to the fact that a = 0.45Å is close to the 1 q0 = 1/2Å maximum established for the mesh size (eq. 3.10). 39 Chapter 3. Results Given the radius of a water molecule (≈ 1Å), l = 0.75Å seems somewhat large, however, in order to simulate systems of any substantial size we shall need to use larger mesh spacings, and so we henceforth use l = 0.75Å. There is no visibly significant effect of the transition length l on the energies for the other four models (Local, Lorentzian models 1 and 2, and water model 1). The only surprise here is that there is no visibe effect for water model 1, however the non-local effects tend to be less pronounced and shorter ranged in the first model than in the second, which may explain this. With this choice of l = 0.75Å, we then have no new constraints on the mesh spacing a, since we have satisfied l≫ a by making l large, rather than a small. 3.4.3 Testing the Finite Size Effects In section 3.2.2 we showed that even for L ≫ d there were still significant finite grid-size corrections due to the circular boundary conditions and the field-containment within the ion channel. Thus we here explore whether this is also the case when significant non-local effects are present, using the second water model, whose non-local effects are most long-ranged. For all other models, the polarization is primarily local beyond 5Å ion separation. Again, we look at only the electric field near the maximum separation to best discerne the effects of the finite grid size. Fig. 3.11 shows that, the differences between the electric fields for different grid sizes L appear to be almost the same for the water and local models, so long as L > 2d. This suggests that these finite-size effects are almost entirely due to local dielectric effects for L > 2d, and so we plot Ex again ’correcting’ for the finite grid-size correction that we know is present for the local model. This corrective factor, which we shall call Ec, is obtained by subtracting equation 3.8 from equation 3.4 and dividing by epsilon. A corrective factor for the energy Uc can be derived similarly. They are given by: Ec = qd ALǫw (3.11) Uc = q2d2 2ALǫw (3.12) As Fig. 3.12 shows, all the local results collapse onto essentially the same line, as expected, and for L > 32 the water models do as well. Thus we conclude that we only need L slightly larger than 2d for the finite grid- size correction to be primarily due to local dielectric effects, which we can 40 Chapter 3. Results -0.002 0 0.002 0.004 0.006 0.008 0.01 3.5 4 4.5 5 5.5 6 6.5 7 7.5 8 U d l = .01 a) Local a = .33 a = .35 a = .4 a = .45 a = .5 -0.002 0 0.002 0.004 0.006 0.008 0.01 3.5 4 4.5 5 5.5 6 6.5 7 7.5 8 U d l = .5 b) Local a = .3 a = .33 a = .35 a = .4 a = .45 a = .5 Figure 3.9: Energy VS ion separation for various mesh sizes a for the second water model. The local results are included for comparison. ǫw = 80, R = 3Å and L = 16Å. The transition length l = 0.01Å in (a) and l = 0.5Å in (b). All energy graphs are shifted by a constant so that each graph’s right-most point has energy zero. 41 Chapter 3. Results -0.002 0 0.002 0.004 0.006 0.008 0.01 3.5 4 4.5 5 5.5 6 6.5 7 7.5 8 U d l = .75 c) Local a = .3 a = .33 a = .35 a = .4 a = .45 a = .5 -0.002 0 0.002 0.004 0.006 0.008 0.01 3.5 4 4.5 5 5.5 6 6.5 7 7.5 8 U d l = 1 d) Local a = .3 a = .33 a = .35 a = .4 a = .45 a = .5 Figure 3.10: Same as Fig. 3.9, but with l = 0.75Å in (c), and l = 1Å in (d). 42 Chapter 3. Results -0.0001 -5e-05 0 5e-05 0.0001 12 13 14 15 16 17 E x d L=30, Local L=32, Local L=36, Local L=40, Local L=30,H2O #2 L=32,H2O #2 L=36,H2O #2 L=40,H2O #2 Figure 3.11: Ex due to the positive ion at the location of the negative ion VS ion separation d for various grid sizes L, using the second water model and the local model. ǫw = 80, l = 0.75Å, a = 0.45Å and R = 3Å. 43 Chapter 3. Results 0 5e-05 0.0001 0.00015 0.0002 0.00025 0.0003 12 12.5 13 13.5 14 14.5 15 15.5 16 E x + E x ,c d L=30, Local L=32, Local L=36, Local L=40, Local L=30,H2O #2 L=32,H2O #2 L=36,H2O #2 L=40,H2O #2 Figure 3.12: Corrected electric field Ex + Ec VS ion separation for various grid sizes L using the second water model and the local model. ǫw = 80, l = 0.75Å, a = 0.45Å and R = 3Å. 44 Chapter 3. Results correct for by adding the terms given by eqs. (3.11) and (3.12). A possible explanation for this is that the amplitude of the oscillations (due to non-local dielectric effects) decreases rapidly with each oscillation, and so the oscillation in the longer direction (d2) around the grid is much smaller than the oscillation in the shorter direction (d1), and so makes little relative difference to the electric field. 3.4.4 The Infinite Channel Now that we have established how to run our simulation and get valid results, we are ready to explore how the pair potential depends on the channel radius R. We still use the configuration shown in Fig. 3.8 and we correct for the finite size effect as just discussed. Water Model 1 We find that for large ion separations, the potential still converges to that of the local model for all different channel widths (Fig. 3.13). We expect the electric field to converge to that of the local model for large ion separations, and to that of water model one without an ion channel for very short ion separations. As Fig. 3.14 shows, this appears to be the case. Finally, we are interested in how the effect of the channel on the electric field (Ewithchan−Ebulkwater) depends on the dielectric model used. Thus we plot Ewithchan−Ebulkwater for various channel widths using the water 1 and local dielectric models (Fig. 3.15). Interestingly, for channel radius 3Å it appears that for ion separations d < 3.5 the electric field is amplified slightly more for the first water model than for the local model, and the reverse is true in the range 3.5 < d < 5, after which the fields of both models appear to converge. For wider channel widths there are no visible non-local effects. Water Model 2 For Water Model 2 we find that the effect of the channel on the electric field (Ewithchan − Ebulkwater) is much more significant than in the local, linear dielectric case (Fig. 3.16). This is still the case for much wider ion channels. Also, for intermediate ion separations the channel modifies the oscillatory nature of the potential such that the location of the local minima and thus of the probable inter-ion distances will be different for ions in an ion channel than for ions in bulk water (Fig. 3.17). 45 Chapter 3. Results -0.003 -0.002 -0.001 0 0.001 0.002 0.003 1 2 3 4 5 6 7 8 U + U c d R=3, Water #1 R=4, Water #1 R=100, Water #1 R=3, Local R=4, Local R=100, Local Figure 3.13: Corrected pair potential U + Uc VS ion separation for various channel radii R using the first water model and the local model. ǫw = 80, L = 20, l = 0.75Å, and a = 0.3Å. All graphs are shifted so that U = 0 for ion separation d = L/2 = 10. Since L = 20, the entire grid is water in the R=100 case. 46 Chapter 3. Results -0.002 -0.001 0 0.001 0.002 0.003 2 3 4 5 6 7 8 E x + E c d R=3, Water #1 R=4, Water #1 R=100, Water #1 R=3, Local R=4, Local R=100, Local 0 0.0001 0.0002 0.0003 0.0004 0.0005 0.0006 2 3 4 5 6 7 8 E x + E c d R=3, Water #1 R=4, Water #1 R=100, Water #1 R=3, Local R=4, Local R=100, Local Figure 3.14: Corrected electric field E +Ec VS ion separation d for various channel radii R, using the first water model and the local model. ǫw = 80, L = 20, l = 0.75Å, and a = 0.3Å. 47 Chapter 3. Results 1e-05 2e-05 3e-05 4e-05 5e-05 6e-05 7e-05 8e-05 9e-05 0.0001 2 4 6 8 10 12 E x ,w ith c ha n - E x ,b ul k w at er ’ d R=3, Water #1 R=4, Water #1 R=3, Local R=4, Local Figure 3.15: Ex,withchan − Ex,bulkwater VS the ion separation d for various channel radii R using the first water model and the local model. For R=4, the results for the local model almost entirely cover those of the first water model. ǫw = 80, L = 20, l = 0.75Å and a = 0.3Å. 48 Chapter 3. Results -0.002 -0.001 0 0.001 0.002 0.003 0.004 0 2 4 6 8 10 12 14 16 E x ,w ith c ha n - E x ,b ul k w at er d R=3, Local R=3, Water #2 R=4, Local R=4, Water #2 -4e-05 -2e-05 0 2e-05 4e-05 6e-05 8e-05 0 2 4 6 8 10 12 14 16 E x ,w ith c ha n - E x ,b ul k w at er d R=6, Local R=6, Water #2 R=9, Local R=9, Water #2 Figure 3.16: Ewithchan−Ebulkwater VS the ion separation d for various chan- nel radii R, using the second water model and the local model. ǫw = 80, L = 35, l = 0.75Å and a = 0.45Å. 49 Chapter 3. Results -0.003 -0.002 -0.001 0 0.001 0.002 0.003 0 2 4 6 8 10 12 14 16 U + U c d R=3, Local R=4, Local R=6, Local R=3, Water #2 R=4, Water #2 R=6, Water #2 Figure 3.17: Corrected Pair Potential U + Uc VS ion separation for various channel radii R for the second water model and the local model. ǫw = 80, L = 35, l = 0.75Å, and a = 0.45Å. 50 Chapter 3. Results 0 1e-05 2e-05 3e-05 4e-05 5e-05 6e-05 7e-05 8e-05 9e-05 1 2 3 4 5 6 7 8 9 10 E x ,c ha n - E x ,b ul k d R=3, Local R=3, Lorentzian #2 R=4, Local R=4, Lorentzian #2 R=6, Local R=6, Lorentzian #2 R=9, Local R=9, Lorentzian #2 Figure 3.18: Ex,withchan − Ex,bulkwater VS the ion separation for various channel radii R using the second Lorentzian model and the local model. ǫw = 80, L = 20, a = 0.45Å and l = 0.75Å. The Lorentzian Models The Lorentzian models do not seem to modify the effect of the ion channel on the electric field from the effects due to the local dielectric model (Fig. 3.18). Discussion The above results show that for the water models, and in particular for the second water model with its more significant non-local effects, the effect of the channel on the pair potential is significantly altered compared to the local case. However, relative to the size of the oscillations for the second water model, these modifications are quite small. Small enough that it is possible that they are a finite grid size correction or a discretization error, however we ran our code several times with different grid sizes L and mesh sizes a and found that, while this had a noticeable affect on the above graphs, it did not change the basic shape of the graph or the conclusions drawn from 51 Chapter 3. Results it. Another possibility is that these results are attributable to the interpo- lation scheme we used for the dielectric profile at the water-cell membrane interface. This is investigated further in the next section on the dielectric barrier of an ion channel. 3.5 The Dielectric Barrier of an Ion Channel In this section we create an ion channel of length less than the dimensions of the grid L, and move an ion from the bulk water outside the channel to inside the channel. To ensure charge neutrality, every vacant node in the lattice has a countercharge of q/(L3 − 1), where q is the charge of the ion. We shall be moving the ion from outside the cell membrane to the center (lengthwise) of the cylindrical ion channel through the cell membrane, as depicted in Fig. 3.19. As discussed in section 1.2.2, the dielectric barrier is determined by tak- ing the minimum over all paths through the ion channel of the maximum self-energy of the ion along each path (and comparing this with the ion’s self-energy in bulk water). However, with non-local electrostatics it is no longer a given that this mini-max occurs at the center of the ion channel. This is why we need to calculate the energy as the ion enters the ion channel, not just at the middle of the channel, and why we must also consider ion positions that are off-center in the y-z plane as well. Once we have found this mini-max, the channel barrier is Ubarrier = Umini−max −Uout, where Uout is the ion’s self energy in bulk water far from the cell membrane. To approximate Uout numerically one could use the ion’s self energy at its initial position in Fig. 3.19 (denoted U(0)), or one could use its self energy in a grid of the same size filled with water only (denoted Ubulkwater). Both of these should result in the calculated dielectric barrier converging to the correct value as the grid size to channel ratio L/M →∞, with the latter giving a greater value for the dielectric barrier than the former for finite values of L/M . We shall use each of these approximations at different times. 3.5.1 The Local Dielectric Approximation As Fig. 3.20 demonstrates, we get the results we would expect for the purely local model as the ion is moved into the ion channel. We next plot the dielectric barrier versus the channel radius (Fig. 3.21), which exhibits decent agreement with our theoretical expectations (eq. 1.13). 52 Chapter 3. Results Figure 3.19: Diagram of the dielectric profile of our grid used to model an ion channel of finite length, and of the trajectory of the ion we are moving. There is still a transition region of length l between the water and the cell membrane, but this has been omitted to unclutter the diagram. 53 Chapter 3. Results 0.00015 0.0002 0.00025 0.0003 0.00035 0.0004 0.00045 0.0005 0.00055 0.0006 0.00065 0.0007 0 2 4 6 8 10 12 U (X ) - U bu lk w at er X R = 2.5, Local R = 3, Local R = 3.5, Local R = 4, Local Figure 3.20: U(X) − Ubulkwater VS the X-position of the single ion present (see diagram 3.19). The local, linear model for polarization is used and the channel radius R is varied. ǫw = 80, L = 24, channel length M = 3R, a = 0.45Å and l = 0.75Å. 54 Chapter 3. Results 0 0.0001 0.0002 0.0003 0.0004 0.0005 0.0006 0.0007 2.6 2.8 3 3.2 3.4 3.6 3.8 4 U B ar rie r R Umax - Ubulk water Umax - U(0) Theor. App. Figure 3.21: Ubarrier = Umax − Uout VS channel radius R, using the local dielectric model. Umax is calculated by taking the maximum of U along the path depicted in Fig. 3.19, and because we are using the local dielectric model we know Umini−max occurs along this path. We use both approxi- mations for Uout (Ubulkwater and U(0)), and the resulting values of Ubarrier sandwich the theoretical prediction eq. (1.13). ǫw = 80, L = 24, M = 3R, a = 0.45Å and l = 0.75Å. 55 Chapter 3. Results 0.0002 0.0003 0.0004 0.0005 0.0006 0.0007 0 2 4 6 8 10 12 U (X ) - U bu lk w at er X U with chan minus U without chan R = 2.5, Local R = 3, Local R = 4, Local R = 2.5, Lorentz #1 R = 3, Lorentz #1 R = 4, Lorentz #1 R = 2.5, Lorentz #2 R = 3, Lorentz #2 R = 4, Lorentz #2 Figure 3.22: Ubarrier = Umax − Ubulkwater VS the X-position of the single ion present (see diagram 3.19). The Lorentzian models for polarization are used and the channel radius R is varied. For the same channel radius, the energies for the local and Lorentzian models match so closely that they are all superimposed on each other. ǫw = 80, L = 24, M = 3R, a = 0.45Å and l = 0.75Å. Given the finite grid-size to channel ratio L/M , the smeared-out counter- charge, and the ambiguity about how to best approximate Uout numerically, exact agreement is not expected. 3.5.2 The Lorentzian Dielectric Approximation As Fig. 3.22 demonstrates, the results for the Lorentzian models are nearly identical to what we would get for the purely local model, but for a slight decrease in the barrier, particularly for narrow channels. We investigate the Lorentzian models no further. 56 Chapter 3. Results -0.003 -0.002 -0.001 0 0.001 0.002 0.003 0.004 0 2 4 6 8 10 12 U (X ) - U bu lk w at er X R = 2.5, Local R = 2.5, Water #2 R = 3, Water #2 R = 3.5, Water #2 R = 4, Water #2 Figure 3.23: Uwithchan−Ubulkwater VS the x-position of the single ion present (see diagram 3.19). The second water model is used, and the channel radius R is varied. ǫw = 80, L = 24, M = 3R, a = 0.45Å and l = 0.75Å. 3.5.3 Water Model 2 The increase in the self energy of an ion as it enters an ion channel is substantially altered from the local result for the second water model (Fig. 3.23). Several interesting features are apparent from this graph: Firstly the potential oscillates as the ion enters the channel, and is not highest in the middle of the channel for narrower channels, suggesting that the potential barrier may not be determined by the potential in the middle of the channel. Secondly, for the narrowest channel there appears to be a potential well in the middle of the channel, meaning that the channel would tend to attract and capture individual ions in dilute ionic solutions. Fig. 3.23 exhibits only the expected minor variations when the mesh size a is reduced and when the grid size L is increased, thus this is not the source of these unusual and unexpected results. The source of these unexpected results appears to lie with the water- membrane interface. To explore these ’edge effects’ the ion was shifted in 57 Chapter 3. Results Figure 3.24: A diagram illustrating the trajectory of the ion we are moving, now with a Y offset. The transition region of length l at the water-cell membrane interface has been omitted to unclutter the diagram. 58 Chapter 3. Results -0.02 -0.01 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0 1 2 3 4 5 6 7 8 9 10 δU X y-offset = 0 y-offset = .45 y-offset = .90 y-offset = 1.35 y-offset = 1.80 y-offset = 2.25 y-offset = 2.70 y-offset = 3.15 Figure 3.25: δU = U(X) − Ubulkwater VS the x-position of the single ion present (see diagram 3.24). The second water model is used, and the y- offset is varied by multiples of the mesh-spacing a. In the case of the last two curves, the ion is actually inside the water-membrane transition region. ǫw = 80, L = 21, R = 3Å, M = 3R, a = 0.45Å and l = 0.75Å. 59 Chapter 3. Results the y-direction by a ’y-offset’ of various multiples of the lattice spacing, before being moved in the x-direction into the ion channel, as shown in Fig. 3.24. This allows us to separate the effects of the narrow channel from the effects of the ions proximity to the water-membrane interface. By varying the y-offset we can move the ion as close as we wish to the interface and even into the water-membrane transition region without changing the channel radius R. Fig. 3.25 shows the resulting energies. The transition region occurs between a y-offset of 2.625 and 3.375 mesh spacings in the y-direction, and so we see that this unexpected potential well is deepest when the ion is still fully inside the ion channel, but very close to the transition region. One possible explanation for why the energy is lowered near the inter- face may lie with the fact that within the cell membrane increasing the polarization (a divergence) is “free”, but maintaining non-zero polarization is energetically costly, whereas in the water-filled ion channel increasing po- larization is costly but maintaining non-zero polarization is cheap because the small κ1 term in the water model. Thus when an ion is situated near the water-membrane boundary, in our model the polarization takes advantage of what is energetically cheap in each medium to lower the overall energy. Specifically, Fig. 3.26 shows how the polarization might increase in the cell membrane, where there is no energy cost for a divergence in the po- larization field, then cross into the water-filled ion channel, where it can minimize the 12(D − P ) 2 + 12κ1P 2 portion of the energy density functional, having already circumvented costs associated with the non-local terms. We shall refer to the kind of polarization field behaviour depicted in Fig. 3.26 as an ’interface effect,’ since it likely results from the implementation of the dielectric profile at the water-membrane interface. To test this theory for why the energy is lowered near the interface, we suppress the kind of polarization field behaviour depicted in Fig. 3.26 in two different ways. First, we set αmembrane = αwater, where previously the gradDiv term αmembrane was zero for the membrane. This makes a diver- gence of the polarization costly in the cell membrane as well, suppressing this ’interface effect’. The result is Fig. 3.27, whose shape is exactly what we would expect for local polarization. However, when we compare the potential of the local and second water models for various different channel radii (Fig. 3.28), we see that, although the shape is quite similar, in scale the potential barrier is far larger for the second water model than for the local model when αmembrane = αwater. The second way in which we suppressed this interface effect is to increase κ1 inside the cell membrane (and thus inside the water-membrane transition 60 Chapter 3. Results Figure 3.26: A diagram with an ion near the edge of the ion channel that illustrates how the polarization may be able to take advantage of what is energetically cheap in each medium to lower the overall energy. 61 Chapter 3. Results 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0 1 2 3 4 5 6 7 8 9 U (X ) - U (0) X y-offset = 0 y-offset = .90 y-offset = 1.35 y-offset = 1.80 y-offset = 2.25 y-offset = 2.70 Figure 3.27: U(x) − U(0) VS the x-position of the single ion present (see diagram 3.24). The second water model is used, amembrane = αwater, and the y-offset is varied by various multiples of the mesh-spacing a. ǫw = 80, L = 18, R = 3Å, M = 3R, a = 0.45Å and l = .01 (a large transition region is no longer necessary). 62 Chapter 3. Results 0.0001 0.0002 0.0003 0.0004 0.0005 0.0006 0.0007 0.0008 0.0009 0.001 0 1 2 3 4 5 6 7 8 9 10 U (X )-U bu lk w at er X R = 2, Local R = 2.5, Local R = 3, Local R = 3.5, Local 0 0.005 0.01 0.015 0.02 0.025 0 1 2 3 4 5 6 7 8 9 10 U (X )-U bu lk w at er X R = 2, Water #2 R = 2.5, Water #2 R = 3, Water #2 R = 3.5, Water #2 Figure 3.28: Uwithchan−Ubulkwater VS the x-position of the single ion present (see diagram 3.19) for various channel radii. The local model is used above and the second water model with αmembrane = αwater (not 0 in the nonlocal case) is used below. ǫw = 80, L = 21, R = 3Å, M = 3R, a = 0.45Å and l = 0.75Å. 63 Chapter 3. Results region as well). In this case we increase κ1 from the default value of 1 (ǫm = 2) to κ1 = 10 (ǫm = 1.1). This strongly supresses all polarization in the cell membrane and throughout the water-membrane transition region. As Fig. 3.29 and 3.28 show, the results are fairly similar for both methods of suppressing this interface effect. Taken together, the results in this section suggest that non-local polar- ization effects may indeed substantially modify the dielectric barrier of a ion channel, but that this modification is highly sensitive to what happens at the water-membrane interface, and we currently have no guidance from molecular dynamics simulations or experiments as to how to model this in- terface. Also, recall that transverse correlations were neglected in equation 1.27 because guidance is needed for how to model them at the interface as well. 64 Chapter 3. Results 0.0002 0.0003 0.0004 0.0005 0.0006 0.0007 0.0008 0.0009 0.001 0.0011 0.0012 0 1 2 3 4 5 6 7 8 9 10 U (X )-U bu lk w at er X R = 2, Local R = 2.5, Local R = 3, Local R = 3.5, Local 0 0.005 0.01 0.015 0.02 0.025 0 1 2 3 4 5 6 7 8 9 10 U (X )-U bu lk w at er X R = 2, Water #2 R = 2.5, Water #2 R = 3, Water #2 R = 3.5, Water #2 Figure 3.29: Uwithchan−Ubulkwater VS the x-position of the single ion present (see Fig. 3.19) for various channel radii and using κ1,membrane = 10 so that ǫm = 1.1, not the default value of ǫm = 2. The local model is used above and second water model is used below. ǫw = 80, L = 21, R = 3Å, M = 3R, a = 0.45Å and l = 0.75Å. 65 Chapter 4 Conclusion In this thesis we have modeled nonlocal polarization effects using an energy density functional for the polarization whose parameters were fit to the dispersion relations obtained from computer simulations. Since predictions vary between computer simulations, we actually implement two different sets of parameter values, which we call the first and second water models, where the second model has more significant and long-range non-local effects. The energy functional is then minimized to find the electric and polariza- tion field configurations, as well as the energy. For comparison, we also used the ’Lorentzian polarization model,’ which simply interpolates between the short-wavelength and long-wavelength limits of the energy functional over a characteristic length-scale of 2π/q0. We try reasonable values of q0 = 2Å and q0 = 3Å, to get a lengthscale of around 2− 3Å, which is about the size of a water molecule. Since only the simplest geometries are tractable analytically, we used an on-grid computational method to calculate this minimum energy configura- tion. The ions were situated on the nodes of the grid, with the displacement and polarization fields situated on the links. Circular boundary conditions were applied at the edges and local minimization moves were applied re- peatedly until the minimum was approached. For the purely local model (with non-local parameters set to 0) we ob- tained excellent agreement with theoretically expected values in bulk wa- ter and in a water filled ion channel through a low-ǫ membrane. For the Lorentzian models and for the first water model we find that the non-local effects are only significant for short ion separations of less than 5Å. For the second water model the non-local effects are very significant and long- ranged, significant on the order of nm. For the second water model we originally found that there was a potential well at the middle of narrow ion channels, however further investigation suggests that this is due to an ’interface effect’ at the water-lipid interface. When this interface effect is suppressed or eliminated, the self energy of an ion moving through an ion channel exhibits a substantially magnified version of the expected behaviour, resulting in a much higher potential barrier. 66 Chapter 4. Conclusion This suggests that an ion’s self energy and thus the dielectric barrier of an ion channel are very sensitive to the implementation of the non-local energy functional at the water-lipid interface, and so guidance from molec- ular dynamics simulations is needed to accurately model this interface. In addition, there are transverse polarization correlations, which we also do not know how to model at the interface. Finally, the non-local polarization energy functionals we have been using to model the dielectric response of water hold in bulk water and are known to hold in ion channels of R ≈ 10Å, however the water confined within a narrow ion channel may not qualify as bulk water, meaning that these energy functionals may not be applicable. Further work with molecular dynamics simulations is needed to investigate these issues. Thus the only conclusion we can draw about the dielectric barrier of an ion channel is that our results suggest it may be significantly modified by non-local dielectric effects, which should therefore certainly not be ignored! 67 Bibliography [1] Kornyshev A. A. and Sutmann G. J. Molecular Liquids, 82:151, 1999. [2] Kornyshev A. A. Bopp P. A. and Sutmann G. Static nonlocal dielectric function of liquid water. Physical Review Letters, 76:1280–1283, 1996. [3] David J. Griffiths. Introduction to Electrodynamics. Prentice-Hall, Inc., third edition, 1999. [4] A.C. Maggs. Simulating nanoscale dielectric response. Physical Review Letters, 96.230603, November 2005. [5] A. A. Kornyshev R. R. Dogonadze, R. Kalman and eds. J. Ulstrup. The Chemica Physics of Solvation. Part A: Theory of Solvation. Elsevier, Amsterdam, 1985. [6] Sofian Teber. Translocation energy of ions in nano-channels of cell mem- branes. J. Stat. Mech., Article Number: P07001, July 2005. 68 Appendix A Verifying Scaling Results and Determining Acceptable Mesh Sizes for the Water Model If we change the mesh spacing a by changing the parameters κ2 and α according to the scaling relation eq. (2.6) and rescale the energy, we expect the energy graphs to be identical, but for differing discretization errors. Since the discretization errors should approach zero as a→ 0, we expect the graphs to converge as a→ 0, which is indeed what we observe in Fig. A.1. It appears that when using larger mesh sizes for the water models, the energy function calculated ‘rounds the corners’ of the actual energy function somewhat. We suspect that this error in the calculated energy occurs be- cause the ratio of the mesh size to the wavelength of maximum susceptibility q0a 2pi is too large, and so the mesh cannot resolve the features of the energy functional correctly. This appears to be confirmed by Fig. A.1, in which the second graph has higher q0 than the first and also requires a smaller mesh size a for convergence than the first. Taken together, the two graphs in Fig. A.1 suggest that for the water model we start to get significant artifacts of the discretization once the mesh spacing a is increased much beyond 1 q0 . For the all of the other models there does not appear to be a very restrictive constraint on a. We expect the energy function to converge to that of local electrostatics at large ion separations, so in this section we test that this does indeed occur and what effect the mesh size a may have on this convergence. As Fig. A.2 shows, this convergence does occur, but much more quickly for larger grid sizes. Thus if we made the mesh size too large we would underestimate the amount of time needed for the energy function to converge to the local energy function. Again this graph also shows that, for mesh sizes that are too large (i.e. a = .8), the energy can deviate so substantially from 69 Appendix A. Verifying Scaling Results and Determining Acceptable Mesh Sizes for the Water Model -0.02 -0.015 -0.01 -0.005 0 0.005 0.01 0.015 0.02 0.5 1 1.5 2 2.5 3 3.5 4 U d a = .25, Local a = .25, Water a = .3, Water a = .4, Water a = .5, Water -0.03 -0.02 -0.01 0 0.01 0.02 0.03 0.5 1 1.5 2 2.5 3 3.5 4 U d a = .25, Local a = .25, Water a = .3, Water a = .4, Water a = .5, Water Figure A.1: Energy VS ion separation for various mesh spacings a, using the water model (eq. 1.27) with A = 3.5 and L = 12Å. In the first graph q0 = 2.5, whereas in the second graph q0 = 4. For comparison, the graph for the local, linear model is shown as well. All energy graphs are shifted by a constant so that each graph’s right-most point has energy zero. 70 Appendix A. Verifying Scaling Results and Determining Acceptable Mesh Sizes for the Water Model -0.0015 -0.001 -0.0005 0 0.0005 0.001 0.0015 2 4 6 8 10 12 U d a = .3, Local a = .3, Water a = .35, Water a = .4, Water a = .5, Water a = .65, Water -4e-05 -3e-05 -2e-05 -1e-05 0 1e-05 2e-05 3e-05 4e-05 7 8 9 10 11 12 U d a = .3, Local a = .3, Water a = .35, Water a = .4, Water a = .5, Water a = .65, Water Figure A.2: Energy VS ion separation for various mesh spacings a. For the water model (eq. 1.27) with q0 = 2.5, A = 3.5, and L = 30Å. For comparison, the energy for the local, linear model is graphed as well. All energy graphs are shifted by a constant so that each graph’s right-most point has energy zero. 71 Appendix A. Verifying Scaling Results and Determining Acceptable Mesh Sizes for the Water Model the correct energy as to make the results entirely useless. 72
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Non local dielectric effects on the dielectric barrier...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Non local dielectric effects on the dielectric barrier of an ion channel Krayenhoff, Willem Bruce 2009
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Non local dielectric effects on the dielectric barrier of an ion channel |
Creator |
Krayenhoff, Willem Bruce |
Publisher | University of British Columbia |
Date Issued | 2009 |
Description | In this thesis we study non-local dielectric models on the sub-nanometer scale. In particular, we focus on the effect of non-local electrostatics on the potential barrier of a water-filled ion channel through a lipid cell membrane, and on the interaction between ions within such an ion channel. A polarization energy functional is used with its parameters calibrated to roughly reproduce the wave-vector-dependent dielectric function ǫ(q) of water. The lipid membrane is still modelled as a local dielectric. This energy functional is discretized onto a lattice and minimized using local moves to find the energy and the electric and polarization field configurations for a given charge distribution and dielectric profile. This method is used to successfully reproduce known theoretical results with and without an ion channel. Necessary, but not necessarily sufficient conditions for obtaining good results are also derived. The dielectric barrier of an ion channel is studied with these non-local dielectric properties of water, and it is found that the results are very sensitive to how the water-membrane interface is modeled. We conclude that more molecular dynamics simulations are needed to provide guidance as to how to implement the polarization energy density functional at this inter- face, as well as whether these energy density functionals, which match the properties of bulk water, need to be modified in order to apply inside narrow ion channels. However, in all instances that we explored the self-energy of the ion traversing the ion channel was substantially modified from its self-energy under the local model, suggesting that non-local dielectric effects may be one significant factor determining the dielectric barrier of an ion channel. |
Extent | 640893 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-03-31 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0067079 |
URI | http://hdl.handle.net/2429/6664 |
Degree |
Master of Science - MSc |
Program |
Physics |
Affiliation |
Science, Faculty of Physics and Astronomy, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2009-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2009_spring_krayenhoff_willem_bruce.pdf [ 625.87kB ]
- Metadata
- JSON: 24-1.0067079.json
- JSON-LD: 24-1.0067079-ld.json
- RDF/XML (Pretty): 24-1.0067079-rdf.xml
- RDF/JSON: 24-1.0067079-rdf.json
- Turtle: 24-1.0067079-turtle.txt
- N-Triples: 24-1.0067079-rdf-ntriples.txt
- Original Record: 24-1.0067079-source.json
- Full Text
- 24-1.0067079-fulltext.txt
- Citation
- 24-1.0067079.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0067079/manifest