Thank you, David !!! Good to visit Cambridge On a ‘typical’ glorious day

1. Thank you, David !!! Good to visit Cambridge On a ‘typical’ glorious day and get away from the snow!

2. How does it work? How do a few atoms control (macroscopic) Biological Function?

3. General Theme Mathematics of Molecular Biology Provides Great Opportunity

4. K+ ~30 Å Ion Channels are Biological Devices Natural nano-valves* for atomic control of biological function Ion channels coordinate contraction of cardiac muscle, allowing the heart to function as a pump Ion channels coordinate contraction in skeletal muscle Ion channels control all electrical activity in cellsIon channels produce signals of the nervous system Ion channels are involved in secretion and absorption in all cells:kidney, intestine, liver, adrenal glands, etc. Ion channels are involved in thousands of diseases and many drugs act on channels Ion channels are proteins whose genes (blueprints) can be manipulated by molecular genetics Ion channels have structures shown by x-ray crystallography in favorable cases *nearly pico-valves: diameter is 400 – 900 picometers

5. Ompf G119D A few atoms make a BIG Difference OmpF 1M/1M G119D 1M/1M OmpF0.05M/0.05M G119D 0.05M/0.05M Glycine replaced by Aspartate Structure determined by Raimund Dutzlerin Tilman Schirmer’s lab Current Voltage relation by John Tang in Bob Eisenberg’s Lab

6. How do a few atoms control (macroscopic) Biological Function? Answer, oversimplified: A few atoms control the electric field Much as they do in transistors

7. + ~3 x 10-9 meters Channels are Devices ValvesandDiodesDifferent Ions Carry Different Signals through Different Channels ompF porin Devices have INPUTS OUTPUTS connected by LAWS involving FLOW from POWER SUPPLIES Analysis of Devices must be NONEQUILIBRIUM with spatially non-uniform BOUNDARY CONDITIONS Thermodynamics, Statistical Mechanics, Molecular Dynamics have No inputs, outputs, flows, or power supplies Power supply=spatially nonuniform inhomogeneous Dirichlet conditions 0.7 10-9 meter = Channel Diameter Flow time scale is 10-4 sec to 1 min Figure of ompF porin by Raimund Dutzler

8. The Electric Field is Strong If you were standing at arm’s length from someone and each of you had  One percent more electrons than protons, the force would lift the Entire Earth! slight paraphrase of third paragraph, p. 1-1 of Feynman, Leighton, Sands (1963)‘The Feynman Lectures on Physics, Mainly Electromagnetism and Matter’ also at http://www.feynmanlectures.caltech.edu/II_toc.html.

9. Channels are Devices Channels are (nano) valves Valves Control Flow Classical Theory & Simulations NOT designed for flow Thermodynamics, Statistical Mechanics do not allow flow Rate and Markov Models do not Conserve Current

10. Maxwell’s Equationsare aboutConservation of Current Conservation of Charge

11. Valves and Devices haveINPUTS OUTPUTS connected by LAWS involving FLOW fromPOWER SUPPLIES so Analysis of Devices must be NONEQUILIBRIUM with spatially non-uniform BOUNDARY CONDITIONS Thermodynamics, Statistical Mechanics, Molecular Dynamics have No inputs, outputs, flows, or power supplies i.e., Power Supply = spatially nonuniform inhomogeneous Dirichlet conditions

12. Thermodynamics, Statistical Mechanics, Molecular Dynamics are UNSUITED for DEVICES Thermodynamics, Statistical Mechanics, Molecular Dynamicshave No inputs, outputs, flows, or power supplies Power supply = spatially nonuniform inhomogeneous Dirichlet conditions Analysis of Devices must be NONEQUILIBRIUM with spatially non-uniform BOUNDARY CONDITIONS

13. Where to start? A different talk! Compute all atoms in a device? Calibrated all-atom simulations are Barely Feasible if they must accurately compute biological function Macroscopic Time & Distance Scales Macroscopic Electric Fields & Gradients Power SuppliesPower Supply = spatially nonuniform inhomogeneous Dirichlet Boundary Conditions Flows, Non-ideal mixtures including [Ca2+]=[10-8↔101 M]

14. Where to start? Biological Adaptation Crowded Charge

15. Active Sites of Proteins are Very Charged 7 charges ~ 20M net charge = 1.2×1022 cm-3 liquidWater is 55 Msolid NaCl is 37 M + + + + + - - - - Selectivity Filters and Gates of Ion Channels are Active Sites Physical basis of function OmpF Porin Hard Spheres Na+ Ions are Crowded K+ Ca2+ Na+ Induced Fit of Side Chains K+ 4 Å Figure adapted from Tilman Schirmer

