Molecular Dynamics simulations of biological ion channels

1. Molecular Dynamics simulations of biological ion channels Titus A. Beu University ”Babeş-Bolyai” Department of Theoretical and Computational Physics Cluj-Napoca, Romania

2. Interest for ion channels • Ion channels - proteins that control the passage of ions across cell membranes. • Responsible for: • generation of action potentials in nerves and muscles • regulation of hormone release from endocrine cells etc. • High selectivity for a particular ion type (Na+, K+, Ca2+, Cl-). • High transport rates (~108 ions per second). • The simplest proteins to which statistical mechanics may be applied. • Highly inhomogeneous electrolyte – challenging aspect. • ”Toolbox” of theoretical models and methods may be validated.

3. The Nobel Prize in Chemistry for 2003:“for discoveries concerning channels in cell membranes.” • Peter Agre:“for the discovery of water channels.” • Roderick MacKinnon:“for structural and mechanistic studies of ion channels.”

4. High-resolution structure of an ion channelD. Doyle, ..., and R. MacKinnon, Science 280, 69-77 (1998). • The KcsA K+ channel • (Streptomyces lividans) • Selectivity filter • High K+ selectivity - adapted to desolvating K ions. • Below and above – fully or partially hydrated K ions. • Inside – binding sites (O atoms) mimic the hydration shell • Gate • Opened by sensor domains

5. Literature Most detailed ion channel model and simulations: • P.S. Crozier et al., Phys. Rev. Lett. 86, 2467 (2001). • P.S. Crozier et al., Biophys. J. 81, 3077 (2001). FFT-accelerated mesh-based Ewald sums (P3M method): • R.W. Hockney and J.W. Eastwood,Computer Simulation Using Particles (IOP, Bristol, 1988). • M. Deserno and C. Holm, J. Chem. Phys. 109, 7678 (1998). • M. Deserno and C. Holm, J. Chem. Phys. 109, 7694 (1998).

6. The model membrane channelP.S. Crozier et al., Phys. Rev. Lett. 86, 2467 (2001). Length = 25 Å, diameter = 10.625 Å, atom-atom distance = 2.5 Å Embedded in a rigid nonpolar membrane (similar to nicotine acetycloline receptor) 388 sites: charges (-0.5e, -0.35e, +0.35e, +0.5e, neutral) + LJ interactions 11 20-atom rings relative rotation 9°

7. 300 K 25 Å 0.02 V/Å 55 Å The simulation cell • The electrolyte – 1M NaCl solution: 600 H2O molecules, 8 Na+ and 8 Cl- • Periodic boundary conditions in all three directions. • Simulation cell: 25 Å x 25 Å x 55 Å

8. Modeling options for water • Molecules made up of atoms subject to constraints • Intermolecular atom-atom potential • Material-point dynamics for atoms • ”Shake”-algorithm - to preserve molecular structure • Rigid molecules: • Intermolecular site-site potential • Rigid-body dynamics for molecules: • Translation of CM – material-point dynamics governed by total force • Rotation about the CM – governed by total torqueQuaternions – alternative to Euler-angles

9. O 0.9572 Å 0.15 Å S 104.52° H H The TIP4P model potential for H2OW.L. Jorgensen et al., Chem. Phys. 79, 926 (1983). • 4 interaction sites:O, H atoms + site S • Electrostatic charges:H atoms + site S • Lennard-Jones interactions:only between O atoms

10. Equations of motion for rigid moleculesMotion of the CM • Total force acting on molecule I: • Newton’s law – equation of motion : • The Verlet propagator – from t to t+∆t:

11. Equations of motion for rigid moleculesDefinition of quaternions • Euler angles – sequence of rotations: • Quaternions – equivalent description – numerically more convenient! • Rotation matrix:

12. Equations of motion for rigid moleculesThe quaternion representation • Angular velocities: • Angular accelerations: Verlet propagator: • Quaternion accelerations:

13. Gaussian thermostat • System in contact with a heat bath (T = const) • Gauss’s principle of least constraint: for the constrained motion • The motion is no longer Newtonian • The constrained equations of motion of the molecules: • The Lagrange multipliers:

