Modeling and simulating complex ionic systems in the presence of dielectric inhomogeneities

[figure: Dezső Boda]

University of Veszprém, Department of Physical Chemistry

[figure]

Dirk Gillespie^{1}, Wolfgang Nonner^{2}, Douglas Henderson^{3}, Bob Eisenberg^{1}

1Rush University Medical Center, Chicago 2University of Miami School of Medicine, Miami 3Brigham Young University, Provo

December 13, 2004, Linz, Austria

Molecular models - two particle system

[figure: u(r_{12}, ω_{1}, ω_{2})] It is usual to divide the model potential into

• repulsive core (hard sphere, exponential, or 1/r^{12} repulsion)

• attractive dispersion force (1/r^{6} or exponential Yukawa attraction)

• elecrostatic forces (Coulomb, dipole, quadrupole,

etc.) Example: SPC/E water

• a Lennard-Jones potential [figure: 4ε [(\frac{σ}{r} )^{12} − (\frac{σ}{r} )^{6}]] on the oxygen

• and partial charges on the atoms

• 6 parameters: 2 LJ parameters, 2 partial charges, distance and angle of hydrogens

Intermolecular pair potential

• depends on the mutual po-

sition (r_{ij})

• and the mutual orientati-

ons (ω_{i}, ω_{j}) of the molecules • It can be calculated from

quantum mechanics or we can use model potentials.

December 13, 2004 Linz, Austria

Molecular models - many particle systems

• The energy is not equal to the sum of the pair potentials: [figure: U ≠ ∑_{i<j} u_{ij}]

• The difference is due to many body forces

• the vicinity of a third particle changes the pair potential

• Solution: using effective pair potentials

• fit the parameters of the potential to experimental results

• For instance: fit vapour-liquid equilibrium curves • How can we relate the microscopic pair potential

to mesaured physical quantities?

[figure]

STATISTICAL MECHANICS

Intermolecular potentials acting between particles

• defines the system on the

microscopic level • determines the energy for

a given configuration of the particles December 13, 2004 Linz, Austria

The relation of the microscopic and the macroscopic level

• What is simulation?

Theory? Experiment?

• Answer: a statistical mechanical tool that gives exact results

• apart from statistical uncertainty and system size problems.

• A theory has to be verified by comparison to simulation.

• A model has to be verified by comparing the simulation results to macroscopic experimental data.

Computer simulations serve as a bridge between real systems and their models December 13, 2004 Linz, Austria

Statistical mechanics 5

Statistical mechanics

• The thermodynamic state has to be specified. For example: canonical ensemble - number of particles (N), volume (V), and temperature (T) are fixed.

• The macroscopic quantities of the system are known if the partition function is known: [figure: Q = ∑_{i} exp(−U_{i}/kT)], where we sum over the possible states of the system and U_{i} is the energy in state i. In the classical theory the sum becomes an integral over phase space (position and velocity).

• For example: Helmholtz free energy: F = −kT ln Q, etc.

• Ensemble average of physical quantity B is [figure: ⟨B⟩ = \frac{1}{Q} ∑_{i} B_{i} exp(−U_{i}/kT)], where B_{i} is the value of B in state i - calculated by Monte Carlo (MC) simulation

• Time average of the physical quantity B is [figure: B̄ = \frac{1}{τ} ∫ B(t)dt]- calculated by molecular dynamics (MD) simulation

• Postulate: the two averages are equal - [figure: ⟨B⟩ = B̄].

December 13, 2004 Linz, Austria Molecular dynamics simulation 6

Molecular dynamics simulation

• We move the particles according to the equations of motions of Newton - a deterministic approach.

• The force acting on each particles is known if the intermolecular potentials are specified.

• Choosing a small Δt time interval, we can calculate the velocities and positions of the particles at time t + Δt if we know the velocities and positions at time t.

• In the simulation, we calculate and store the various physical quantities and we calculate the time average at the end.

• Appropriate to simulate non-equilibrium systems and dynamical properties (e. g. diffusion).

• Problems: native ensemble is NV E (microcanonical).

• Application to other ensembles: sort of mixing of MD and MC.

• Problems with thermostating.

December 13, 2004 Linz, Austria Monte Carlo simulations 7

Monte Carlo simulations

A simplified picture • We randomly sample the phase space - hence the name

• First idea: generate random possible states uniformly and estimate Q = ∑_{i} exp(−U_{i}/kT)

• The ensemble average of B would be ⟨B⟩ = ∑_{i} B_{i}p_{i}.

• where p_{i} = \frac{1}{Q} exp(−U_{i}/kT) is the probability of state i.

• Problem: the number of possible states is huge.