16. Working Hypothesis Crucial Biological Adaptation is Crowded Ions and Side Chains A different talk: Biology occurs in concentrated >0.3 M mixtures of spherical charges NOT IDEAL AT ALL Solutions are extraordinarily concentrated >10M where they are most important, near DNA, enzyme active sites, and channels and electrodes of batteries and electrochemical cells. (Solid NaCl is 37M) Ideal Gas, zFV+RT ln C, fits NOTHING!!

17. Crowded Active Sitesin 573 Enzymes Jimenez-Morales,Liang, Eisenberg

18. Don’t worry! Crowded Charge is GOOD It enables SIMPLIFICATION by exploiting a biological fact (an adaptation) Charges are Crowded where they are important!

19. Where do we begin? Crowded Charge enables Dimensional Reduction* to a Device Equation Essence of Engineering is knowing What Variables to Ignore! WC Randels in Warner IEEE Trans CT 48:2457 (2001) *Dimensional reduction = ignoring some variables

20. Where do we begin? • Crowded Charge • has • HUGE electric fields • Poisson Equation, i.e., Conservation of Charge • and • LARGE steric repulsionFermi distribution • Stay tuned: • for a Poisson-Fermi treatment

21. The Electric Field is Strong If you were standing at arm’s length from someone and each of you had  One percent more electrons than protons, the force would lift the Entire Earth! slight paraphrase of third paragraph, p. 1-1 of Feynman, Leighton, Sands (1963)‘The Feynman Lectures on Physics, Mainly Electromagnetism and Matter’ also at http://www.feynmanlectures.caltech.edu/II_toc.html.

22. Where do we begin? Crowded Charge has HUGE electric fields Poisson Equation, i.e., Conservation of Charge Forces and Potential ALWAYS depend on charge

23. Motivation and Assumption for Fermi-Poisson Largest Effect of Crowded Charge is Saturation Saturation cannot be described at all by classical Poisson Boltzmann approach and is described in a uncalibrated way by present day Molecular Dynamics when Mixtures and Divalents are Biologically Important in Concentrations of 10-8 to 101 M • Simulating saturationby interatomic repulsion (Lennard Jones) is a singular mathematical challengeto be side-stepped if possible, • particularly in three dimensions, Eisenberg, Hyon and Liu (2010) JChemPhys 133: 104104

24. Motivation Natural Description of Crowded Charge is a Fermi Distribution because it describes Saturationin a simple way used throughout Physics and Biophysics, where it has a different name! Simulating saturationby interatomic repulsion (Lennard Jones) is a significant mathematical challengeto be side-stepped if possibleEisenberg, Hyon and Liu (2010). JChemPhys 133: 104104

25. Fermi Description usesEntropy of Mixture of Spheresfrom Combinatoric Analysis W is the mixing entropy of UNEQUAL spheres with N available NON-UNIFORM sites Connection to volumes of spheres and voids, and other details are published in 5 papers J Comp Phys (2013) 247:88 J Phys Chem B (2013) 117:12051J Chem Phys (2014) 141: 075102 J Chem Phys, (2014) 141: 22D532 Physical Review E (2015) 92:012711 Expressions in other literature are not consistent with this entropy

26. A Nonlocal Poisson-Fermi Model for Electrolyte Solutions Jinn Liang Liu 劉晉良 Jinn-Liang is first author on our papers • J Comp Phys (2013) 247:88 J Phys Chem B (2013) 117:12051J Chem Phys (2014) 141: 075102 J Chem Phys, (2014) 141: 22D532 Physical Review E (2015) 92: 012711Chem Phys Letters (2015)637: 1J Phys Chem B (2016) 120: 2658

27. Voids are Needed (Electro)Chemical Potentialand Void Volume It is impossible to treat all ions and water molecules as hard spheres and at the same time have Zero Volume of interstitial Voidsbetween all particles

28. Consistent Fermi Approach is NovelConsistent Fermi approach has not been previously applied to ionic solutionsas far as we, colleagues, referees, and editors know Previous treatments* have inconsistent treatment of particle size They do not reduce to Boltzmann functionals in the appropriate limit Previous treatments often do not include non-uniform particle size Previous treatments* are inconsistent with electrodynamics and nonequilibrium flows including convection Details Previous treatments do not include discrete water or voids. They cannot deal with volume changes of channels, or pressure/volume in general Previous treatments do not include polarizable water with polarization as an output *Previous treatments Bazant, Storey & Kornyshev,. Physical Review Letters, 2011. 106(4): p. 046102. Borukhov, Andelman & Orland, Physical Review Letters, 1997. 79(3): p. 435. Li, B. SIAM Journal on Mathematical Analysis, 2009. 40(6): p. 2536-2566. Liu, J.-L., Journal of Computational Physics 2013. 247(0): p. 88-99. Lu & Zhou, Biophysical Journal, 2011. 100(10): p. 2475-2485. Qiao, Tu & Lu, J Chem Phys, 2014. 140(17):174102 Silalahi, Boschitsch, Harris & Fenley, JCCT2010. 6(12): p. 3631-3639. Zhou, Wang & Li Physical Review E, 2011. 84(2): p. 021901.