14. Periodic boundary conditions Minimizing surface effects for bulk phase • Simulation region replicated to infinity to fill the space. • "Image" particles move solidary with the "real" ones. Correcting particle coordinates • When a particle exits the simulation region – an image enters through the opposite boundary. • Number of particles is conserved. Interactions for short-range potentials • For distances rij > Rcut – interactions can be ignored. • For box size > 2Rcut – particle i lies within Rcut with at most one of all periodic images of another particle j. Minimum image criterion • Among all images of a particle, interactions are considered only with the closest one.

15. The Ewald Sum • Long-range potentials - interactions with distant images cannot be neglected. • Electrostatic energy for all charges and their periodic images: • Slowly decaying - straightforward summation is impracticable. • The trick: split the problem: • f(r)/r - negligible for r > Rcut (rapidly converging in the real space). • [1 - f(r)]/r - slowly varying (rapidly converging in the k-space).Fourier transform can be represented by only a few k vectors.

16. The Ewald SumThe smeared charge density • Equivalent interpretation:Charge density - very sharpperiodic function: • The Fourier representation would never converge • We "smear" the charges - replace them with Gaussian functions: • Smearing parameter a – (inverse length) tunes the relative weight of the real and reciprocal space contributions. • Large a - sharply defined charges.

17. The Ewald SumThe electrostatic energy • The Ewald formula for the electrostatic energy: • Real space contribution: • Reciprocal space contribution: • Self energy correction:

18. The Ewald SumThe electrostatic forces • The Ewald formula for the electrostatic forces: • Real space contribution – interaction of smeared charges: • Reciprocal space contribution – interaction of point and smeared charges:

19. The Ewald SumThe P3M FFT-accelerated method of Hockney and Eastwood • Mesh-based charge density and charge assignment function of order P :

20. The Ewald SumThe P3M FFT-accelerated method of Hockney and Eastwood • Mesh-based charge density and charge assignment function of order P : • Optimal influence function – computed only once: • Mesh-based electrostatic field: • Reciprocal-space contributions to energy and forces:

21. Simulation details1M NaCl solution: 600 H2O molecules, 8 Na+ and 8 Cl- Initial state: Equilibration time: 0.25 ns Time step: 2.5 fs Total simulation time: 200 ns Data storage interval: 2.5 ps

22. Computer run snapshot

23. Electrostatic potential • Uniform electric field - applied in the z-direction to produce the membrane potential: 0.02 V/Å • Potential due to particles and channel sites - Poisson‘s equation for ensemble-averaged mesh-based charge distribution in reciprocal space - using FFT

24. Ion passages • Ion passages are accompanied by an increase of the water polarization, followed by relaxation • Polarizationangle q –between the water dipole and the electric field (channel axis) • Polarization in the channel: • large fluctuations • quick relaxation • reverses sign after ion passage • lower on the average

25. Average density distributions Structured channel – H2O molecules form boundary layers Spatial density profiles depend little on the applied magnetic field

26. Statistics of ion transitions

27. Net current • The magnetic fields • cause aslight increase of the ion current (up to 10%), not a decrease • enhance ion transport indirectly, by enhancing water polarization.

28. Conclusions • The ion channel • structured – water forms boundary layers in the channel and at the membrane walls • high transport rates ~3x108 ions/s (agree with experiments), ion currents ~50 pA • high selectivity – Na+ passages are ~60 times more probable than Cl- passages. • The magnetic fields (1-20 T) and gradients (~100 T/m) technologically available • cause aslight increase of the ion current (up to 10%), not a decrease • enhance ion transport indirectly, by enhancing water polarization. • The channel selectivity is not affected by magnetic fields. • Ion passages cause a pronounced water polarization - importance of water model. • Water polarization in magnetic fields- enhanced in reservoirs, unchanged in the channel. • Cost of “experiment“: • ~1.25 ns/day on up-to-date Compaq workstations or PC (3 GHz) • ~1600 ns simulated  ~1300 CPU days (~3.5 CPU years)