• We sample states that are far from equilibrium: U_{i} is large; exp(−U_{i}/kT) is small.

• Their contribution to Q is small - they are not important.

• Solution: importance sampling - we sample states that are important

Monte Carlo simulations

Importance sampling • Let us generate states with p_{i} probability (Boltzmann sampling).

• We have to bias in the ensemble average with the Boltzmann factor:

[figure] ⟨B⟩ = \frac{∑_{i} B_{i} exp(−U_{i}/kT)/p_{i}}{∑_{i} exp(−U_{i}/kT)/p_{i}} = \frac{∑_{i} B_{i}}{M}

where M is the number of generated states.

• This is a simple mean average.

• Realization of this sampling: Metropolis sampling

N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, E. Teller, Chem. Phys. 1953, 21, 1087. (simulation of hard disks - Los Alamos)

• There are a lot of other possibilities.

Monte Carlo simulations

Metropolis sampling

• Probability of transition i → i + 1 is

[figure] p_{i→i+1} = p_{i+1}/p_{i} = exp[−(U_{i+1} − U_{i})/kT] = exp(−ΔU/kT)

• State i + 1 is accepted by relating to the previous state i - a Markov chain.

• The new state is accepted with the probability [figure: min[1, exp(−ΔU/kT)]].

• The MC step is always accepted if ΔU < 0.

• When ΔU > 0, it is accepted with probability exp(−ΔU/kT)

• This is a Boltzmann sampling.

December 13, 2004 Linz, Austria

Monte Carlo simulations

Other ensembles • NpT: fixed pressure - volume change

• μV T grand canonical: fixed chemical potential - particle insertion/deletion

• Gibbs ensemble - two simulation boxes with particle and volume exchange - simulation of phase equilibria

• Biased moves: to improve the efficiency of sampling and accelerate simulations

• Example: particle exchange between a small volume subsystem (ion channel) and a large volume subsystem (bath).

• Normally, such moves are rare and sampling of occupancy of the channel is inefficient.

Periodic boundary condition

• imitation of macroscopic (infinite) sample

• the Universe is filled with the periodic images of the central simulation cell

• if a particle moved out of the central simulation cell, its periodic image is coming back in from the other side

• surface effects are eliminated • other problems: enforced

periodicity, system size dependence prevails

December 13, 2004 Linz, Austria

Mixed boundary conditions

double layer (DL) geometry: PBC in x and y dimensions, confined by walls in z dimension • PBC is applied in dimensions

where the system is homogeneous

• in the dimension of inhomogeneity the system can be confined by walls

• Inhomogeneous systems: the fluid at interfaces is in equilibrium with a bulk system.

• The bulk is a part of the simulation cell (far from walls) or specified by a chemical potential (grand canonical scheme).

December 13, 2004 Linz, Austria

Problems with large scale MD simulations

• If we know detailed 3D structure of the system (e. g. the KcsA potassium channel on the left), we can perform all atom MD simulation.

• LJ potentials and partial charges on atoms and ions, SPC/E for water.

• For example:

T. W. Allen, S. Kuyucak, S. H. Chung: Molecular dynamics study of the KcsA potassium channel, Biophysical Journal, 1999, 77, 2502.

• Problems. At least, in my view.

D. Doyle et al., Science 1998, 280, 69.

December 13, 2004 Linz, Austria

Problems with large scale MD simulations

• Uncertainties of models and molecular parameters - models and parameters that are good in a given situation might not be good in other situation.

• Example: water - some models are good for thermodynamics, others good for dielectric, interfacial, etc. properties.

• results should be verified by comparison to experimental data (current-voltage).

• long enough simulations cannot be run to get a macroscopic sample

D. Doyle et al., Science 1998, 280, 69.

December 13, 2004 Linz, Austria

Continuum models

Summarized : molecular simulations of systems modelled in full atomic detail might be problematic:

• time consuming - many particles (macroscopic sampling is not possible)

• unknown 3D structures (proteins)

• uncertainties in estimation of potential parameters "Coarse-graining"

• important particles are treated explicitely (e. g. ions)

• the rest is smoothed into a continuum (e. g. solvent, proteins, membrane, electrode)

• the various components with different polarizabilities are modelled as continuums with different dielectric coefficients

• it can be described by an inhomogeneous dielectric coefficient: ε(r) Basic task: solving the Poisson equation [figure: −ε_{0}∇ · [ε (r) ∇ψ (r)] = ρ (r)]

December 13, 2004 Linz, Austria

Macroscopic electrostatics 16

Macroscopic electrostatics

The total charge is split into two parts:

• source (free or explicitely treated) charges: [figure: g(r) = ρ(r)/ε_{0}]and