29. ChallengeCan Simplest Fermi Approach • Describe ion channel selectivity and permeation? • Describe non-ideal properties of bulk solutions? There are no shortage of chemical complexities to include, if needed! Classical Treatments of Chemical Complexities

30. Evidence (start)

31. Poisson Fermi Approach toBulk Solutions Same Fermi Poisson Equations, different model of nearby atoms in Hydration Shells

32. Bulk Solution How well does the Poisson Fermi Approach for Bulk Solutions? Same equations, different model of nearby atoms Occupancy is 6 + 12 Waters*held Constant in Model of Bulk Solution in this oversimplified Poisson Fermi Model Liu & Eisenberg (2015) Chem Phys Ltr 10.1016/j.cplett.2015.06.079 *in two shells: experimental Data on OccupancyRudolph & Irmer, Dalton Trans. (2013) 42, 3919 Mähler & Persson, Inorg. Chem. (2011) 51, 425

33. Activity CoefficientsNa+ Cl-‘normalized’ free energy per mole

34. Activity CoefficientsCa2+ Cl2¯‘normalized’ free energy per mole

35. Gramicidin A Unusual SMALL Bacterial Channel often simulated and studied Margaret Thatcher, student of Nobelist Dorothy Hodgkin Bonnie Wallace leading worker Validation of PNP Solvers with Exact Solution following the lead of Zheng, Chen & Wei J. Comp. Phys. (2011) 230: 5239.

36. Three Dimensional TheoryComparison with Experiments Gramicidin A

37. Steric Effect is Large in (crowded) Gramicidin PNPF spheresvs PNP points Points Water Occupancy Spheres Current vsVoltage Spheres Points K+ Occupancy Spheres Points Three Dimensional Calculation Starting with Actual Structure

38. Cardiac Calcium Channel CaV.n Binding Curve Lipkind-Fozzard Model Liu & Eisenberg J Chem Phys 141(22): 22D532

39. Signature of Cardiac Calcium Channel CaV1.n Anomalous* Mole Fraction (non-equilibrium) Na Channel Ca Channel *Anomalous because CALCIUM CHANNEL IS A SODIUM CHANNEL at [CaCl2] 10-3.4 Ca2+ is conducted for [Ca2+] > 10-3.4, but Na+ is conducted for [Ca2+] <10-3. Liu & Eisenberg (2015) Physical Review E 92: 012711

40. More Detail COMPUTING FLOW

41. Poisson-Fermi Analysis is NON-Equilibrium Flows are Essential in Devices & Biology Structure is Essential in Devices & Biology Implemented fully in 3D Code to accommodate 3D Protein Structures Flows cease only at death • PNPF uses treatment by Santangelo 20061 used byKornyshev 20112 of near/far fields crudely separated by fixed correlation length • PNPF introducessteric potential3,4 so unequal spheres are dealt with consistently • PNPF force equation reduces3,4 to pair of 2nd order PDE’s and Appropriate boundary conditions that are consistent and allow Robust and Efficient Numerical Evaluation 4) PNPF combines Force Equation and Nernst-Planck Description of Flow 1PhysRev E (2006)73:041512 2PhysRev Ltrs (2011) 106:046102 3JCompPhys (2013) 247:88 4J PhysChem B (2013) 117:12051

42. What is PNPF? PNPF = Poisson-Nernst-Planck-Fermi Implemented fully in 3D Code to accommodate 3D Protein Structures Flow Three Dimensional computation is facilitated by using 2nd order equations Force approximates dielectric of entire bulk solution including correlated motions of ions, following Santangelo 20061 with Liu’s corrected and consistent Fermi treatment of spheres.2,3,4 We introduce3,4two second order equations and boundary conditions That give the polarization charge density 1PhysRev E (2006)73:041512 2PhysRev Ltrs (2011) 106:046102 3JCompPhys (2013) 247:88 4J PhysChem B (2013) 117:12051