• polarization charges [figure: h(r) = ρ_{pol}(r)/ε_{0}](to be determined)

• constitutive relation: ρ_{pol}(r) = −∇ · P(r), where

• polarization: P(r) = −ε_{0}χ(r)∇ψ(r)

• dielectric susceptibility: χ(r) = ε(r) − 1

The corresponding Poisson equation: ∇^{2}ψ(r) = −\frac{1}{ε_{0}} [ρ(r) + ρ_{pol}(r)]

The potential is also split into two parts:

•

[figure]

ψ(r) = ψ_{e}(r) + ψ_{i}(r) = ∫_{D} G(r − r′)g(r′)dr′ + ∫_{D} G(r − r′)h(r′)dr′ • the Green-function satisfies ∇^{2}G(r − r′) = −δ(r − r′)

• G(r − r′) = 1/(4π|r − r′|) for the boundary condition ψ(r) → 0 if r → ∞

Variational formalism

Solution of the Poisson equation minimizes the functional

• I [ψ] = \frac{1}{2} ∫_{D} ∇ψ · ∇ψdr − ∫_{D} ψ [g + \frac{1}{2}∇ · (χ∇ψ)] dr

• In minimum the value −ε_{0}I[ψ] is the electrostatic energy

• Substituting the Green function form of ψ(r), the functional can be expressed in terms of the charges: I[g, h]

• At a fixed g(r), we have to find h(r) that minimizes I[g, h].

• This is equivalent to minimizing the h-dependent part of the functional I_{2} [h] for a fixed g.

So, the task is to solve the extremum condition  \frac{δI_{2}[h]}{δh(r)} = 0

•

R. Allen, J.-P. Hansen, and S. Melchionna, Phys. Chem. Chem. Phys. 3, 4177 (2001).

• Allen et al. solved it with the steepest descent method.

Solution of the extremum condition 18

Solution of the extremum condition

\frac{δI_{2}[h]}{δh(r)} = 0

We present a different solution resulting in a matrix equation. The following integral equation can be obtained from the extremum condition:

h(r)ε(r) − ∫_{D} h(r′)∇_{r}ε(r) · ∇_{r}G(r − r′)dr′

= ∇_{r}ε(r) · ∇_{r}ψ_{e}(r) − \frac{1}{ε_{0}} [ε(r) − 1] ρ(r)

• This equation should be discretized and solved for h.

• It is general and valid for arbitrary ε(r) and ρ(r).

• We can make our life easier by making some assumptions for ε(r) and ρ(r).

D. Boda, D. Gillespie, W. Nonner, D. Henderson, B. Eisenberg, Phys. Rev. E 2004, 69, 046702.

December 13, 2004 Linz, Austria

Point charges as source charges

• Suppose that the source charges are point charges: [figure: ρ(r) = ∑_{k} z_{k}eδ(r − r_{k})]

• Polarization charges of magnitude [figure: −z_{k}e\frac{ε(r_{k}) − 1}{ε(r_{k})}δ(r − r_{k})]are induced around the point charges (localized on them).

• The electric potential produced by the point charges and the polarization charges induced on them: [figure: ψ_{e}(r) = \frac{e}{4πε_{0}} ∑_{k} \frac{z_{k}}{ε(r_{k})|r − r_{k}|}] where ε(r_{k}) is the dielectric coefficient at the position of the kth source charge.

• From now on, h(r) denotes only the polarization charges induced on the dielectric interfaces. In this case, the basic equation reduces to : h(r)ε(r) − ∫_{D} h(r′)∇_{r}ε(r) · ∇_{r}G(r − r′)dr′ = ∇_{r}ε(r) · ∇_{r}ψ_{e}(r)

• This is valid for arbitrary ε(r).

December 13, 2004 Linz, Austria Sharp dielectric interfaces 20

Sharp dielectric interfaces

• When the dielectric is piecewise uniform, ε(r) jumps from one value to another along a boundary.

• In this case, the volume integrals over D become surface integrals over the dielectric boundary surfaces B.

• The basic equation then reduces to:

h(s)ε(s) − Δε(s) ∫_{B} h(s′)∇_{s}G(s − s′) · n(s)ds′ = Δε(s)∇ψ_{e}(s) · n(s)

where

• h(s) is now surface charge distribution

• ε(s) is the mean average of the dielectric constants on the two sides of the surface • Δε(s) is their difference, and

• n(s) is a normal vector of the surface

December 13, 2004 Linz, Austria

The matrix equation

The previous equation is valid for any discrete value of s_{α}. Use the notations • normal vector: n_{α} = n(s_{α})

• mean dielectric constant: ε_{α} = ε(s_{α})

• jump in the dielectric constant: Δε_{α} = Δε(s_{α}) The surface B is discretized into surface elements. Let us suppose that the induced surface charge is constant on surface element β denoted by h_{β}. Then our equation becomes

∑_{β} h_{β} [ε_{β}δ_{αβ} − Δε_{α} ∫_{a_{β}} ∇_{s_{α}}G(s_{α} − s_{β}) · n_{α}ds_{β}] = Δε_{α}∇ψ_{e}(s_{α}) · n_{α}.

where the integral expresses the polarization of the surface element β by the induced charge at s_{α}, and vice versa. The integral is taken for the βth surface element.

The matrix equation

If s_{α} is the center of the αth surface element, the previous equation can be written in the matrix form: [figure: Ah = c], where

• the [figure: A_{αβ} = ε_{β}δ_{αβ} − Δε_{α} ∫_{a_{β}} ∇_{s_{α}}G(s_{α} − s_{β}) · n_{α}ds_{β}]matrix depends only on the geometry of the dielectric

• h is the vector of induced charges

• [figure: c_{α} = Δε_{α}∇ψ_{e}(s_{α}) · n_{α}]depends on the source charges

• If the geometry does not change in a simulation, the matrix need be inverted only once at the beginning of the simulation.

• Moving the source charges (the ions) changes the vector c, from which the induced charges are calculated by the matrix-vector multiplication h = A^{−1}c.

Calculation of the integral

Regarding the calculation of the integral I_{αβ} = ∫_{a_{β}} ∇_{s_{α}}G(s_{α} − s_{β}) · n_{α}ds_{β} there are two approaches

1. Treating the surface charge h_{β} as a point charge of magnitude h_{β}a_{β} at the center of the βth surface element of area a_{β}

[figure: I_{αβ} = ∇_{s_{α}}G(s_{α} − s_{β}) · n_{α}a_{β}]

for β ≠ α and 0 otherwise [\underline{Boda et al., Phys. Rev. E, 2004, 69, 046702.}]. This matrix has been developed before from the boundary condition D_{1n} = D_{2n}, E_{1t} = E_{2t}. [\underline{Hoshi et al. JCP, 1987; Lu and Green, PCPS, 1997}]

2. Treating h_{β} as a constant surface charge over the βth surface element, and calculating I_{αβ} numerically after an appropriate parametrization of the surface [\underline{Nonner and Gillespie, Biophys. J. in preparation}]. This approach is especially important if the surface is curved, therefore, we use the name curvature correction.

December 13, 2004 Linz, Austria

Calculation of the energy

• The source charge - source charge interaction energy: [figure: W_{e} = \frac{1}{2} ∑_{j} ez_{j}ψ_{e}(r_{j})]

• The source charge - induced charge interaction energy:

[figure: W_{i} = \frac{e}{8π} ∑_{j} z_{j} ∑_{β} h_{β} ∫_{a_{β}} \frac{ds}{|s − r_{j}|}]

Similarly to the calculation of the integral I_{αβ}, this energy can also be calculated on two levels of approximation:

1. Treating the induced charge as a point charge of magnitude h_{β}a_{β}, the integral in the above equation becomes \frac{a_{β}}{|s_{β} − r_{j}|}.

2. Treating the induced charge as a surface charge with the constant value h_{β} over the βth surface element, the integral should be calculated numerically. This is a \underline{time consuming step in a simulation} because the |s − r_{j}| distances have to be evaluated once the jth ion is displaced.

December 13, 2004 Linz, Austria

Ion is modelled as a hard sphere with a point charge in the center 25

Ion is modelled as a hard sphere with a point charge in the center

• The dielectric constant inside the hard spheres is the same as outside.

• The ions cannot overlap with the dielectric boundaries and cannot leave their host dielectric.

• Otherwise, new dielectric boundaries would appear - matrix A should be inverted again.

One flat interface: ε_{1}|ε_{2} geometry - ions on the right side 26

One flat interface: ε_{1}|ε_{2} geometry - ions on the right side

Comparison of density profiles obtained from the ICC and the image charge method for various resolutions of the grid. [\underline{Boda et al. PRE, 2004, 69, 046702.}]

Two flat interfaces: ε_{1}|ε_{2}|ε_{3} geometry - electrical double layer 27

Two flat interfaces: ε_{1}|ε_{2}|ε_{3} geometry - electrical double layer

It is usual to assume that the inner layer has a dielectric constant different from that of the diffuse layer or the elecrode because of the need for empirical adjustments to obtain agreement with experiments. The elecrostatic consistency of the approach should be studied.