43. Computational Problems Abound and are Limitingif goal is to fit real data Scientists must grasp, …….not just reach, if we want devices to work and models to be transferrable It is very easy to get results that only seem to converge, and are in fact Not Adequate approximations to the converged solutions Jerome, J. (1995) Analysis of Charge Transport. Mathematical Theory and Approximation of Semiconductor Models. New York, Springer-Verlag. Markowich, P. A., C. A. Ringhofer and C. Schmeiser (1990). Semiconductor Equations. New York, Springer-Verlag. Bank, R. E., D. J. Rose and W. Fichtner (1983). Numerical Methods for Semiconductor Device Simulation IEEE Trans. on Electron Devices ED-30(9): 1031-1041. Bank, R, J Burgler, W Coughran, Jr., W Fichtner, R Smith (1990) Recent Progress Algorithms for Semiconductor Device Simulation Intl Ser Num Math 93: 125-140. Kerkhoven, T. (1988) On the effectiveness of Gummel's method SIAM J. Sci. & Stat. Comp. 9: 48-60. Kerkhoven, T and J Jerome (1990). "L(infinity) stability of finite element approximations to elliptic gradient equations." Numer. Math. 57: 561-575.

44. Computational Electronics has solved these problems over the last 40 years in thousands of papers used to design our digital devices Devices and calculations work Models are transferrable Vasileska, D, S Goodnick, G Klimeck (2010) Computational Electronics: Semiclassical and Quantum Device Modeling and Simulation. NY, CRC Press. Selberherr, S. (1984). Analysis and Simulation of Semiconductor Devices. New York, Springer-Verlag. Jacoboni, C. and P. Lugli (1989). The Monte Carlo Method for Semiconductor Device Simulation. New York, Springer Verlag. Hess, K. (1991). Monte Carlo Device Simulation: Full Band and Beyond. Boston, MA USA, Kluwer. Hess, K., J. Leburton, U.Ravaioli (1991). Computational Electronics: Semiconductor Transport and Device Simulation. Boston, Kluwer. Ferry, D. K. (2000). Semiconductor Transport. New York, Taylor and Francis. Hess, K. (2000). Advanced Theory of Semiconductor Devices. New York, IEEE Press. Ferry, D. K., S. M. Goodnick and J. Bird (2009). Transport in Nanostructures. New York, Cambridge University Press. It is very easy to get results that only seem to converge, and are in fact not adequate approximations to the converged solutions. Jerome, J. W. (1995). Analysis of Charge Transport. Mathematical Theory and Approximation of Semiconductor Models. New York, Springer-Verlag.

45. Keys to Successful Computation • Avoid errors by checking against analytical solutions of Guowei and collaborators • Avoid singularities (i.e., acid/base charges) on protein boundaries that wreck convergence • Use a simplified Matched Interface Boundary sMIB method of Guowei and collaborators modified to embed Scharfetter Gummel SG criteria of computational electronics (extended to include steric effects). Scharfetter Gummel is REQUIREDto ENSURE CONTINUITY OF CURRENT Charge Conservation is not enough Scharfetter and Gummel, IEEE Trans. Elec. Dev.16, 64 (1969) P. Markowich, et al, IEEE Trans. Elec. Dev. 30, 1165 (1983). Zheng, Chen, and G.-W. Wei, J. Comp. Phys. 230, 5239 (2011). Geng, S. Yu, and G.-W. Wei, J. Chem. Phys. 127, 114106 (2007). S. M. Hou and X.-D. Liu, J. Comput. Phys. 202, 411 (2005). J.-L. Liu, J. Comp. Phys. 247, 88 (2013). • Modified Successive Over-relaxation SOR for fourth order PNPF

46. More Detail INSIDE CHANNELS

47. Steric Effect is Significant Gramicidin is Crowded Shielding is Substantial Electric Potential Steric Potential Shielding Shielding has been ignored in many papers, where Results are often at one concentration or unspecified concentration,as in most molecular dynamics Channel is often described as a potential profileThis is inconsistent with electrodynamicsas in classical rate models Shielding

48. GramicidinTwo K+ Binding Sites OUTPUTS of our calculations Binding sites are prominent in NMR measurements & MD calculations BUT they VARY with conditions in any consistent model and so cannot be assumed to be of fixed size or location

49. Inside Gramicidin Water Density Dielectric Function an OUTPUT of model Liu & Eisenberg J Chem Phys 141: 22D532

50. Inside the Cardiac Calcium ChannelCaV1.n Dielectric Function An Output of this Model Water Density Liu & Eisenberg (2015) Phys Rev E 92: 012711 Liu & Eisenberg J Chem Phys 141(22): 22D532