A dielectric sphere ε_{1} = 80 embedded in a dielectric ε_{2} = 2 28

A dielectric sphere ε_{1} = 80 embedded in a dielectric ε_{2} = 2

• \underline{Radius of the sphere}: R = 5 Å.

• \underline{Distance of a source charge e from the center of the sphere}: d = 4 Å

• \underline{Grid}: spherical coordinates ϑ and ϕ evenly divided into intervals Δθ = Δϕ, N is the number of surface elements.

• \underline{Analytic solution}: series in terms of Legendre polynomials.

December 13, 2004 Linz, Austria

A dielectric sphere ε_{1} = 80 embedded in a dielectric ε_{2} = 2 29

The curvature correction is important!!!

December 13, 2004 Linz, Austria

The effective force acting between two charged dielectric spheres

The effective pair potential is defined as [figure: φ(r) = W (r) − W (r → ∞)]. The energy W (r → ∞) corresponds to the interaction energy between the charge and the induced charge on its own sphere. The effective force is [figure: F (r) = −dφ(r)/dr].

MC simulation for a model calcium channel

The selectivity filter of the calcium channel is known to be lined by 4 glutamate groups (EEEE locus). They are represented by 8 half charged oxygen ions that are confined to the filter but can move freely inside. A 0.1M NaCl solution used and a few Ca^{++} ions. Question: whether the Na^{+} or the Ca^{++} enters the filter with a higher probability?

December 13, 2004 Linz, Austria

MC simulation for a model calcium channel

Results for uniform dielectric constant - ε = 78.5 everywhere, no vestibule.

D. Boda et al., Molec. Phys., 2002 100, 2361.

• The concentration of sodium ions is kept constant in the bath at 0.1 M. CaCl_{2} is added to the system.

• The calcium ions replace sodium ions even when their bath concentration is 0.01 times smaller than that of the sodium ions.

• Explanation: charge-space competition (CSC) mechanism

• The calcium ions provide two times the charge to neutralize the glutamate groups than sodium ions using about the same volume.

December 13, 2004 Linz, Austria

Preliminary results for inhomogeneous dielectric

Preliminary results for inhomogeneous dielectric 34 December 13, 2004 Linz, Austria Preliminary results for inhomogeneous dielectric 35 December 13, 2004 Linz, Austria

Problem of absence of charge neutral bulk

The previous figure with a different scale showing the bulk. There is no an electroneutral bulk.

Explanation: • “Fixed charge ensemble”

• Normal procedure in homogeneous ε simulations: the number of ions is chosen so that the system is electroneutral.

• The electroneutral bulk is automatically formed far from the walls (inhomogeneities).

• Here, polarization charges are induced on the dielectric boundaries and prevent the formation of charge neutral bulk.

Simulation by imposing Dirichlet or Neumann boundary conditions

Idea of Nonner and Gillespie (still in preparation): Define confining surface of the system (they call it electrode). Apply an appropriate surface charge σ(r) on the electrode to fulfill BC.

• Dirichlet BC: prescribed potential on the surface: V (s)

• Neumann BC: prescribed normal electric field on the surface: E_{n}(s)

• As the ions move, both the polarization charge h(s) and the electrode charge σ(s) are changing.

• Electrode charge induces polarization charge and vice versa, so their calculation has to be coupled.

Coupling induced and electrode charge calculations 38

Coupling induced and electrode charge calculations

Equation for the dielectric boundaries B:

[figure]

4π \frac{ε(s)}{Δε(s)}h(s) + ∫_{B} h(s′)\frac{(s − s′) · n(s)}{|s − s′|^{3}}ds′ + ∫_{S} σ(s′)\frac{(s − s′) · n(s)}{|s − s′|^{3}}ds′
= − ∑_{k} q_{k} \frac{(s − r_{k}) · n(s)}{ε(r_{k})|s − r_{k}|^{3}}

Equation for the electrode S:

[figure] ∫_{B} \frac{h(s′)}{|s − s′|}ds′ + ∫_{S} \frac{σ(s′)}{|s − s′|}ds′ = 4πε_{0}V (s) − ∑_{k} \frac{q_{k}}{ε(r_{k})|s − r_{k}|}

To solve the problem, both surfaces B and S have to be discretized resulting in a matrix equation. Work is in progress...

December 13, 2004 Linz, Austria

Acknowledgement

• Dr. Mónika Valiskó (University of Veszprém)

• Dr. David Busath (Brigham Young University)

• Tímea Nagy (University of Veszprém)

• Tibor Varga (University of Veszprém)

Download: http://egsys.hu/^{∼}dezo/Citations/citations.html

Thank you for your attention!