University of Southern California Dissertations and Theses

Capturing the complexity of ion channels: simulations of long-time dynamics, conformational changes and the effect of the membrane potential

(USC Thesis Other)

Capturing the complexity of ion channels: simulations of long-time dynamics, conformational changes and the effect of the membrane potential

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content CAPTURING THE COMPLEXITY OF ION CHANNELS: SIMULATIONS OF LONG-TIME DYNAMICS, CONFORMATIONAL CHANGES AND THE EFFECT OF THE MEMBRANE POTENTIAL by Anatoly Dryga A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulllment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (CHEMISTRY) May 2012 Copyright 2012 Anatoly Dryga Acknowledgements I would like to thank Professor Arieh Warshel for leading this work and for his active participation and help during all phases of my graduate studies. Prof. Warshel's passion for science and his intuition have helped me to understand how to conduct relevant research. I also want to thank Professor Aiichiro Nakano who introduced me to the world of high-performance computing and shared his vision of how to write robust and scalable scientic programs. Professor Nakano participated in my screening and qualication exams and served on my defense committee. He also graded my nu- merous projects and assignments. Of course many past and current group members of Prof. Warshel's lab helped me enormously during my struggles with science and programming. I greatly en- joyed numerous talks about all aspects of our lives. I would like specically to thank Prof. Kamerlin and Dr. Pisliakov for help in the very beginning of my grad- uate work and for showing me what good team work means. Special thanks to Dr. Burykin and Dr. Kato for their work on ion channels, which was the basis for all my work. Dr. Vicatos and Anna Rychkova helped me with coarse-grained model. And of course our senior programmer, Dr. Chu, solved numerous tasks in coding Molaris. Thank you all! And nally I would like to thank my family for their love and support. Kate, Max and Alex keep me busy all the time! I am very glad to have you in my life. ii Table of Contents Acknowledgements ii List of Tables v List of Figures x Abstract xi Chapter 1. Introduction 1 1.1 Ion channels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Molecular Structures . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Simulation of Ion Channels . . . . . . . . . . . . . . . . . . . . . . . 5 Chapter 2. Multi-Scale Simulations 6 2.1 General Consideration . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.2 Coarse-Graining and Dynamics . . . . . . . . . . . . . . . . . . . . 7 Chapter 3. Long-time Simulations 9 3.1 Dynamics for Microsecond Simulations . . . . . . . . . . . . . . . . 9 3.2 Gramicidin A: Method . . . . . . . . . . . . . . . . . . . . . . . . . 13 3.3 Gramicidin A: Results . . . . . . . . . . . . . . . . . . . . . . . . . 17 3.4 Barrier Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 Chapter 4. Conformational Changes 35 Chapter 5. Voltage-activated Ion Channel 38 5.1 Motivation and System . . . . . . . . . . . . . . . . . . . . . . . . . 38 5.2 Coarse-Grained Model for the Channel . . . . . . . . . . . . . . . . 42 5.3 Eect of Electrolyte and Membrane Potential . . . . . . . . . . . . 56 5.4 Model Systems and Validation . . . . . . . . . . . . . . . . . . . . . 74 iii 5.5 Simulation of Kv1.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 5.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 Chapter 6. Conclusions 95 6.1 Main Achievements . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 6.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 Bibliography 99 iv List of Tables Table 5.1 Parameters for the G vdw side term. In the table 1 is used for interactions involving residues of dierent polarity whereas 1 + 2 is used for interactions involving residues of the same polarity as well as interactions with ionized residues . . . . . 52 Table 5.2 Parameters for the U phipsi mm term. The angular values in the table are given in degree. . . . . . . . . . . . . . . . . . . . . 53 Table 5.3 Self Energy Parameters for the Ionizable Residues . . . . . . 53 Table 5.4 Self Energy Parameters for the Polar Residues. Equations identical to 5.10, 5.11 and 5.14 were used for polar residues (but the parameters used are dierent and listed in the table) 54 Table 5.5 Self Energy Parameters for the Hydrophobic Residues from polar and membrane neighbors. Equations identical to 5.11 and 5.14 were used for the hydrophobic (but the parameters used are dierent and listed in the table) . . . . . . . . . . . 54 Table 5.6 The parameters used in equation 5.12 for calculation of the number of neighboring residues. d mem refers to the membrane grid spacing. . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 Table 5.7 Parameters for Hydrophobic energy term . . . . . . . . . . . 55 Table 5.8 The parameters for the hydrogen bond function . . . . . . . 55 Table 5.9 The partial charges of ionized residues in a monomer of the Kv1.2 channel (only residues with dierent charges for dier- ent conformations are included ) . . . . . . . . . . . . . . . . 69 Table 5.10 The free energy at dierent regions of the conformation/voltage landscape (The eect of the change in potential in each struc- ture was adjusted, subtracting the capacity work in the open structure from the energy of all the other structures.) . . . . 88 v List of Figures Figure 1.1 Molecular representation of Kv1.2 channel (side view). . . . 3 Figure 1.2 Molecular representation of Kv1.2 channel (top view). . . . 4 Figure 3.1 Demonstrating the nature of response to a 1-D force in a 2- D system. The Figure is aimed at clarifying the fact that nonequilibrium type approaches can never provide reliable PMF with short simulations and large forces. In such case, the system will pass through path 3 and will not collect infor- mation about the actual least energy path. Of course running an innite number of short trajectories with strong force will also rarely sample path 1, but this would not oer an ecient sampling strategy. . . . . . . . . . . . . . . . . . . . . . . . 12 Figure 3.2 Benchmark system of a truncated gR with a sodium ion. . . 16 Figure 3.3 PMFs for the high barrier model (from ref 5) and the low barrier (9.3 kcal/mol) model of Figure 3.4 with charges q = -0.3. The PMF of the system with the low barrier is shifted to have the same minimum as in the system with the high barrier. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 Figure 3.4 Truncated gR with four additional charges at the center of the system. The change of the charges provides a systematical way to alter the original PMF. . . . . . . . . . . . . . . . . 19 Figure 3.5 Two trajectories obtained with the explicit model for dierent pulling forces, K. The simulations were done for the high barrier case, and K is given in kcal/mol/Angstrom 2 . . . . 20 Figure 3.6 Typical MD and LD ion trajectories for the explicit and im- plicit models, respectively (A), as well as the velocity auto- correlation functions (B) for the high barrier model and for K = 0.0001 kcal/ mol/Angstrom 2 . . . . . . . . . . . . . . 21 vi Figure 3.7 First passage time for explicit and implicit models for dier- ent force constants. The data shown are for two models: A high barrier model (from ref [51]) and a low barrier model (truncated gR with four negative charges). . . . . . . . . . 22 Figure 3.8 Long time trajectories of the explicit and implicit models for the system with a low barrier (9.3 kcal/mol) and without a pulling force. The Figure demonstrates that we can repro- duce the long time behavior of the explicit model. . . . . . 24 Figure 3.9 Original PMF for the explicit model and the modied PMFs for the system with four additional positive charges. g 6= is dened as a dierence in the barrier height of the system with four xed charges and the high barrier system (without four xed charges). . . . . . . . . . . . . . . . . . . . . . . 25 Figure 3.10 Relative rst passage time as a function of the g 6= (dened in Figure 3.9) of the explicit model for three dierent pulling forces. The estimated g 6= is independent of the magnitude of the pulling force. The relative time is dened as the ratio of the rst passage time (with a given pulling force) for the system with a modied PMF to the rst passage time for the original high barrier model. . . . . . . . . . . . . . . . . . . 27 Figure 3.11 Comparison of the g 6= obtained by the renormalization approach and the exact results obtained with the careful anal- ysis of the explicit model. . . . . . . . . . . . . . . . . . . . 28 Figure 3.12 Comparison of the PMF of the explicit model (black) and the PMF (red) obtained by using the renormalization approach (using LD simulations that reproduce the rst passage time for dierent segments of the explicit model) as well as esti- mation of the barrier height. . . . . . . . . . . . . . . . . . 30 Figure 4.1 Qualitative renormalization study of the conformational change in ADK for the full-atom, coarse-Grain, and 2-D LD models. It describes the timed appendance of the response to force of the three models. . . . . . . . . . . . . . . . . . . . . . . . 37 vii Figure 5.1 Illustrating the free energy balance in a voltage-activated ion channels. (A) The gure considers schematically the activa- tion process, where the protein moves from the closed to open state due to the allosteric eect of the external potential (V). When the external potential is negative (V < 0), the closed conformation is more stable and the channel is blocked. Ap- plication of a positive potential (V > 0) stabilizes the open conformation and the system moves to its ion conduction state. (B) A tentative landscape for activation of the Kv1.2 system in the two-dimensional space dened by the protein closed to open conformational coordinate and the external potential. Note that the potential coordinate is a coordinate of an external constraint and as such the motion from nega- tive to open potential does not have to be downhill (because it does not include the external power source). Note also that the surface was generated with empirical valence bond-type parabolas and from the data of free energy table. Thus the increase in energy, upon moving from the minima toward the boundary of the map, does not re ect the real surface. . . . 40 Figure 5.2 From full atom representation to Coarse Grained representa- tion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 Figure 5.3 A schematic description of the treatment of the simulation system. The regions considered are (i) the protein membrane system (region I); (ii) the region of the solution and the ions that is treated with explicit ions and Langevin dynamics sim- ulation (region II), which is not shown because this region is optional and is not being used in the present study; (iii) the region with explicit grid points (region III); (iv) the bulk re- gion. The region between the bulk and the electrodes (region IV). The indices L and R stand for the left and right sides, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . 57 Figure 5.4 A schematic representation of the membrane/electrodes sys- tem (the membrane can also contain a protein). . . . . . . 63 viii Figure 5.5 Clarifying the considerations that led to the treatment ofV i ext . The gure describes a division of the distance between the grid pointZ start pro (where we have the "ideal potential") to the position of the given residue into two steps, which are then being used to evaluate the potential on the residues using a distance-dependent dielectric. . . . . . . . . . . . . . . . . 66 Figure 5.6 The contribution from the crucial Arg residues to the ener- getics of the activation of the Kv1.2 channel. . . . . . . . . 72 Figure 5.7 The potential around a charge obtained by our simulation and by the Debye Huckel (DH) theory as well as Coulombic (unscreened) potential. . . . . . . . . . . . . . . . . . . . . 75 Figure 5.8 The GouyChapman potential. . . . . . . . . . . . . . . . . 76 Figure 5.9 (a) The electrolyte potential (in Volt) and (b) the electrolyte charge distribution along Z-axis for a system with two elec- trodes (external potential 100 mV). We have shown the re- sults for two dierent electrolyte concentrations: 10 mM (red) and 100 mM (blue). Fig. 5(a) also shows the linear elec- trode potential (in absence of electrolyte) as reference (black). These calculations have been performed with (17X17) peri- odic images in the XY plane (see Fig. 7 for convergence study). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 Figure 5.10 (a) The electrostatic potential (in kcal/mol) and (b) elec- trolyte charge distribution along Z-axis for a system with two electrodes (external potential 100 mV), 40 A thick mem- brane at the center of the box (marked by red vertical lines) and electrolyte with concentration 100 mM. The present case required the use of (17x17) periodic images (in the XY plane) to reach convergence. For this system, the computed capaci- tance is 2:2 10 3 CV 1 m 2 and the analytical capacitance is 4:3 10 3 CV 1 m 2 . . . . . . . . . . . . . . . . . . . . . . 79 Figure 5.11 A simplied model for exploring the gating charges. The model involves a membrane and positive charges that are moved from one side to the other. . . . . . . . . . . . . . . 80 ix Figure 5.12 The accumulative integration of the charge proles for the model system described in Fig. 8 when the+16e charge re- sides at Z=+10 (black) and Z=-10 (red). Here we dene the gating charge as the change in the accumulated electrolyte charges (the dierence of the maxima in this gure), on ei- ther side of the membrane, upon movement of the (protein) charge inside the membrane (see Eq. 42). The gating charge for this system is found to be 5. . . . . . . . . . . . . . . . 81 Figure 5.13 The potential in the Kv1.2 system in the open and closed states. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 Figure 5.14 The gating charge for the Kv1.2 system. The gure presents the integrated charge . The dierence in the maximum ac- cumulated charge between two states is taken as the gating charge, which is approximately 10 in this gure. . . . . . . 86 Figure 5.15 Simulating the gating uctuations. The gure describes the results of the Langevin dynamics simulations on the tentative surface of Fig. 5.1 forV = 0. The simulations were done with a friction of 150 ps 1 . . . . . . . . . . . . . . . . . . . . . . 90 Figure 5.16 A conceptual gure of the activation of the Kv1.2 and related channels. The gure considers the allosteric eect of the po- tential, where the change of position of the positive gating charges reduces the repulsion between theses charges and the positive elements of the gating region, allowing the attached nonpolar region to open and to reduce the desolvation penalty of a transferred ion. . . . . . . . . . . . . . . . . . . . . . . 91 x Abstract The relationship between membrane potential and the gating of voltage-activated ion channels, as well as ion current through pores, is currently a problem of great interest. For example, simulations of conformational changes during the voltage- induced transition and ion current represent a major challenge that cannot be over- come at present by brute-force molecular dynamics approaches. To progress on this front, we developed and validated a renormalization approach and extended our previously developed coarse-grained (CG) model to describe the eect of the membrane potential and electrolyte solutions. Here we apply our model in realis- tic simulations of the energetics of the Kv1.2 channel. The simulations reproduce the observed experimental trend for both the gating charge and the eect of the membrane potential on the stability of the closed and open conformations. This in- dicates that our CG modeling oers a powerful tool for simulating voltage-activated channels. xi Chapter 1. Introduction 1.1 Ion channels The smallest biological object with all properties required for life is a cell. The membrane was developed in the early stages of evolution as a mechanism to create conditions inside the cell that are dierent from those outside the cell. Phospholid membranes solve the problem of isolating the components of the cell from the external world, but because of the need in living species to exchange energy and components with the rest of the world, special molecules were needed to facilitate these processes. Ion channels are one of the important classes of molecules for the exchange of ions and the distribution of information through the multicellular organism. These complex proteins allow selective movement of ions [103, 8, 89, 52, 3] through the membrane. Ions move from regions in which they are highly concentrated into regions of lower concentration thus it is a spontaneous manner. The selective movement of ions is also regulated by the gating part of the ion channels. Most of the channels can exist at least in two dierent conformations: one conformation allows ions to pass through pore while the other blocks ion movement. The transition between distinct conformation states is due to some external stimuli, e.g. ligand binding, changes in the membrane potential, etc. Here we concentrate on voltage-gated ion channels. 1 This classss of ion channels changes its conformation state due to changes in the membrane potential. The molecular representation of voltage-gated potassium ion channel [65] is shown in Figure 1.1 and 1.2. The main characteristic feature of any ion chan- nel is a pore through which ions can move [87]. The channel also spans the whole membrane and very often has an intra- or extracellular part which can play an important role in the regulation of the function of the channel. 1.2 Molecular Structures As was described in the previous section, ion channels are a class of pore-forming transmembrane proteins. The unique location of the channels partially in the mem- brane and partially in the solution causes extreme problems for current structure determination. Crystals required for X-ray characterization of ion channels are no- toriously dicult to obtain in the stable form due to problems with crystallization. For many years the only available structure was the structure of gramicidin A peptide which is used as a topical antibiotic for gram-positive bacteria. The rst resolved X-ray of the real biological structure was that of the KcsA channel [25]. The structure of this channel made possible signicant progress in understanding the selectivity of ion channels [16, 15]. Later, several other structures were resolved by X-ray methodology [47, 118, 46]. Although substantial eort has been made to determine structures by X-ray, there are still only 160 dierent solved structures of membrane-bound proteins. 2 Figure 1.1: Molecular representation of Kv1.2 channel (side view). 3 Figure 1.2: Molecular representation of Kv1.2 channel (top view). 4 1.3 Simulation of Ion Channels The atomistically based simulation of ion channels was restricted for many years because of the absence of reliable structures. The only available structure was gramicidin A , which has many features of biological channels and about which vast amount of experimental information has been written. The notable achievement was the correct barrier estimation [6, 4]. With the availability of KcsA and later aquaporin structures, the eld has our- ished [114, 42, 95, 85, 5, 58, 23, 99, 4, 86]. As timescales in the ion movement through pores is usually in the micro- or millisecond range, brute force molecular dynamics (MD) does not suce to answer many questions. For this reason, many multiscale approaches were developed to tackle the problem. Most notably, a care- ful prole estimation with Langevin Dynamics simulation gave arguably the rst molecular simulation of the ion current and selectivity [16] of the KcsA channel. Although the problem of selectivity is well understood, the gating of the voltage- gated ion channels is still a source of debate[115, 70, 7, 37, 60, 102, 48, 90, 28, 24, 79]. 5 Chapter 2. Multi-Scale Simulations 2.1 General Consideration Computer modeling of the function of macromolecules [49, 114, 9] and related sys- tems presents a problem of enormous complexity. Although available computer power has increased rapidly, there are, nevertheless, still many cases for which the use of brute-force simulations is clearly not the best approach. Furthermore, there exist many systems whose nature was correctly elucidated even before the emer- gence of the current level of computing power.There exist many examples of cases in which one cannot (and perhaps should not) progress without the use of a simplied model, and, in fact, recent years have witnessed a growing appreciation of the fact that the simulation of complex systems, and, in particular, the modeling of biologi- cal function, can greatly benet from physically sound simplications [31, 68]. The idea of such coarse graining rst appeared in protein folding studies and has since become a common, well-accepted, and powerful strategy. We try to emphasize that often focusing on minute details is not the best way to model a complex system. We also emphasize the importance of capturing the relevant physical features. In general, one may start from Einsteins advice to "make everything as simple as pos- sible, but not simpler." Here, one should of course dene what the question is, and what level of understanding is desired. It is also important to keep in mind the available computer power and its ability to give a stable (converging) result at that given moment in the history of the eld. Sometimes a dangerous assumption is that the availability and capability of running long trajectories (such as a millisecond 6 trajectory for an ion channel) represent a major breakthrough that will yield great insight into a system. This can be quite problematic, for example, because under- standing selectivity requires running multiple simulations and exploring the eect of dierent parameters on the overall current (e.g., [15]). The same holds true for a single long trajectory that explores protein folding [33]. Of course, with increas- ing computational power, it will be possible to run multiple long simulations using all-atom models. However, by that time, many of the key problems will probably have already been resolved by the use of simpler models. 2.2 Coarse-Graining and Dynamics Here we want to present general idea behind our Coarse-grain and complementing renormalization approaches, detailed description may be found in the next chap- ters. One of the most interesting problems in CG modeling is the challenge of simulating long-timescale events. The most obvious approach is the use of a fric- tion term to re ect the eect of the implicit thermal bath, using the Einstein or Wang & Uhlenbeck formulations. In fact, the use of frictional models to describe biological and chemical problems has a long history (e.g., [1, 116]). However, when focusing on obtaining similar physics in the full and CG systems, and on studying long-timescale biological processes, we cannot safely use the frictions obtained by standard frictional models. Our main approach for studying long-timescale processes involves taking an ex- plicit all-atom potential model and the corresponding simplied model (with its 7 simplied free-energy landscape) and running MD simulations and Langevin dy- namics (LD) on each model, respectively. This is done while imposing a series of constraints on both models, which force these systems to move along a given reac- tion coordinate on dierent timescales. The optimal friction is obtained by using the same constraints for the simplied and explicit models and adjusting the fric- tion in the simplied model until the timescale for the motion (for each constraint) becomes equivalent in the two models. However, the main point of the renormal- ization process is that we can force the explicit model to undergo large structural changes within reasonable computational time by using a large constraint, but it is essential to also interpolate the results to cases without constraints (which would require enormous computational time). In this respect, considering the diculties with attempts to estimate eective friction using a standard treatment, we insist on the idea that there must be some correspondence between the full and CG models and that the best way to obtain this correspondence is by simultaneously adjusting the friction and PMF of the CG model until the best agreement between the two models is reached. The power of the renormalization approach has also been demonstrated when exploring the selectivity of the KcsA ion channel, It is also useful to mention that the long-timescale behavior of the conformational coordinates has been eectively explored by the use of CG models. 8 Chapter 3. Long-time Simulations 3.1 Dynamics for Microsecond Simulations One of the long-standing problems in biomolecular simulations is getting meaningful time-dependent information on long-time functional processes. That is, whereas meaningful simulations of the function of biological systems on the picoseconds time range have been feasible for a long time,[109, 111] simulating processes in the millisecond and second time range has been a remarkable challenge. Of course, using free energy perturbation (FEP) and umbrella sampling (US) approaches (e.g., refs [110] ) has been extremely useful for evaluating activation barriers and then using transition state theory, but exploring true dynamical behavior on a long time scale, however, has remained a major challenge. Another long-standing challenge is involved in the eort to obtain converging free energy proles by FEP and US approaches within a reasonable computer time (e.g., see our discussion in ref [51]). For example, the insightful work of Jarzyn- ski [45, 44] that has been examined in exploring alternative to FEP approaches (using many short time simulations) is considered by many to provide somehow a way to circumvent the fundamental convergence problems. However, although the Jarzynskis identity is formally correct, it does not mean that the corresponding expression converges faster than standard FEP approaches. Furthermore, although related strategies [76] provided impressive analysis and formal prescriptions for using nonequilibrium simulations, we are not aware of quantitative demonstrations that such approaches provide major convergence acceleration. We are also not aware of 9 clear analysis that would provide a guide for choosing an eective pulling force. Sim- ilarly, the promising replica exchange approach[41, 97] has not been demonstrated to be much more eective than using repeated FEP calculations with random initial conditions. As far as the estimate of free energy proles is concerned, it is tempting to ex- amine so-called steered molecular dynamics (SMD) approaches [43, 76]. At present, the problem with SMD is that it is not clear how to convert the results obtained with large external force to the corresponding results that would be obtained with- out external force at equilibrium (see below). Therefore, whereas such an approach may tell us about the path between two known structures, it is hard to determine the eect of the bias on the path. The use of SMD for estimating activation barriers is even more problematic. Before describing our strategy for obtaining long time properties and for esti- mating free energy barrier, it is useful to clarify a general concern about current approaches that use SMD or related nonequilibrium approaches. Our point is best illustrated with the help of Figure 3.1, which presents a typical case where we like to know the PMF for moving from the initial point to the nal point and to do so in an ecient way. Now, a formal examination of this problem, from the perspective of nonequilibrium approaches (e.g., approaches that are related to the Jarzynski equality), is that we can get G(t) for the given move by considering the corre- sponding nonequilibrium work. Unfortunately, if we use a strong force, then we will go directly through path 3, which will very rarely or never sample points along the 10 least energy path (path 1) and thus will never give the correct PMF, unless we have an innitely long run. This indicates that we will have to follow the same strategy as in regular PMF calculations and to spend the same amount of computer time. Of course, one can argue that this is well known. Here we develop an approach [27] that tries to get the best correspondence be- tween the free energy and dynamics of a full explicit model [108, 112, 61, 50, 106, 78, 80, 63] (all-atom molecular dynamics, MD) and an implicit Langevin Dynamics [15, 49, 62, 91, 75, 92] (LD) model of reduced dimensions. The correspondence is obtained by using constraints of dierent strengths in both the implicit and the ex- plicit models and then adjusting the friction in the LD treatment to maximize the agreement between the time dependent responses of both models. Using external constraint allows us to run all-atom molecular dynamics simulations in a reasonable computer time (from several hours to several days using serial code) and then to obtain key information for LD simulations of long-time processes. The renormaliza- tion approach had already been used very eectively in studies of the selectivity of ion channels [15] in simulating proton transport, and in exploring protein dynamics. In fact, this approach has been validated on some level in ref [16]. However, up to now, we have not provided a systematic validation and a clear rationale for our strategy, and this is done in the present work. 11 Figure 3.1: Demonstrating the nature of response to a 1-D force in a 2-D system. The Figure is aimed at clarifying the fact that nonequilibrium type approaches can never provide reliable PMF with short simulations and large forces. In such case, the system will pass through path 3 and will not collect information about the actual least energy path. Of course running an innite number of short trajectories with strong force will also rarely sample path 1, but this would not oer an ecient sampling strategy. 12 3.2 Gramicidin A: Method Our aim is to nd a practical way to perform long-time simulations (in the micro- or milliseconds time range) of complex systems while still capturing the correct physics. We look for an approach that still considers the explicit response ofhQi to dierent forces (dierent K), but then generates the corresponding friction and eective potential in an LD treatment, where Q can be a set of generalized coordinates (in the case of the system described in the present work, it is the z component of the ion coordinate), and K is a force constant in the present case. In other words, following our renormalization approach[16, 15, 81] we can write and represent the full system by a simplied system that is propagated by an LD simulation using m q =m _ q @G @q @H 0 @q +A 0 (t) (3.1) where runs over the x, y, and z coordinates of the ion, m and are the mass and friction of the ion, and G is the eective free energy function that includes mainly the electrostatic prole but also the steric restriction term used in ref 15. A(t) is a random force that satises the uctuation-dissipation theorem and H 0 is a perturbation, which in general can be written as: H 0 =KQ (3.2) .Integrating [56] the LD equation, we can obtainhQ 0 (t)i K and compare it with the corresponding coordinate in the full system (with the corresponding constraint). 13 Such a procedure allows us to nd the best friction constant (or if needed to use an even more sophisticated memory kernel) and use it in long-time LD simula- tions. This approach is described in more details below for the case of gramicidin A channel. Our challenge is to simulate the dynamics in the space dened by the eective coordinates in a way that it would represent the long-time behavior of the real system. This is done by adopting a hierarchical renormalization model, where we force a model of a limited dimension with a coordinate setq (which will be described here by three coordinates) to represent an explicit all-atom model with a coordinate set r. We start by mapping the full model along the reaction coordinate and generating a 1-D eective free energy surface. In the next step we force the two models to have equivalent dynamics. This is done by starting with a regular MD simulations of the time dependence of moving between two sets of coordinates (z 1 and z 2 ) of the explicit model under the in uence of a constraint (or pulling) potential U con =K(zz 0 ) 2 (3.3) where K is a force constant and where a large value of K is needed to obtain this transition in a reasonable time. The value of z0 (z0 = 20000) was chosen to ensure that the force (jFj 2Kz 0 ) on the ion is (almost) constant during the pulling process, where the constraint potential aects directly only the ion in the simulation system. Next, we move to the implicit model and evaluate the time 14 dependence of moving between two coordinate sets q 1 and q 2 (that correspond to z 1 and z 2 ) using a LD treatment as outlined in eq 3.1 m q z =m _ q z @G @q z 2Kz 0 +A 0 (t) (3.4a) m q =m _ q @G steric @q +A 0 (t) (3.4b) where qz is the z coordinate of the ion, runs over the x and y coordinates of the ion, and G steric is the change in free energy associated with the shift from the channel z axis (see ref [27] for details). The key issue here is to force the dynamics of the implicit (simplied) model to correspond to that of the explicit (full atomistic) model. Therefore, we force the time dependence for moving between the two generated sets of coordinates in each model (e.g., the time required to move between two signicantly dierent coordi- nate sets in the explicit model should correspond to the time of moving between the equivalent coordinates in implicit model) for dierent values of the applied con- straint. We then rene the friction, , by comparing the passage time and overall time dependence obtained from the simulations in the dierent models. It is also possible to validate and ne-tune the dynamical properties of our renormalization approach by requiring that the velocity autocorrelation of z(t) or other related properties have similar behavior in the two models. Once the model is calibrated by the above procedure, we can use it to explore the behavior of the system on a very long-time scale by simply running the calibrated implicit model for a long time 15 Figure 3.2: Benchmark system of a truncated gR with a sodium ion. and exploring long time properties, as was done, for example, in our study of the coupling between conformational and chemical dynamics [81]. All of the simulations with the explicit model were performed with the MO- LARIS program package. The parameters used in the present work are exactly the same as those in ref [51]. The LD simulations were done with the CHANELIX [16] module of MOLARIS. We use the truncated gramicidin A (gR) in the gas-phase model depicted in Figure 3.2. The initial structure (pdb 1JNO) was mutated to all glycine residues, and position constraint of 10kcal=Angstrom 2 on allC atoms was used to create a stable benchmark (this was done because the membrane is removed and the simplied model channel cannot keep a reasonable shape by itself). The system 16 used in the current work as well as the direction of the pulling force can also be seen in Figure 3.2. 3.3 Gramicidin A: Results The initial step in the validation of the renormalization approach requires one to obtain the potential of mean force (PMF) for the given simulation system. In the present case, we already have a reliable PMF for the simulation system because the specic PMF had already been obtained in ref [51]. We also considered changes in the PMF, and these were obtained by the same procedure used in ref [51]. At any rate, the PMFs obtained for the explicit model (Figure 3.3) was used as the eective free energy prole in the LD model. We also generated a model that allowed us to modify the PMF and to be able to see how the ion moves through the channel without a pulling potential. To achieve this goal, we altered our simulation system by adding four xed charges outside the ion channel with a distance of 8 from the center of the channel (Figure 3.4). Changing the charge on the atoms allowed us to change the activation barrier in a parametric way without modifying the friction. Because the low barrier system has a slightly dierent minimum along the PMF prole, we dene the average time required for the ion to pass the channel as the time for moving the ion from z = -6 to 6 . Next, we started to generate thehz(t)i as well as the average time required for the ion to pass through the channel .This is dened as the time for moving ion form z=-8 to 8 through the barrier of the explicit model. The calculations 17 Figure 3.3: PMFs for the high barrier model (from ref 5) and the low barrier (9.3 kcal/mol) model of Figure 3.4 with charges q = -0.3. The PMF of the system with the low barrier is shifted to have the same minimum as in the system with the high barrier. 18 Figure 3.4: Truncated gR with four additional charges at the center of the system. The change of the charges provides a systematical way to alter the original PMF. were done by using the gR models with dierent K values. Typical trajectories of the explicit model for dierent force constants are depicted in Figure 3.5, whereas Figure 3.6 compares trajectories and autocorrelation functions for the explicit and implicit models. As shown in Figure 3.6, we could easily force the LD model (by adjusting the friction coecient) to have similar trajectories and similar velocity autocorrelation function to that of the explicit model. Obviously, using the same friction for all force constants is not fully justied, but this seems to be reasonable for the models studied to date. At any rate, the results obtained for our systems are summarized in Figure 3.7, where we provide the overall dependence of the rst passage time on the force constant of the explicit and implicit models. In this case, we determine the average time by using 16 or 32 19 Figure 3.5: Two trajectories obtained with the explicit model for dierent pulling forces, K. The simulations were done for the high barrier case, and K is given in kcal/mol/Angstrom 2 20 Figure 3.6: Typical MD and LD ion trajectories for the explicit and implicit models, respectively (A), as well as the velocity autocorrelation functions (B) for the high barrier model and for K = 0.0001 kcal/ mol/Angstrom 2 21 Figure 3.7: First passage time for explicit and implicit models for dierent force constants. The data shown are for two models: A high barrier model (from ref [51]) and a low barrier model (truncated gR with four negative charges). trajectories for each force constant. Typically, the total time is 5-10 times larger than the average time required for the ion to pass the channel (the rst passage time) for a given force constant. As seen from the Figure, we obtained an excellent agreement between the behaviors of the implicit and explicit models over a very large range of the rst passage time. 22 After establishing our ability to obtain the same rst passage times in the two models, we also evaluated a very long time trajectory in the absence of any force (but with a constraint that prevents the ion from going out of the channel). The simulation [27] of the corresponding test system with a low barrier and without the pulling potential took several CPU days and yields reasonable agreement between the explicit and implicit models (see Figure 3.8). Remarkably, the long time results described in Figure 3.7 were obtained without any readjustment of the friction coef- cient (using the value of 8 ps 1 , obtained with the large barrier and strong pulling force). Considering the above results, we believe that we establish that a renormal- ized ion channel system can be used in reliable LD simulations to investigate ion current and selectivity. 3.4 Barrier Estimation The renormalization approach not only allow us to simulate long-time dynamics but also provides a powerful way for evaluating PMFs [51, 27, 76, 84, 57, 104, 71]. This non-equilibrium strategy is also developed and demonstrated here. We examine here options for accelerating calculations of activation free energy. That is, after calibrating the friction term in the simplied LD equation, we can consider cases when the potential of the explicit model is changed (lets say because of a mutation) and then try to get the same passage time andhz(t)i for the new explicit modied surface and the LD eective potential. This is done while keeping the friction unchanged and adjusting the eective potential until we reproduce the same rst 23 Figure 3.8: Long time trajectories of the explicit and implicit models for the sys- tem with a low barrier (9.3 kcal/mol) and without a pulling force. The Figure demonstrates that we can reproduce the long time behavior of the explicit model. 24 Figure 3.9: Original PMF for the explicit model and the modied PMFs for the system with four additional positive charges. g 6= is dened as a dierence in the barrier height of the system with four xed charges and the high barrier system (without four xed charges). passage time andhz(t)i as in the explicit system. Once this is accomplished, we expect the eective potential of the simplied model to give a reliable estimate of the PMF of the explicit system. To do so, we needed a simple and controllable way of changing the activation barrier in the explicit model, without drastic change of the factors that control the friction. This was accomplished (as discussed above) of ref [51] by using the system of Figure 3.4. Here we considered several cases with known change in the PMF 25 (shown in Figure 3.9) and evaluated the eect of strong pulling force on the passage time of these test systems, pretending that we do not know the activation barriers. We then forced the LD simulations (with a friction coecient that corresponds to the system without four xed charges) to reproduce the rst passage time of the explicit system by adjusting the eective potential for the LD runs. The adjustment was done by adding a Gaussian shape potential with width of 2 (Figure 3.9) to the initial known potential. As shown in Figures 3.9-3.11, the potentials that allowed us to obtain the best match between the implicit and explicit passage times reproduce the actual activation barriers in the explicit system. The next challenge has to determine the complete PMF when the shape is not known. In this case, we had to determine the PMF from several segments and such an approach was explored using the system with the high barrier (see Figure 3.12). The actual simulations were done on two segments and then extended by exploiting the symmetry of the system. The determination of the PMF from separate segments required the use of a dierent force constant to achieve a passage time that is not too fast and not too long (a rst passage time in the range of hundreds of picoseconds was used). It should also be noted that our approach requires one to obtain the PMF of the segments by running uphill trajectories because downhill trajectories with and without external force have a similar rst passage time. However, this limitation can be easily overcome by applying pulling potential in the opposite direction. Next, we explored dierent strategies for estimating the barrier height for the explicit system. The starting point in these studies was the crucial need of having 26 Figure 3.10: Relative rst passage time as a function of the g 6= (dened in Figure 3.9) of the explicit model for three dierent pulling forces. The estimated g 6= is independent of the magnitude of the pulling force. The relative time is dened as the ratio of the rst passage time (with a given pulling force) for the system with a modied PMF to the rst passage time for the original high barrier model. 27 Figure 3.11: Comparison of the g 6= obtained by the renormalization approach and the exact results obtained with the careful analysis of the explicit model. 28 reasonable eective friction. In principle, we can obtain a reasonable approximation from a renormalization study of regions with low activation barriers, but here we used the already known optimal friction and tried rst to estimate the barrier height by using eective potentials with a Gaussian shape. The corresponding results are presented in Figure 3.12, and as seen from the Figure, these results are only qualitatively correct. We obtain a deviation of 5 kcal/mol from the correct barrier of the explicit model. Therefore, we adopted a second approach where we know the shape of the potential (this can be done with the approach considered below). In this strategy, we generated potentials that have the same minima and barrier as the explicit potential and systematically scaled the barriers by a factor (in the range of 0.6 to 1.6). Now we obtained an excellent agreement between the renormalized PMF and the corresponding explicit result (see Figure 3.12). The above validation emphasized the advantage of knowing the general shape of the potential in conducting a renormalization process. This highlighted the advan- tage of simple models in obtaining the general shape of the potential or of using the above renormalization approach in studies of mutational eects when the changes in the potential are small. Of course, we also have to have a reasonable idea of the best eective friction, and this is another crucial advantage of the renormalization idea. In this respect, we would like to note that our approach for obtaining the PMF depends on having the optimal friction. However, as mentioned above, the 29 Figure 3.12: Comparison of the PMF of the explicit model (black) and the PMF (red) obtained by using the renormalization approach (using LD simulations that reproduce the rst passage time for dierent segments of the explicit model) as well as estimation of the barrier height. 30 optimal friction can usually be estimated without the full renormalization treat- ment, by performing renormalization studies of regions with low activation barriers of the given system. At any rate, the advantage of the present approach is most apparent when we can use a relatively strong force constant and short simulation time. This can be extremely useful in assessing the rate of biological processes with large barriers. Of course, if we like to evaluate the exact details of the PMF, in particular, in segments where the PMF is shallow and where the reaction coordinate does not follow a simple path, then we need to use smaller constraint forces and longer simulation times. At this point, it is useful to discuss the simulation times of the dierent models. Obviously, the simulation time depends on the force constant. The longest simula- tions reported for the model with the high barrier (from ref [51]) were around 0.2s, which corresponds to several days of computer time using one core of Dual Quad- core AMD Opteron 2.3 GHz. The corresponding calculations for the implicit model took several hours. The PMF simulations for the model with the high barrier (from ref [51]) took, with averaging for backward and forward runs, about 16 h, whereas evaluating the PMF with the renormalization approach requires 2 h when a large force constant was used for determination of the barrier height. We expect greater speedup in solution, and this issue is under investigation. Of course, the main time saving in our model is not in the PMF calculations but in the ability to run long time simulations on the implicit model. Here we can obtain a microsecond (0.1 s) simulation in Figure 3.8 with the implicit model (in several hours), whereas we will 31 need for this several days with the explicit model (even for our small simulation system). For larger systems, the saving is much more remarkable. 3.5 Discussion The renormalization method is developed and validated on a system that re ects the main complexity of realistic ion channels. The renormalization approach is found to be reliable and arguably presents what is perhaps the rst approach that allows one to exploit SMD in quantitative studies. The renormalization approach has already been used in studies of several impor- tant problems [81, 15], but the present work provided the rst systematic validation of this method. It is established that we can reproduce the long time behavior of large complex systems by using LD simulations of a renormalized implicit model. This is done without spending the enormous time needed to obtain such trajectories in the explicit system. The present study also provided a promising breakthrough in accelerated evaluation of general shape of PMFs. This is done by adjusting the eective potential in the implicit model to reproduce the same passage time as that obtained in the explicit model under the in uence of an external force that provides a way to extract the PMF without inverting the time need from regular PMF calcu- lations. This approach, which is illustrated here in realistic calculations, is expected to provide a major help in the study of complex landscapes. It also should provide a very useful way for assessing activation free energies of macromolecular processes. The main point in our approach is the idea that there must be some correspondence 32 between the full and the implicit models and the best way to get this correspon- dence is to adjust the friction and PMF of the implicit model simultaneously until the best agreement between the two models is obtained. The approach that is probably the closest in spirit to our approach is the so- called surrogate process approximation (SPA) of Calderon and coworkers [18, 17]. This approach that appeared recently independently from our previous renormal- ization approach uses SMD to construct PMFs. Calderons approach was actually validated in calculating the PMF of the gR channel and obtained encouraging re- sults . At any rate, the SPA approach is focused on evaluating PMF, whereas our approach is more directed toward providing the long time properties of the system, which cannot be estimated by other approaches. In discussing alternative strategies, it is important to state that several ap- proaches outlined eective ways for evaluating PMFs with a known friction con- stant. The problem is that the proper friction constant is not known. Here our idea of obtaining both the friction and PMF from SMD-type treatments presents a crucial advantage. In this respect, we like to reemphasize that the friction obtained from standard treatments (namely, the Einsteins or related formulation) may not provide the optimal eective friction. It must be claried here that although the nature of frictional constant has been understood by giants like Chandrasekhar, [21] Kubo,[56] and others, the actual value depends on the model used (e.g., ref [105]). Furthermore, the friction in the eective LD model may re ect factors that are only considered implicitly. For 33 example, if the PMF missed some corrugation, then it will be re ected in the change of the optimal friction. The gR system chosen might be considered by some to be a simple 1-D system. However, in fact, we have here a system with the main complexity of ion chan- nels. That is, the reorganization of the channel dipoles upon transfer of the cation provides a complex solvent coordinate with the typical solvent relaxation dynam- ics. Capturing the coupling between the solvent and solute motions is a challenge that has to be addressed by the renormalization approach. Here the friction and the PMF must re ect the behavior of the correct solute/solvent system with its eective 2-D features. The demonstration of the power of the renormalization approach is likely to open new applications. For example, this very promising approach is likely to advance in simulating the selectivity and activation of voltage-activated ion channel in a way that is still reliable and yet feasible with the available computer time. 34 Chapter 4. Conformational Changes Considering the eectiveness of the renormalization approach in the case of ion chan- nels, we are also interested in the performance of this approach in other types of problems, such as conformational transitions in proteins. In such cases, we are deal- ing with more complex landscape, and the renormalization procedure can be more demanding. As an example, we consider the conformational change of adenylate ki- nase (ADK), which has been used as a generic model in our study of conformational changes[81]. In this case, we are dealing with a large conformational change that has been explored by several workers (e.g., ref [66] ). In our preliminary study, we considered the exact details of the conformational transition to be highly irrelevant because our study was about the general issue and not about any specic enzyme. Here we considered the ADK type renormalization problem in more detail, exam- ining three models: the full explicit model, the CG model where each side chains is replaced by single interaction centers [73], and a 2-D model with one conformational coordinate and one chemical coordinate. We started the renormalization by apply- ing the same forces to the three models while adjusting the corresponding frictions in the CG and 2-D models. This was done, however, by using a simple single barrier potential in the 2-D model rather than by evaluating rst the PMF in the full or CG model . In the current case, we applied weaker forces than those used in ref [81], and the corresponding results are presented in Figure 4.1. As seen from this Figure, we can nd friction constants that satisfy the renormalization requirement; however, we did not try to establish the uniqueness of our results. Furthermore, 35 calculation of the autocorrelation of the conformational coordinate shows that we can reproduce the autocorrelation of the explicit model with very dierent frictions (namely, 0.5 and 100 ps-1). In fact, attempts to simulate small peptides by LD approaches (e.g., ref [32]) seem to give similar results with frictions ranging from 0.5 to 50 ps 1 . In our opinion, the renormalization approach should provide unique results or at least a consistent way for modeling the long time motions and free energy proles in studies of protein conformational changes. However, this should be demonstrated in studies of relevant systems. The most obvious rst step would be to estimate the actual PMF of the protein conformational change and to rene the renormalization results. Such a study can start by considering the correspondence between the CG and the 2-D models rather than starting from the full model (the CG model, in fact, may provide a more reliable PMF than the full model). Another useful strategy can be provided by using relatively weak forces that would allow one to pull the system along one small segment at a time. This approach can also be used to construct the PMF simultaneously. The performance of the renormalization approach in such cases would be validated in a subsequent study. However, we also like to emphasize that we do not see any fundamental problem in using the renormalization strategies in studies of protein conformational changes, except the likely requirement to perform the renormalization on segments of the overall reaction coordinate change to obtain more reliable results. 36 Figure 4.1: Qualitative renormalization study of the conformational change in ADK for the full-atom, coarse-Grain, and 2-D LD models. It describes the timed appen- dance of the response to force of the three models. 37 Chapter 5. Voltage-activated Ion Channel 5.1 Motivation and System Understanding the detailed mechanism of the activation of voltage- gated ion chan- nels is a problem of great current interest. Reliable molecular simulations of voltage eects present a major challenge because meaningful converging microscopic simula- tions are not yet available and macroscopic treatments involve major uncertainties regarding the dielectric constant used and other key features. The current work reports methods to overcome some of the above challenges by using coarse-grained (CG) model in simulating the activation of the Kv1.2 channel. The CG model has allows to explore problems that cannot be addressed at present by fully microscopic simulations, while providing insights on some features that are not usually consid- ered in continuum models, including the distribution of the electrolytes between the membrane and the electrodes during the activation process and thus the nature of the gating current. Furthermore, the clear connection to microscopic descriptions combined with the power of CG modeling oers a powerful tool for exploring the energy balance between the protein conformational energy and the interaction with the external potential in voltage-activated channels. The elucidation of the structure of voltage-activated ion channels and biophys- ical studies (e.g., refs. [39, 25, 64, 19, 20]) have provided key information about the relationship between the membrane voltage and the gating process. However, despite these great advances, we still do not have a clear picture of the correspond- ing structurefunction correlation. Furthermore, although there has been signicant 38 progress in the computational modeling of the energetics of ion channels (e.g., refs. [16, 100, 67, 12, 36]), the understanding of the voltage activation process is rather limited. Not only have the exact structural changes not been fully determined, but the energetics of the conformational transition and the coupling to the external voltage are far from being understood. To clarify the current challenges, it is useful to follow the general concept de- picted in the schematic diagram of Fig. 5.1(A). As shown in the gure, in the resting state (with a negative potential), the channel is in its closed state, where the free en- ergy of moving to the open state is positive. The application of a positive potential stabilizes the open state, leading to a conformational transition to this state, where the channel allows ion transport. The overall process can be captured once we have the relevant free energy landscape, which should follow the trend of Fig. 5.1(B). This landscape depicts a path from a region of closed channel at negative potential to a region of open channel at positive potential. The fundamental challenge is to capture the energetics of the conformational transition and to reproduce the eect of the external potential that would overcompensate the conformational energy and move the system to the open state. Further challenge is associated with capturing the microscopic nature of the gating charge and the barrier for the conformational transition, as well as the uctuations during the gating process. One of the major diculties in simulating voltage-activated ion channels is the treatment of the membrane potential and its eect on the protein congurations and the electrolyte distribution. These issues have been explored by interesting 39 Figure 5.1: Illustrating the free energy balance in a voltage-activated ion channels. (A) The gure considers schematically the activation process, where the protein moves from the closed to open state due to the allosteric eect of the external potential (V). When the external potential is negative (V < 0), the closed con- formation is more stable and the channel is blocked. Application of a positive potential (V > 0) stabilizes the open conformation and the system moves to its ion conduction state. (B) A tentative landscape for activation of the Kv1.2 system in the two-dimensional space dened by the protein closed to open conformational coordinate and the external potential. Note that the potential coordinate is a co- ordinate of an external constraint and as such the motion from negative to open potential does not have to be downhill (because it does not include the external power source). Note also that the surface was generated with empirical valence bond-type parabolas and from the data of free energy table. Thus the increase in energy, upon moving from the minima toward the boundary of the map, does not re ect the real surface. 40 PoissonBoltzmann (PB) studies ([36, 83]), but the results re ect a macroscopic perspective and problematic dielectric constants (for a review, see ref. [114]). In principle, one may try to use molecular dynamics simulations with all-atom explicit models ([53]), but at present we believe that such an approach is not likely to provide converging free energies (due to both the challenge of obtaining stable solvation free energies in protein interiors and the diculties in capturing the response of the ionic atmosphere). Probably the most eective current strategy should involve the use of coarse-grained (CG) models with proper electrostatics [26]. In view of the above considerations, we have recently developed a unique strategy based on our existing CG model, while adding the crucial eect of the ions in solution and the membrane potential. This model enables us to simulate membrane proteins in the presence of external potential and electrolyte solutions. In the present work, we focus on modeling the function of the Kv1.2 voltage-activated channel. The simulations provide interesting insights about the energy balance, gating charges, and gating uctuations in voltage-activated channels. Here we brie y discuss our CG model and validation, detailed description of the methodology is given is subsequent sections. Our CG treatment is aimed at modeling the protein/membrane system and its interaction with the external potential. The modeling of the protein/ membrane system has been described in detail also in [88]. Here we outline some of key points about the description of the electrolyte solution and the membrane potential. Our simulation system for studies of membrane potential is described in Figure 41 5.3. The system includes a simulation box that explicitly includes the membrane containing the protein (region I), an optional region with explicit ions (region II, which is not considered in the present study), and a grid representing the electrolyte solution (region III). We also add a "bulk region" far away from both the membrane and electrode surfaces, as a specialized way for spanning the space between the membranes to the electrodes, without using an enormous grid. The explicit grid model re ects a compromise between the fully microscopic models of the electrolyte (in an implicit solvent) to a macroscopic model The main features of our approach have been validated recently by reproduc- ing the trends in DebyeHuckel and GouyChapman models, as well as the trend in the more challenging problem of two electrodes and electrolyte solution . We have also simulated the expected capacitance of a neutral membrane in an electrode- electrolyte system and demonstrated that the ions in the solution provide almost a complete screening of the external potential, which starts to increase only when it reaches the membrane boundaries. This behavior is similar to what is expected from macroscopic considerations. 5.2 Coarse-Grained Model for the Channel The present work uses a CG model that describes the main-chain atoms by an explicit model and represents the side chains by a simplied united atom model, as shown on gure 5.2. This CG model provides a more advanced treatment of the electrostatic eects than most current CG models. More specically, the present 42 Figure 5.2: From full atom representation to Coarse Grained representation. model, which is a modied version of our recent works [49, 73, 88, 82, 26], expresses the overall free energy (in kcal/mol) as G total = G main + G mainside + G side (5.5) The main-chain atoms are treated explicitly with implicit solvent corrections, and the main-side interaction involves van der Waals and screened electrostatic 43 terms [73]. The major and most relevant part of our treatment comes from the G side , which is given by G side = G vdw side + G elec side + G hyd side (5.6) where the rst term describes the eective van der Waals interactions between simplied side chains, which is described in detail in ref. [73]. The second term represents the electrostatic interaction between the ionizable residues (see below) and the last term represents the hydrophobic contributions which are not included implicitly in the rst term (see below). The G elec side side term is given by the following equation (see also below for recent modication): G elec side =2:3 RT X i Q i (pK w a;i pH) + G QQ + G self (5.7) wherei runs over the proteins' ionized residues,pK w a;i is thepK a of thei th residue in water andQ i is the charge of thei th residue in the given ionization state. G QQ is the charge-charge interaction free energy, which is given (in kcal/mol) by: G QQ = 332 X i<j Q i Q j r ij eff (5.8) where the distances and charges are expressed in A and electronic charge units, and eff is the eective dielectric for charge-charge interaction, which re ects the idea established in many of our earlier works (e.g. [114, 82]) that the optimal value 44 is large even in protein interiors (namely eff > 20). This type of dielectric has been found to provide very powerful insight in recent studies of protein stability (see [82, 107]). The ionization state of the protein residues were determined by the Metropolis Monte Carlo approach of [73] for the given pH. A key element in our approach is the treatment of the self energy, G self , associated with charging each ionized group in its specic environment. This term is given by: G self = X i U np self (N np i ) +U p self (N p i ) +U mem self (N mem i ) (5.9) where U designates eective potential, i runs over all ionized residues, U np self , U p self and U mem self are the contributions to the self-energy from non-polar (np) residues, polar (p) residues and membrane (mem) atoms (more precisely, membrane grid points as claried below), respectively. Here N np i , N p i and N mem i are, respectively, the number of non-polar residues, polar residues and membrane atoms in the neigh- borhood of thei th residue. Note that the non-polar contribution for the membrane is taken into account separately in the hydrophobic term (described below). The empirical functions U p self and U np self are given by: U np self (N np i ) = 8 > > > < > > > : B self np exp[0:05(N np i 15) 2 ]; 0<N np i 15 B self np ; N np i > 15 (5.10) 45 and U p self (N p i = 8 > > > < > > > : B self p exp[0:05(N p i 10) 2 ]; 0<N p i 10 B self p ; N p i > 10 (5.11) The number of non-polar residues neighboring the i th ionized residue is expressed by the analytical function: N np i = X j(np) F (r ij ) (5.12) with F (r ij = 8 > > > < > > > : 1 r ij r np exp[ np (r ij r np ) 2 ] r ij >r np (5.13) Similar equations were used for the number of polar residues neighboring the i th ionized residue (with the parameters r p and p ) and for number of membrane grid points neighboring thei th ionized residue (with the parameters mem and mem ). The relevant parameters are given in tables below. Our model was also applied to membranes which are represented by grids of unied atoms (as we have done in our previous studies (e.g. see ref [88]) and this grid is used in evaluatingN mem by the equivalent of equation 5.12. The membrane grid is used in calculating U self;0 mem where we have U mem;0 self = 8 > > > < > > > : B self mem exp[0:05(N mem i 28) 2 ]; 0<N mem i 28 B self mem ; N mem i > 28 (5.14) 46 where U self;0 mem is given by the same expression as equation 5.12. Additionally, we have used U self mem to the form U mem self = 8 > > > < > > > : U self;0 mem exp[( R min 18 9 ) 2 ]; R min 18 U self;0 mem ; R min > 18 (5.15) where R min is the distance to the closest solvent molecule, which is determined by the grid around the system, and using the distance to the closest grid point. The membrane grid was treated with continuous derivatives in order to reduce the need for generating a new grid when the protein is displaced or changes its structure. This was done by building a continuous membrane (instead of deleting membrane points that appears in direct contact with the protein). Accounting for the fact that the membrane grid should be deleted upon contact with the protein atoms we replaced the standard van der Waals interaction between the protein and the membrane by U vdw = X i<j h A ij ( +r 6 ij ) 2 B ij +r 6 ij i (5.16) whereA ij andB ij are vdW parameters for the interaction of the i th andj th atoms, r ij is the distance between the two atoms, and is a vdW cuto parameter. The corresponding parameters are given in ref [73]. The eect of zwitterionic membrane head groups was simulated by placing pos- itive and negative charges on the outer and the subsequent layer of the membrane grid, respectively. The corresponding electrostatic interactions with the protein 47 charges were treated with eff = 20, which has been justied by earlier studies of the eld from membrane head groups [113]. The term G elec side of equation 5.7 has been recently rened following further studies of protein stability and is now given by: G elec side =1:38 X i hQ i i (pK i a pK w a ) (5.17) +1:38 X i h (Q uf i hQ i i) (pK w a pH) (1e jpK i a pK w a j ) i +V QQ where,hQ i i is the MC averaged charge of the i th residue and = 0:1. The parameters and the functions used in our CG treatment are being rened continuously and the most recent renements are outlined in the following sections. The above electrostatic treatment involves a self-consistent treatment of the interdependent self-energy, chargecharge interaction, and the external pH (where the ionization state is determined by a Monte Carlo treatment of the energetics of Eq. 5.5). The hydrophobic eect was treated originally by adopting the same model used in the self-energy calculations. More recently, we have rened this model by trying to get a more explicit description of the water exclusion eect and the dierence between small and large hydrophobic residues. The new model replaces the previous implicit estimate of the missing solvent molecules by a more explicit estimate. This estimate was obtained by generating a grid of a radius 7 and spacing of 3 around each hydrophobic side chain, excluding the sites within a certain van der Waals cuto distance from the protein or membrane atoms and then collecting 48 all the grid points in the surface between R in and R out = R in + 2, and obtaining the number of water molecules in the ring (or more precisely in the spherical shell). This number, N Ring i was then used in G hyd = X i [U hyd np (N Ring i )+U hyd p (N p i )+U hyd mem (N mem i )]+B bonds N bonds +C hyd (N hyd ) 2:5 (5.18) where i runs over all nonpolar (hydrophobic) residues. The function used for evalu- atingU hyd np (N Ring i ) np is given below, whereasU hyd p andU hyd mem mem were evaluated as before by the same treatment used in the self-energy calculations. The parameters for these terms are given in table below . TheN bonds term of the equation represents the compensation of the hydrophobic eect by entropic factors and other eects and is represented by the number of single bonds times a factor B bonds (B bonds = 0:15). N hyd represent the number of hydrophobic residues, and the corresponding term accounts for the nonadditivity of the hydrophobic contributions (Chyd = 0:0005). The function U hyd np is given by U hyd np (N Ring i ) =C np i (1exp[a u np (N water i Ring i )]=N water i ) (5.19) where C np i i values are given in table below, N water i is the maximum number of water molecules (this value depends only on the residue type), and N Ring i is the number of water molecules in the spherical shell with Rin < R < Rout. The relevant parameters are given in table below, a u np = 1:4. 49 In describing movement of the protein from water to membranes (and to some extent to the proteins), it is important to consider the change of solvation of the main chain (see, e.g., refs. [11] and [22]). In order to model this eect, we consider the solvation of each C atom using U i ;np = 8 > > > < > > > : 1 0N np 10:3 exp[1:4(N np 10:3) 2 ] N np > 10:3 (5.20) and U i ;mem = 8 > > > < > > > : 1 0N mem 8 exp[0:1(N mem 8) 2 ] N mem > 8 (5.21) where i runs over the C atoms. We further dened a mixing rule U i =U i ;np U i ;mem (5.22) With these relationships, we approximated the main-chain solvation contribu- tion by G main solv =2 X i U i (5.23) The current CG treatment also modied the expression used for the hydrogen bond- ing between the main-chain atoms using the following: 50 U HB = 8 > > > < > > > : A rr 0 Aexp[(rr 0 ) 2 ] r>r 0 (5.24) where the relevant parameters are given in the next sections. Furthermore, to re ect the fact that it should be easier to break the hydrogen bond (HB) in water (due to solvation) than in the membrane, we have scaled the above HB energy by considering its exposure to the solvent and the chance for hydrogen bonding to the solvent, using U cor HB =U HB [3(1 U i U j )]=4 (5.25) Further renements and calibrations of the model have been subject to the con- straint of keeping the calculated absolute folding energies close to the observed value. Here we provide the current set of our CG model. This set is constantly being rened. 51 Residue R 0 i ( A) 1 (kcal/mol) 2 (kcal/mol) Polarity A 4.26 0.03 0.05 nonpolar C 4.79 0.06 1.17 polar D 4.45 0.13 1.17 polar E 4.68 0.17 1.17 polar F 4.55 0.1 0.13 nonpolar H 4.75 0.2 1.17 polar I 4.68 0.13 0.18 nonpolar K 4.74 0.16 1.17 polar L 4.72 0.13 0.18 nonpolar M 4.76 0.12 0.17 nonpolar N 4.59 0.13 1.17 polar P 4.72 0.23 0.31 nonpolar Q 4.79 0.17 1.17 polar R 4.47 0.24 1.17 polar S 4.54 0.06 1.17 polar T 4.62 0.1 1.17 polar V 4.7 0.1 0.14 nonpolar W 4.48 0.25 0.35 nonpolar Y 4.51 0.27 1.17 polar Table 5.1: Parameters for the G vdw side term. In the table 1 is used for interac- tions involving residues of dierent polarity whereas 1 + 2 is used for interactions involving residues of the same polarity as well as interactions with ionized residues 52 region A i (kcal/mol) i 0 ! i ;0 i 0 ! i ;0 - helix -1 -95 8 -5 16 - sheet -10 -150 10 175 16 `-turn -3 -80 6 65 5 L- - helix -5 60 5 40 16 Table 5.2: Parameters for the U phipsi mm term. The angular values in the table are given in degree. B self;ion p Equation B self;ion np Equation B self;ion mem Equation ARG Ionizable -0.3 5.10 2.5 5.11 10 5.14 LYS Ionizable -0.3 5.10 2.5 5.11 3 5.14 GLU Ionizable -0.4 5.10 4 5.11 10 5.14 ASP Ionizable -0.4 5.10 4 5.11 10 5.14 HIS Ionizable -0.3 5.10 2.5 5.11 10 5.14 Table 5.3: Self Energy Parameters for the Ionizable Residues 53 B self;pol p Equation B self;pol np Equation B self;pol mem Equation THR Polar -0.37 5.10 0.37 5.11 0.13 5.14 SER Polar -0.69 5.10 0.69 5.11 0.23 5.14 CYS Polar -0.03 5.10 0.03 5.11 -0.01 5.14 ASN Polar -1.27 5.10 1.27 5.11 0.43 5.14 TYR Polar -0.03 5.10 0.03 5.11 -0.35 5.14 GLN Polar -1.15 5.10 1.15 5.11 0.38 5.14 Table 5.4: Self Energy Parameters for the Polar Residues. Equations identical to 5.10, 5.11 and 5.14 were used for polar residues (but the parameters used are dierent and listed in the table) B self;np p Equation B self;np np Equation MET Hydro 0.55 5.11 -0.34 5.14 ILE Hydro 0.55 5.11 -0.56 5.14 VAL Hydro 0.55 5.11 -0.23 5.14 LEU Hydro 0.55 5.11 -0.63 5.14 ALA Hydro 0.55 5.11 0.25 5.14 PRO Hydro 0.55 5.11 0.07 5.14 PHE Hydro 0.55 5.11 -0.86 5.14 TRP Hydro 0.55 5.11 -1.04 5.14 Table 5.5: Self Energy Parameters for the Hydrophobic Residues from polar and membrane neighbors. Equations identical to 5.11 and 5.14 were used for the hydrophobic (but the parameters used are dierent and listed in the table) 54 neighbors type type r type polar p 0.1 5 nonpolar np 0.1 7 membrane mem 6 d mem 2:5 Table 5.6: The parameters used in equation 5.12 for calculation of the number of neighboring residues. d mem refers to the membrane grid spacing. residue C np i N water i R in ALA -0.86 26 3 VAL -1.85 26 3.5 LEU -2.56 38 4 ILE -1.7 54 4.5 PRO -2.72 26 3.5 MET -1.7 26 3.5 TRP -3.7 24 6 PHE -2.83 54 5 Table 5.7: Parameters for Hydrophobic energy term parameter value dimensionality A 1 kcal 15 1/Angstrom r 0 2 Angstrom Table 5.8: The parameters for the hydrogen bond function 55 5.3 Eect of Electrolyte and Membrane Potential In trying to introduce the electrolytes, it is important to retain the option of multi- stage modeling by describing a key part of the system in a fully discrete model. For example, we can treat the immediate region near the protein by a primitive model (see gure 5.3) with the explicit ions moving by Langevin dynamics treatment [16]. While keeping this option in mind, we will focus below on the next layer of surround- ing where the electrolytes are described in a simplied way. Now, the surround- ing electrolytes can be treated macroscopically going back to the ideas behind the Debye-Huckel and the Gouy-Chapman models and subsequent analytical [30, 59] and numerical Poisson-Boltzmann (PB) models (e.g., ref [72] ). One may also follow the microscopic strategy that was used in MC studies [14, 101, 38, 54, 13, 69] and in more recent Langevin dynamics (LD) and molecular dynamics studies [117] . Here, however, we try to take a compromise between full MC and grid-type approaches. That is, we would like to start conceptually by placing ions on grid points (with an assumed separation) and then to use a MC model of type we and others applied in determining the ionization states in proteins [94, 10]. However, for practical purposes, we will try to introduce some additional sim- plications (which can also be removed). The rst possible simplication involves creating a grid for the ion positions and placing k1 positive ions and k2 negative 56 Figure 5.3: A schematic description of the treatment of the simulation system. The regions considered are (i) the protein membrane system (region I); (ii) the region of the solution and the ions that is treated with explicit ions and Langevin dynamics simulation (region II), which is not shown because this region is optional and is not being used in the present study; (iii) the region with explicit grid points (region III); (iv) the bulk region. The region between the bulk and the electrodes (region IV). The indices L and R stand for the left and right sides, respectively. 57 ions on n grid points by a MC procedure, leading to the following expression for the average charge at ith grid point: hq i i = P m q (m) i expG (m) P m expG (m) (5.26) wherehq i i is the average charge at the ith grid point, m designates the congura- tions of the charges on the grid points, and G (m) is the free energy of the mth conguration. The next simplication is similar to the TanfordRoxby model for ionization states in proteins [98] , using the self-consistent iterations with the po- tential from the other grid points instead of the MC congurational averaging [94]. This model is obtained by considering congurations with +1.0 and 0.0 charges for the positive ions, and 1.0 and 0.0 for the negative ions. Of course, we can still use the hybrid approach where the nearest neighbors are treated by MC. Next, we further simplify our approach by placing each ion pair on the same grid point (this simplication has been removed in some treatments). The above approximations basically led to the semi-macroscopic strategy applied in our previous electrostatic modeling [61], which is similar to the approach introduced originally by Klein and Pack [55] but it retains a more microscopic view. The grid spacing is taken here as and the volume element = 3 centered at the ith grid point contains a residual charge q g i determined by q g i =q + i +q i (5.27) 58 where q i = N box e ( i ) A = Q box e ( i ) A (5.28) where q + i and q i are, respectively, the positive and negative fractional charges that are assigned to the ith grid point, is the ion charge of the electrolyte ions in atomic units (namely,1 for the 1:1 electrolyte used in our calculations), N box box is the total number of cations/anions in the simulation box, Q box box is the total charge of cations/anions in the simulation system given byQ box = N box , i is the electrostatic potential (times unit charge) at the ith grid point, and = (kT ) 1 . A is a normalization constant. In the initial step of the calculation N box box is obtained as follows for a 1:1 electrolyte: N box =n bulk N box grid (5.29) where n bulk is number density (number of ions/Angstrom 3 ), which is connected to the molar concentration (C) by n bulk = C 1666 (5.30) boxN grid box is the number of grid points within the simulation box (outside membrane and protein). 59 Using the normalization condition P q i on Eq. 5.28, we obtain q i = N box e i P i e i = Q box e i P e i (5.31) Our normalization treatment, which is dierent than the common use of the bulk density, helps to obtain a clearer physical picture. However, for a very large simulation system and at the limit of ! 0, we obtain: A = X e i !N grid box (5.32) The potential times unit charge (which is actually potential of mean force) can be expressed as i = 332 X j q P i gp eff r ij + 332 X j q g i wat r ik +V ext i (5.33) whereV ext i represents the external potential on ith grid point. Here,q P j is the charge of the jth protein residue (these charges are evaluated by MC procedure described above) andq g k is the point charge at the kth grid point (representing the excess net charge of the kth volume element). The eect of the regions that are not included in the simulation system was represented by using a partially periodic (quasiperiodic) treatment. That is, we used a periodic boundary condition with explicit images of the simulation system in the nearest neighbor cells and with simplied replicas for farther images, where 60 the total charge of the slab at each Z coordinate is centered on the middle point of the slab (with coordinate (0;0;Z)). In principle, we could have used a treatment similar to that by Torrie and Valleau (18), but we felt more comparable with the current treatment (of course after validating it). The dielectric gp eff e represents chargecharge dielectric for interaction between the protein and the electrolyte grid points. This parameter can be approximated by a large number between 40 and 80 or by the distance dependent function (see ref. 5) = 1 + 60[1exp(0:1r ij )] (5.34) With the self-consistently evaluated set of the grid charges, we evaluate the eect of the ionic strengthG is by G is = 332 X ij q P i q g i gp eff r ij + 332 X k>j q g j q g j wat r ik (5.35) . In order to model the eect of the external potential, one can consider formally the membrane/ protein/water system as a capacitor. In this case, it is possible to use the well-known macroscopic capacitor model. That is, in principle we can use the well-known macroscopic capacitor model (e.g., refs. 32 and 33), where the external potential induces surface charges (f ), whose value will be dened below, 61 and creates the corresponding displacement vector D 0 (see ref. 32 and Fig. S1). In this case, we have D 0 = 4 f (5.36) where D 0 is not aected by the medium. We also have the relationship E =D 0 = =D 0 4P (5.37) where is the macroscopic local dielectric constant (which is not equal to the above ij ), E is the macroscopic eld, and P is the macroscopic polarization. Our rst task is to determine the membrane potential and the electrolytes charges, so we can evaluate the free energy of the protein charges in the presence of this potential, which is done by expressing the external potential as V i ext = Z z i z 0 D 0 z =(z)dz (5.38) whereZ0 is theZ coordinate at the left electrode (in the current work, we dene the left side as the side with smaller value for the Z coordinates and the right side with the larger Z value). Now, if we consider the potential on the electrolyte grid points or on the membrane grid points, we can dene what is considered here as the ideal potential (namely, the potential in regions without the protein) by 62 Figure 5.4: A schematic representation of the membrane/electrodes system (the membrane can also contain a protein). 63 V out i =V ideal i = 8 > > > > > > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > > > > > > : (zz min )D 0 = wat z <z ig 1 (z ig 1 z min )D 0 = wat + (zz ig 1 )D 0 = mem z ig 1 <=z <=z ig 2 (z ig 1 z min )D 0 = wat + (z ig 2 z ig 1 )D 0 = mem + (zz ig 2 )D 0 = wat z >z ig 2 (5.39) where z ig 1 is a left boundary of the membrane, z ig 2 is a right boundary of the membrane, andz min is the lowest (left-most) Z coordinate of the grid in the system (for this purpose we assume " at" membrane). Here wat and mem are the dielectric constants of water solution and the membrane/protein system, respectively. When we deal with the protein charged groups, we evaluate the relevant external potential by starting from the closest point, where we know the ideal potential of Eq. 5.39. In principle, because the potential is a state function, we can choose any path for evaluating the potential at the site of the given side chain. Here we follow the path depicted in Fig. 5.5, where we move up to z start pro and then evaluate the contribution to the potential at point i by replacing the integration over closely spaced points by a sum of the contributions from two points. With this approximation and our estimate of the dielectric at the protein sites, we obtain 64 V i ext = 8 > > > < > > > : V i ideal i +V shift D 0 (zz start pro ) 1 site (r wat ) 1 ideal jr wat center r wat j< 3 V (z start pro ) + D 0 (zcenterz start pro ) (r wat center ) + D 0 (zzcenter ) (r wat ) jr wat center r wat j>= 3 (5.40) where the dierent distances are dened in Fig. 5.5 and where the closest distance from ther center = (x;y;z center ) to the closest electrolyte grid point isr center wat . Here, the distance-dependent dielectric is dened as (r) = 34e (x=5) 2 + 6 (5.41) The shift potential, V shift , is based on the calculation of the grid potential V shift = V right V left 2 (5.42) Our next task is to have a clear microscopic (or semi-microscopic) idea of the energetic of the protein/membrane system in the presence of the electrolytes and the external potential, which is obtaining by writing G tot = G CG + 332 X ij q P i q g i gp eff r ij + 332 X k>j q g j q g j wat r ik + X j V ext q P j (5.43) where the leading term is the free energy in the absence of the external potential and the electrolytes. The use of equations above with a uniform dielectric is valid 65 Figure 5.5: Clarifying the considerations that led to the treatment of V i ext . The gure describes a division of the distance between the grid point Z start pro (where we have the "ideal potential") to the position of the given residue into two steps, which are then being used to evaluate the potential on the residues using a distance- dependent dielectric. 66 when we do not have a membrane between the electrodes because the eective dielectric can be taken as the water dielectric. However, the proper dielectric in the membrane region is less obvious. Thus we have used the D0 approach (with periodic electrolytes) for membrane-type systems. In such cases, one may wonder what should be the eect of the dielectric that will reproduce the potential from the electrolytes at the membrane. Our model has been validated by several test cases and subsequently explored by simulating a system of a neutral membrane introduced in the electrodeelectrolyte system, where we could reproduce the capacitance of the membrane. The calculations have demonstrated that the ions in the solution provide almost complete screen- ing of the external potential, which starts to increase only when it reaches the mem- brane boundaries. This behavior is similar to what is expected from macroscopic considerations. The present CG model focuses on a reliable description of the electrostatic free energy of the system, which includes re ecting the changes of the self-energy of the ionizable groups in the interiors of proteins and membrane. It is well-known that the self-energy can be reduced signicantly in nonpolar regions (as has been rst quantied by our early electrostatic studies, ref. 11) and that such desolvation penalty plays a major role in dierent processes, including proton transport (3437), ion translocations (7, 38), and other similar processes (5). The role of the desol- vation eects (and the corresponding controversy) in the case of ionized Arginines in membranes and membrane proteins has been widely discussed (e.g., refs. 3, and 67 3941). However, it has been established by Garca- Moreno and coworkers (42, 43) that the solvation penalty is less than what is usually assumed (re ecting, in fact, our long-standing view that, because of protein reorganization and water penetra- tion, the local dielectric in proteins is much larger than what is widely assumed; refs. 5, 11 and 44). These reduced solvation penalty has been carefully considered in implementing and calibrating our CG model, considering observed absolute folding energies and observed energetics of internal groups (6) and also considering (e.g., ref. 3) the results from simulations of water penetration eects in membranes and membrane protein interfaces (e.g., refs. 40 and 41). We would also like to point out that our CG ionization energy does not always correspond to the pKa obtained at the folded structure because of rather complex considerations described in ref. 4. Nevertheless, in most cases, the ionization state obtained by the CG model re ects the ionization state of the folded system and the calculated average eective charges are given in Table S5. The gating charge concept dates back to the Hodgkin and Huxley idea [40]. The simple form of this model assumes that the application of the membrane potential, V , increases the probability of the open conformation due to the free energy contri- bution Q gate V , associated with the displacement of the gating charge Qgate across the full membrane potential. In this case, we can write (see ref [39]) O C =e (G 0 +QgateV ) (5.44) 68 Residue # Residue Type hQ closed i hQ int i hQ open i 81 E*2 1.00 0.75 1.00 89 E*2 1.00 0.67 0.56 93 E*2 1.00 1.00 0.54 97 E*2 0.75 0.67 0.55 105 E*2 1.00 1.00 0.67 111 E*2 1.00 1.00 0.67 113 E*2 1.00 1.00 1.00 126 E*2 1.00 1.00 0.75 159 D*2 0.56 0.57 1.00 163 D*2 0.60 0.75 0.75 189 D*2 0.75 1.00 0.75 242 E*2 0.67 0.52 0.52 245 E*2 0.67 0.63 0.56 319 E*2 0.50 0.50 0.00 324 D*2 0.60 0.67 1.00 332 D*2 0.60 0.67 0.60 348 D*2 0.50 0.67 0.57 Table 5.9: The partial charges of ionized residues in a monomer of the Kv1.2 channel (only residues with dierent charges for dierent conformations are included ) 69 where O and C stand for the fraction of open and closed forms, respectively, and where G 0 is the free energy dierence between open and closed conformations at zero external potential. Of course, this model is oversimplied, but it will serve as a starting point. Now let us consider the case when a real xed (displacement) charge (Q d ) moves from point r op to r cl and the electrostatic potentials in corresponding points are V (r op ) and V (r op ) , respectively. In this case, we have O C =exp[G 0 (V (r op )Q d V (r cl )Q d )] =exp[G(V )] (5.45) where G(V ) corresponds to the free energy dierence between the closed and open congurations at the given applied potential. Now we can write G(V ) =G 0 + (V (r op )Q d V (r cl )Q d ) (5.46) If we deal with a real protein with many charges, we can only calculate G(V ) explicitly. However, if we assume a linear potential, we can write G(V ) =G 0 +V P i (z op i z cl i )q i L (5.47) Now using Eq 5.46 for zero potential and for potential V =V 0 , we can write G(V 0 ) G(0) =Q gate V 0 (5.48) and using Eq 5.47, we obtain 70 Q gate V 0 = P i (z op i z cl i )q i L (5.49) Obviously, such a picture would be justied if the potential across the protein/membrane system were obtained by converging microscopic simulations. However, in such cases, it is unlikely that we will have the same potential in dierent sites with the same Z value. Thus such a treatment does not provide a microscopic description even if a few of the relevant quantities (e.g., the average displacements) are eval- uated microscopically. A more microscopic strategy (in terms of the so-called G route of ref [83]) involves calculating the free energy change associated with each residue in the two congurations under the eect of the external potential. This approach does provide a fully microscopic treatment of the eect of the potential on the protein. However, it is not clear how accurate is the corresponding G-route free energy treatment because the convergence of applying such procedure to volt- age eects has never been established by reproducing well dened observed eects (namely pKa changes, redox changes, or spectral shifts). The diculty is to capture the equilibration of the electrolytes by the microscopic calculations. Furthermore, calculation of the charging free energy at dierent sites (e.g., after and before the conformational change) is extremely challenging, as has been demonstrated in our earlier studies of pKa in nonpolar protein interiors [29, 34]. At any rate, instead of using a linear potential in the evaluation of the gating charge, we used our CG model and evaluated the contributions of key residues; the corresponding results are summarized in Fig. 5.6. 71 Figure 5.6: The contribution from the crucial Arg residues to the energetics of the activation of the Kv1.2 channel. 72 As claried in the main text, our aim is to evaluate the gating charge and to have the option of evaluating the gating current by considering the electrolytes. In this respect, we note that Q gate V 0 = Z 0 (q=t)dt (5.50) where is the time where the gating current becomes zero (before the ion conduction starts) and where q=t is the amount of charge that passes the bulk region toward the electrode. In principle, we can calculate q=t by considering the change in the electrolyte charges in a response to the channel uctuations (evaluated by LD). Such a treat- ment may be accomplished by considering a time constant for this response in terms of the equilibration of the grid distribution upon change of the interaction with the protein potential. However, with the reasonable assumption that the electrolytes reequilibrate faster than the time for the change in the average protein structure, we can evaluate q(Z) at the beginning and the end of the protein motion, and the eective current will be given by the amount of charge that passes the bulk region in the time . In this case, we will get Q gate V 0 = Z 0 (q=t)dt Q = Z Z 0 inf (q grid (Z;V )=Z)dZ (5.51) 73 where q grid (Z;V ) is the dierence between the accumulative sum of q grid (Z;V ) of the open and closed channel. 5.4 Model Systems and Validation In order to establish the validity of our model we have examined it on several levels. First we have examined our ability to reproduce the DebyeHuckel (DH) analytical approximations for the charge distribution around an ion. More specically, as in the classical DH model we have represented the ion explicitly by a charge (+1 charge in this case) and the rest of the system was represented by our grid model. In order to get results which are independent of grid spacing and position we have generated several shifted and sparse grids (with a spacing of 5 between the grid points) and then averaged over the corresponding potential and charge distribution. The simulation results, which are summarized in Fig. 5.7, have reproduced the DH analytical results quite well. Next we explore our ability to model the charge distribution and potential of the electrolytes near a charged surface at a given ionic strength, which is usually described by the Gouy-Chapman (GC) approximation [66]. This type of problem has been the subject of extensive studies including PB, MC and MD [4150,52]. For this system we have built an explicit sheet of charges on a simulated electrode surface and generated our electrolyte grid on one side of the charged surface. The grid was then allowed to interact with the explicit charges of the electrode as well as with the charges of the dierent grid points, and the corresponding calculated 74 Figure 5.7: The potential around a charge obtained by our simulation and by the Debye Huckel (DH) theory as well as Coulombic (unscreened) potential. 75 Figure 5.8: The GouyChapman potential. potential is given in Fig. 5.8. The calculations have established the ability of our model to reproduce the GC screened electrostatic potential and the condition that the net charge of the electrolytes is equal to the total surface charge. We have also reproduced the correct decay of the potential at the Debye length (30 for the system under consideration). More specically the electrostatic potential becomes zero around 150 which corresponds to the analytical result (being approximately equal to ve times the Debye length). We have also explored the performance of our approach on a system composed of two electrodes and electrolytes. The corresponding results (Fig. 5.9) did provide 76 further major validation, reproducing the trend obtained in MC simulations of the related systems [46]. Next we have explored the performance of our model for a system with two electrodes, electrolyte solution and a neutral membrane. The calculated proles for the potential along the Z coordinate,where the external potential is 100 mV, are depicted in Fig. 5.10 a. We have also evaluated the behavior of the charge distribution and presented in Fig. 5.10 b the corresponding results. An important requirement for the validation of our (or any other) model is the ability to reproduce the capacitance of the membrane system. Examining the performance of the present model we found out that the computed capacitance of our membrane system agrees quite well with the corresponding analytical expression of Lauger [40]. Next we consider the eect of moving charges in the membrane region, as a preparation for exploring the corresponding displacement charges in Kv1.2. This was done by rst applying an external voltage of 0.1 V, placing a charge of +16e near the surface of the membrane (at the coordinate (0, 0,10)) and letting the electrolytes equilibrate separately at each side (Zb0 and ZN0) without any transfer across the membrane. Subsequently we moved the charge (+16e) to the other side of the membrane (at the coordinate (0, 0, 10)) and let the electrolytes equilibrate in the same way as before. We dene the gating charge as the change in the accumulated electrolyte charges, on either side of the membrane, upon movement of the (protein) charge inside the membrane (a more quantitative denition will be 77 Figure 5.9: (a) The electrolyte potential (in Volt) and (b) the electrolyte charge distribution along Z-axis for a system with two electrodes (external potential 100 mV). We have shown the results for two dierent electrolyte concentrations: 10 mM (red) and 100 mM (blue). Fig. 5(a) also shows the linear electrode potential (in absence of electrolyte) as reference (black). These calculations have been performed with (17X17) periodic images in the XY plane (see Fig. 7 for convergence study). 78 Figure 5.10: (a) The electrostatic potential (in kcal/mol) and (b) electrolyte charge distribution along Z-axis for a system with two electrodes (external potential 100 mV), 40 A thick membrane at the center of the box (marked by red vertical lines) and electrolyte with concentration 100 mM. The present case required the use of (17x17) periodic images (in the XY plane) to reach convergence. For this system, the computed capacitance is 2:2 10 3 CV 1 m 2 and the analytical capacitance is 4:3 10 3 CV 1 m 2 . 79 Figure 5.11: A simplied model for exploring the gating charges. The model involves a membrane and positive charges that are moved from one side to the other. given in Eq. (42)). The computed results for such a system is described in Fig. 5.12 and the gating charge for the described conditions are found to be 5. 5.5 Simulation of Kv1.2 As stated in the previous sections, our challenge is to model the energy balance in voltage activation channels. Thus we took as a test system the Kv1.2 voltage- activated channel and started our study by evaluating the potential in the open and closed models of the protein/membrane system. The calculations started with the structural models built by Rosetta (with structural information on the open state; ref. [64]) and used by Pathak et al. [77] for both the open and closed structures In each case, we started by determining the ionization state of the protein ionizable groups using the CG model and our constant pH simulation approach (see refs. [26] 80 Figure 5.12: The accumulative integration of the charge proles for the model system described in Fig. 8 when the+16e charge resides at Z=+10 (black) and Z=-10 (red). Here we dene the gating charge as the change in the accumulated electrolyte charges (the dierence of the maxima in this gure), on either side of the membrane, upon movement of the (protein) charge inside the membrane (see Eq. 42). The gating charge for this system is found to be 5. 81 Figure 5.13: The potential in the Kv1.2 system in the open and closed states. and [49]). We started the calculations with initial Monte Carlo (MC) procedure for determination of the initial ionization states of the protein. We then introduced the electrolytes and the external potential and allowed the complete system to equilibrate in terms of the electrolyte distribution and the ionization states (which was done by performing the MC procedure at equal intervals of the electrolyte equilibration steps). The calculated potential prole is depicted in Fig. 5.13 for the closed and open channel. 82 The change in potential occurs mainly in the membrane protein region, establish- ing the fact that our electrolyte model provides basically innite dielectric screening in the bulk region. The fact that we obtained the screening (without assuming it) is encouraging because we are using a model with discrete features, where it is not guaranteed that we will get the physically correct trend. The gating charge is a crucial parameter that represents the shift of the rela- tive free energy dierence between the closed and open congurations, due to the change in an external potential. This parameter, which was initially postulated by Hodgkin and Huxley [40], provides a qualitative explanation of the coupling of the external potential to the channel activation. The evaluation of the gating charge is usually done in an indirect way, using reasonable but not necessarily microscopic assumptions. One may simply determine the relative population of open and closed channels as a function of the applied potential and determine the Boltzmann prob- ability for the voltage-induced structural change [39]. This treatment is equivalent to the assumption that the energy needed to move the gating charge, Q gate , in the applied electric eld is equal to the work of moving the protein charges between the two congurations under the membrane electric eld. This assumption can be formulated as Q gate V = G cl!op (5.52) where G cl!op is the contribution of the membrane potential to the work of moving from the closed to open conguration, and V is the change in the electrostatic potential between the initial and nal position of the protein eective charge . Under 83 the assumption of linear membrane potential, it is simple to calculate Q gate if the structures of the open and closed states are known, and a useful related insight has been obtained from macroscopic studies [77, 35], . An attempt to evaluateQ gate with the philosophy of Eq. 5.52 using microscopic simulations has also been reported [53] . However, there are some problems with the current strategies, including the fact that the potential is not really linear or uniform in the protein membrane system (in contrast to the implicit assumption of one of the versions of the treatment of ref [53]). At any rate, we must consider the fact that the real observables are G cl!op and the actual measured gating current (obtained as discussed in, e.g., ref. [93]). In other words, despite the remarkable concep- tual importance of the gating charge obtained from Eq. 5.52, modeling approaches should arguably focus on reproducing G cl!op and the experimentally measured in- tegrated gating current. Thus we focused on the actual measured quantity, namely, the gating current, rather than on its interpretation. Our explicit calculations of the gating charge was done by rst applying the potential in the closed structure, letting the electrolytes equilibrate, and then xing the number of positive and neg- ative ions on both sides of the membrane. Next we allowed the protein to move from the closed to the open structure, while allowing the electrolytes to equilibrate with the crucial constraint that they cannot pass through the channel. In considering the gating current or its integral over time, we note that the move- ment of the positively charged protein residues to theZ direction leads to movement of negative solution ions toward the membrane protein system (this polarization is 84 clearly captured by the PB treatments). However, the measured current is actu- ally due to the movement of the compensating positively charged ions toward the electrode. Obviously, this current has the same magnitude (with opposite sign) as the negative current that moved toward the membrane, but the physics is better described by considering the actual gating current. Thus we evaluate the gating charges by Q gate (V )int Z 0 inf (q grid (Z;V ))dZ (5.53) , where q grid (Z;V ) is the dierence between the accumulative sum of q grid (Z;V ) of the open and closed channels, and Z 0 is the point to the left of Z 0 , where the electrolyte charge distribution near the membrane changes sign. At this point, the integrated charge reaches a plateau and then starts to decrease. Note that the in- tegral evaluates the charges generated by the current after it equilibrates on the left side of the membrane, but before it actually penetrates the membrane. The is shown in Fig. 5.14. The value of Q gate appears to depend on the dielectric con- stant used for the interaction between the electrolyte grid points and the protein charges. Here we have found that using an eective dielectric constant of 40 and 80 gives gating charge (Q gate ) of 10e and 5e, respectively, whereas the observed value is 1214e [93, 2]. A more careful comparison to the experimentally determined value of Q gate should involve simulation of the integrated current at dierent potentials and consideration of the eect of the landscape on the conformational transition. 85 Figure 5.14: The gating charge for the Kv1.2 system. The gure presents the integrated charge . The dierence in the maximum accumulated charge between two states is taken as the gating charge, which is approximately 10 in this gure. 86 Let's now concentrate on energetics of the voltage activation process. At present, it is extremely challenging to explore the conformational landscape by explicit sim- ulations. Alternatively, we expect to obtain semiquantitative information by our CG model (e.g., see our recent work in ref. [74]). Here we have focused on the energetics at the end points (namely, closed and open structures) at dierent values of the external potential. We have also examined the energetics of an intermediate structure, which was done by pushing the system from the closed to open struc- ture, using targeted molecular dynamics, locating intermediate structures, and then evaluating the CG free energy at the corresponding region. The calculated results are summarized in Table 1. The most important nding that has emerged from the calculations is that we succeed in reproducing the eect of moving the channel from its closed to open conguration upon change of the external potential from a nega- tive to positive value. The overall change in energy is about 30 kcalmol (+14 and 16 kcalmol for 200 and +200 mV, respectively), re ecting both the conformational energy and the eect of the external potential. Most signicantly, it appears that our model gives stable results and that the trend obtained in Table 5.10 is retained even when we change the parameters in the model in a reasonable range. For ex- ample, changing the grid protein dielectric from 80 to 40 results in overall change of energy (upon change of potential from 200 to +200 mV) of about 40 instead of 30 kcalmol. It is also important to note that fully microscopic calculations are expected to produce (at present) errors of more than 100 kcalmol (see discussion of 87 the related simulations of F1-ATPase, ref. [74]). Omitting the hydrophobic contri- butions lead to an even smaller change in trend than the one obtained with change of dielectric. Voltage, volt Closed Intermediate Open -0.2 -203 -201 -189 -0.1 -201 -165 -193 0.0 -201 -130 -200 0.1 -201 -6 -208 0.2 -202 -65 -218 Table 5.10: The free energy at dierent regions of the conformation/voltage land- scape (The eect of the change in potential in each structure was adjusted, sub- tracting the capacity work in the open structure from the energy of all the other structures.) The challenge of evaluating the barrier for the conformational transition has been addressed here in a preliminary way, considering only the intermediate conguration discussed above. As seen from Table 5.10, the barrier at the intermediate structure increases when the membrane potential becomes positive. This trend appears to be stable (with regard to the dielectric used) and associated with contributions from many residues rather than just a few. The height of the free energy barrier, however, is signicantly larger than the experimental estimate of about 5 kcalmol at V = 0 [96]. The calculated barrier can be reduced by several factors, including conformational pathway where the gating Arg residues move closer to the protein and even by a path where the four subunits of the tetramer move sequentially, 88 rather than moving together. The actual height of the barrier and the exact nature of landscape should be evaluated by more detailed mapping and this can be done by our renormalization approach [49]. The general trend of the calculations was used to generate the tentative surface by constructing an empirical valence bond-type conformational surface (e.g., ref. [81]) with the calculated G cl!op at dierent voltage values and with the experi- mental estimate for the activation barrier (5 kcalmol at V = 0), modulated by a dependence on the potential which follows the calculated trend (which is, however, drastically scaled down). This drastic modication of the barrier re ects the need for more careful mapping at the intermediate region, including the need to allow the system to nd the least energy path with the applied potential. One of the most challenging issues in modeling ion channels is the ability to provide a microscopic description of the time dependence of the gating charge. Here one can use the relevant experiments in constructing eective models (e.g., ref [96]). However, such models cannot be related directly to the underlying structures and also are not likely to be unique. Although our study has not yet produced the complete landscape, it is instructive to demonstrate here the potential of the use of a calculated landscape in simulating the time-dependence relationship between the gating current and the applied voltage. Such a demonstration is provided here by running Langevin dynamics simulations based on the eective landscape of Fig. 5.15 (at V = 0), while using the eective friction determined by our renormalization approach in a similar system [81]. The calculated time dependence of the eective 89 Figure 5.15: Simulating the gating uctuations. The gure describes the results of the Langevin dynamics simulations on the tentative surface of Fig. 5.1 for V = 0. The simulations were done with a friction of 150 ps 1 . conformational coordinate can then be translated to current by assuming that the electrolytes respond rapidly to the protein charge displacement . The simulated result presented in Fig. 5.15 provides a rst glimpse of the origin of the uctuations of the gating current. The phenomenological pioneering studies of ref [96] have addressed the same issue, but have not tried to use structure-based simulations. The above analysis can be further advanced in a more realistic direction, which can range from establishing a clear correlation between the reaction coordinate and the combined eective position of the protein charges, to simulating the time dependence of the uctuations of the CG model. Further insight can be obtained 90 Figure 5.16: A conceptual gure of the activation of the Kv1.2 and related channels. The gure considers the allosteric eect of the potential, where the change of posi- tion of the positive gating charges reduces the repulsion between theses charges and the positive elements of the gating region, allowing the attached nonpolar region to open and to reduce the desolvation penalty of a transferred ion. by using our full renormalization approach [27] to simulate the time dependence of the system under fast changing potential pulse. 5.6 Discussion This work used a recently developed CG model for studies of the eect of external potentials on membrane proteins and explored the utility of such a model in simu- lating the eect of membrane potential on ion channels. The model used involves our early CG model of the protein membrane system and a dierent model of the electrolyte solution and the external electrodes. The electrolyte model allows one to navigate between the more microscopic MC modeling to faster mean eld models and helps in providing some light on various elements that are needed in order to understand the nature of the external potential. Applying the CG model to study 91 the energetics of the voltage-activated Kv1.2 channel provided several advances. These include the evaluation of the balance between the protein conformational energy and the applied potential, while obtaining the correct trend in this balance without the need of any specically adjusted parameter. The CG simulations also provided a clearer nonstandard description of the gating current and gating charge, allowing one to look directly at the change in the electrolyte charges rather than the interaction between the linearized external potential and the protein charges. Finally, our study has introduced a landscape-based simulation of the gating uc- tuations and basically the gating current. At this point, it may be useful to point out that, although it is interesting to understand the nature of the gating charge and the corresponding contribution from dierent residues, it is arguably more important to analyze the energy contribution of each residue . That is, to un- derstand the interplay between the external potential and the protein/membrane conformational energy, we have to see what are the free energy contributions of the dierent residues. Now, although the main eects are electrostatic, they re ect the interactions between the residue charges and their surroundings (plus the external electrode potential), rather than just the interaction between the charges and the linearized potential. One of the most important questions regarding the present model is its stability in terms of the crucial ionization states and the overall ener- getics. We will only mention here that we believe that the current model provides one of the most reliable, comprehensive, and stable treatments of the electrostatic energies in protein/ membrane systems, and we consider this important issue as 92 well as the ionization states of the protein residues in the SI Text. The CG model may appear to some to represent an oversimplication that overlooks the many substates of the real surface and thus the origin of the uctuating current. In fact, our renormalization approach generates an eective friction that should represent the uctuations between the dierent substates (see ref. 18). Furthermore, the eective CG free energy does include the entropic contributions to the overall free energy. Thus we expect that more careful landscape mapping will be able to ad- dress questions about the role of entropic eects (e.g., ref. 24). Another important issue is the control of ion selectivity, which is most probably due to change in ionion interaction in multiion conductance (6). Here it would be important to explore the eect of the uctuations in the substates of the open state. The present work also provides important insight on the overall function of voltage-activated channels. That is, as described in Fig. 6, the opening of the channel can be viewed as moving a positive charge from a region near the gate to a region farther away, in response to the applied potential. In fact, viewing the voltage activation as electrostatic allosteric process of interplay between the interaction between the applied potential and the protein positive groups, which are in turned interacting (indirectly) with the transferred charge, seems to be a useful generalization. That is, in the closed state, the positive protein charges are closer to the site where the conducted ions have a higher barrier, where the actual barrier for the ion transfer is due to the local nonpolar environment in the closed gate (see ref. 10). When the positive charges move away, following activation, the indirect interaction with the transferred ion is 93 reduced. This change in interaction involves the opening of the nonpolar region, which is, however, due to the change of its interaction with the positive charges. In other words, the force that moves the nonpolar region to its location is the in- teraction between the gate region and the gating charges (see a related case in the function of hemoglobin in ref. 26). Regardless of the general qualitative insight, the key point of the present study is in highlighting the ability to use initial structural models and to explore their consistency with the functional constraints about the energy balance and the remarkable sets of experimental constraints. In particu- lar, a more rened CG structure- based energy landscape should be very useful in exploring the molecular meaning of the current uctuations with dierent applied voltage pulses and their relationship to the biological function of voltage-activated channels. Furthermore, the model should provide an eective way for extracting detailed molecular information from dierent mutational studies. 94 Chapter 6. Conclusions 6.1 Main Achievements The results of current work are discussed in the following papers: Arieh Warshel, Anatoly Dryga. Simulating electrostatic energies in proteins: Perspectives and some recent studies of pKas, redox, and other crucial functional properties Proteins, Vol. 79, No. 12. (2011), pp. 3469-3484. Shina C. L. Kamerlin, Spyridon Vicatos, Anatoly Dryga, Arieh Warshel. Coarse- Grained (Multiscale) Simulations in Studies of Biophysical and Chemical Systems. Annual Review of Physical Chemistry, Vol. 62, No. 1. (2011), pp. 41-64. Anatoly Dryga, Suman Chakrabarty, Spyridon Vicatos, Arieh Warshel. Realistic simulation of the activation of voltage-gated ion channels. Proc Natl Acad Sci USA. 2012 Feb 13 (in press) Anatoly Dryga, Suman Chakrabarty, Spyridon Vicatos, Arieh Warshel Coarse grained model for exploring voltage dependent ion channels. Biochim Biophys Acta. 2012 Feb;1818(2),pp. 303-17. Anatoly Dryga and Arieh Warshel. Renormalizing SMD: The Renormalization Approach and Its Use in Long Time Simulations and Accelerated PMF Calculations of Macromolecules. J. Phys. Chem. B, 2010, 114 (39), pp. 1272012728. The main focus of the current dissertation is a simulation of ion channels and de- velopment of a methodology to overcome current issues in modeling of such type of systems. More specically gramicidin A was used as a model channel and Kv1.2 was 95 used as a large voltage-activated channel. To overcome current diculties of brute- force molecular dynamics a renormalization approach was developed and validated for gramicidin A channel. It was shown that the renormalization approach is eec- tive for the system with a known reaction coordinate and is capable of accelerating the simulation up to two orders of magnitude. The renormalization approach also allowed us to investigate the eect of mutations on free energy proles in channels thereby enabling simulations of challenging conformation changes in many systems. The current work extends CG model developed in [73]. The description of interaction of the membrane [88] and protein is improved and parameters are rened. A semi-microscopic treatment of electrolyte solution is introduced and validated for the system with known results. And interaction of the external potential and protein charges are described on semiquantitative level. These developments allowed us to simulate the gating process in the Kv1.2 channel based on the molecular model and arguably presents the rst demonstration of stability calculations for open and closed conformations of the channel. Combining our renormalization approach and CG model allowed us to obtain preliminary results for the gating mechanism. We performed renormalization stud- ies of the movement from open to closed structures for a stepwise and concerted mechanism. Our simulations indicates that stepwise mechanism is more likely to occur than the concerted mechanism. The nature and the process of the gating is assumed to be due to S4 Arginines interactions with charged residues in S5 and S6 units. 96 The overall current work demonstrated the power of renormalization and the CG approach for the studies of the complex biomolecular systems and particularly ion channels. We feel that combining these two powerful approaches opens the way to investigate many properties which cannot be performed at present or in the near future by brute-force approaches. 6.2 Future Work The model developed in the current work opens numerous possibilities of future progress in molecular-based simulations of ion channels and proteins. In particular the renormalization approach allows us to simulate ion current and selectivity of large channels with a reasonable CPU time. This may signicantly improve our understanding of the selectivity mechanism in dierent classes of ion channels. It is also possible to investigate the mechanism of the gating to get better un- derstanding of it and to devise ways to aect the gating process. With the current version of the CG model it is possible to simulate the eect of the binding of small molecules and its eect on the stability of the ion channel conformations. The devel- oped renormalization approach and rened CG model also enable the investigation of many intricate biophysical processes. We feel that the methodology will be of great use for translocon systems and for conformational changes in many complex systems. It should also be emphasized that methodology developed here does re- quire modest computational resources and the results of simulations are very often 97 of comparable quality to those of fully-atomistic simulations but without the need for extreme computational power to overcome convergence problems. 98 Bibliography [1] SA Adelman. Fokker-Planck equations for simple non-markovian systems. Journal of chemical physics, 64(1):124{130, 1976. [2] SK Aggarwal and R MacKinnon. Contribution of the S4 segment to gating charge in the Shaker K+ channel. Neuron, 16(6):1169{1177, 1996. [3] DJ Aidley and PR Staneld. Ion Channels: Molecules in Action. Cambridge University Press, 1 edition, 1996. [4] TW Allen, OS Andersen, and B Roux. Energetics of ion conduction through the gramicidin channel. Proceedings of the national academy of sciences of the United States of America, 101(1):117{122, 2004. [5] J Aqvist and V Luzhkov. Ion permeation mechanism of the potassium channel. Nature, 404(6780):881{884, 2000. [6] J Aqvist and A Warshel. Energetics of ion permeation through membrane channels - solvation of Na+ by gramicidin-A. Biophysical journal, 56(1):171{ 182, 1989. [7] WL Ash, MR Zlomislic, EO Oloo, and DP Tieleman. Computer simulations of membrane proteins. Biochimica et biophysica acta-biomembranes, 1666(1- 2):158{189, 2004. [8] FM Ashcroft. Ion Channels and Disease. Academic Press, 1 edition, 1999. [9] GS Ayton, WG Noid, and GA Voth. Multiscale modeling of biomolecular systems: in serial and in parallel. Current opinion in structural biology, 17(2):192{198, 2007. [10] D Bashford and M Karplus. Pkas of ionizable groups in proteins - atomic detail from a continuum electrostatic model. Biochemistry, 29(44):10219{ 10225, 1990. [11] N BenTal, A BenShaul, A Nicholls, and B Honig. Free-energy determinants of alpha-helix insertion into lipid bilayers. Biophysical journal, 70(4):1803{1812, 1996. 99 [12] S Berneche and B Roux. Energetics of ion conduction through the K+ channel. Nature, 414(6859):73{77, 2001. [13] D Boda, KY Chan, and D Henderson. Monte Carlo simulation of an ion- dipole mixture as a model of an electrical double layer. Journal of chemical physics, 109(17):7362{7371, 1998. [14] D Boda, WR Fawcett, D Henderson, and S Sokolowski. Monte Carlo, density functional theory, and Poisson-Boltzmann theory study of the structure of an electrolyte near an electrode. Journal of chemical physics, 116(16):7170{7176, 2002. [15] A Burykin, M Kato, and A Warshel. Exploring the origin of the ion selectivity of the KcsA potassium channel. Proteins-structure function and genetics, 52(3):412{426, 2003. [16] A Burykin, CN Schutz, J Villa, and A Warshel. Simulations of ion current in realistic models of ion channels: The KcsA potassium channel. Proteins- structure function and genetics, 47(3):265{280, 2002. [17] CP Calderon. On the use of local diusion models for path ensemble aver- aging in potential of mean force computations. Journal of chemical physics, 126(8):78{92, 2007. [18] CP Calderon, L Janosi, and I Kosztin. Using stochastic models calibrated from nanosecond nonequilibrium simulations to approximate mesoscale infor- mation. Journal of chemical physics, 130(14):781{792, 2009. [19] WA Catterall. Structure and function of voltage-sensitive ion channels. Sci- ence, 242(4875):50{61, 1988. [20] A Cha, GE Snyder, PR Selvin, and F Bezanilla. Atomic scale movement of the voltage-sensing region in a potassium channel measured via spectroscopy. Nature, 402(6763):809{813, 1999. [21] S Chandrasekhar. Stochastic problems in physics and astronomy. Reviews of modern physics, 15(1):0001{0089, 1943. [22] A Chetwynd, CL Wee, BA Hall, and MSP Sansom. The Energetics of Transmembrane Helix Insertion into a Lipid Bilayer. Biophysical journal, 99(8):2534{2540, 2010. 100 [23] WF DeGrado, CM Summa, V Pavone, F Nastri, and A Lombardi. De novo design and structural characterization of proteins and metalloproteins. Annual review of biochemistry, 68:779{819, 1999. [24] S Dorairaj and TW Allen. On the thermodynamic stability of a charged arginine side chain in a transmembrane helix. Proceedings of the national academyofsciencesoftheUnitedStatesofAmerica, 104(12):4943{4948, 2007. [25] DA Doyle, JM Cabral, RA Pfuetzner, AL Kuo, JM Gulbis, SL Cohen, BT Chait, and R MacKinnon. The structure of the potassium channel: Molec- ular basis of K+ conduction and selectivity. Science, 280(5360):69{77, 1998. [26] A Dryga, S Chakrabarty, S Vicatos, and A Warshel. Coarse grained model for exploring voltage dependent ion channels. Biochimica et biophysica acta- biomembranes, 1818(2):303{317, 2012. [27] A Dryga and A Warshel. Renormalizing SMD: The Renormalization Approach and Its Use in Long Time Simulations and Accelerated PMF Calculations of macromolecules. Journal of Physical Chemistry B, 114(39):12720{12728, 2010. [28] SR Durell, YL Hao, and HR Guy. Structural models of the transmembrane region of voltage-gated and other K(+) channels in open, closed, and inacti- vated conformations. Journal of structural biology, 121(2):263{284, 1998. [29] JJ Dwyer, AG Gittis, DA Karp, EE Lattman, DS Spencer, WE Stites, and B Garcia-Moreno. High apparent dielectric constants in the interior of a protein re ect water penetration. Biophysical journal, 79(3):1610{1620, 2000. [30] CT Everitt and DA Haydon. Electrical capacitance of a lipid membrane separating 2 aqueous phases. Journaloftheoreticalbiology, 18(3):371{&, 1968. [31] JD Faraldo-Gomez, GR Smith, and MSP Sansom. Setting up and optimiza- tion of membrane protein simulations. European biophysics journal with bio- physics letters, 31(3):217{227, 2002. [32] M Feig. Kinetics from implicit solvent simulations of biomolecules as a func- tion of viscosity. Journal of chemical theory and computation, 3(5):1734{1748, 2007. [33] PL Freddolino, F Liu, M Gruebele, and K Schulten. Ten-microsecond molec- ular dynamics simulation of a fast-folding WW domain. Biophysical journal, 94(10):L75{L77, 2008. 101 [34] GB Goh, Garcia-Moreno EB, and CL Brooks III. The high dielectric constant of staphylococcal nuclease Is Encoded in its structural architecture. Journal of the american chemical society, 133(50):20072{20075, 2011. [35] M Grabe, HC Lai, M Jain, YN Jan, and LY Jan. Structure prediction for the down state of a potassium channel voltage sensor. Nature, 445(7127):550{553, 2007. [36] M Grabe, H Lecar, YN Jan, and LY Jan. A quantitative assessment of models for voltage-dependent gating of ion channels. Proceedings of the national academy of sciences of the United States of America, 101(51):17640{17645, 2004. [37] J Gullingsrud and K Schulten. Gating of MscL studied by steered molecular dynamics. Biophysical journal, 85(4):2087{2099, 2003. [38] D Henderson and D Boda. Insights from theory and simulation on the elec- trical double layer. Physical chemistry chemical physics, 11(20):3822{3830, 2009. [39] B Hille. Ion Channels of Excitable Membranes. Sinauer Associates, 3 edition, 2001. [40] AL Hodgkin and AF Huxley. A quantitative description of membrane cur- rent and its application to conduction and excitation in nerve. Journal of physiology-London, 117(4):500{544, 1952. [41] K Hukushima and K Nemoto. Exchange Monte Carlo method and application to spin glass simulations. Journal of the physical society of Japan, 65(6):1604{ 1608, 1996. [42] W Im, M Feig, and CL Brooks. An implicit membrane generalized born theory for the study of structure, stability, and interactions of membrane proteins. Biophysical journal, 85(5):2900{2918, 2003. [43] B Isralewitz, M Gao, and K Schulten. Steered molecular dynamics and mechanical functions of proteins. Current opinion in structural biology, 11(2):224{230, 2001. [44] C Jarzynski. Equilibrium free-energy dierences from nonequilibrium mea- surements: A master-equation approach. Physical review e, 56(5, Part a):5018{5035, 1997. 102 [45] C Jarzynski. Nonequilibrium equality for free energy dierences. Physical review letters, 78(14):2690{2693, 1997. [46] YX Jiang, A Lee, JY Chen, M Cadene, BT Chait, and R MacKinnon. The open pore conformation of potassium channels. Nature, 417(6888):523{526, 2002. [47] YX Jiang, A Lee, JY Chen, V Ruta, M Cadene, BT Chait, and R MacKinnon. X-ray structure of a voltage-dependent K+ channel. Nature, 423(6935):33{41, 2003. [48] V Jogini and B Roux. Dynamics of the Kv1.2 voltage-gated K(+) channel in a membrane environment. Biophysical journal, 93(9):3070{3082, 2007. [49] SCL Kamerlin, S Vicatos, A Dryga, and A Warshel. Coarse-grained (mul- tiscale) simulations in studies of biophysical and chemical systems. Annual Review Physical Chemistry, Vol 62, 62:41{64, 2011. [50] M Karplus and JA McCammon. Molecular dynamics simulations of biomolecules. Nature structural biology, 9(9):646{652, 2002. [51] M Kato and A Warshel. Through the channel and around the channel: Vali- dating and comparing microscopic approaches for the evaluation of free energy proles for ion penetration through ion channels. Journal of physical chem- istry b, 109(41):19516{19522, 2005. [52] JNC Kew and CH Davies. Ion channels: from structure to function. Oxford University Press, 1 edition, 2010. [53] F Khalili-Araghi, V Jogini, V Yarov-Yarovoy, E Tajkhorshid, B Roux, and K Schulten. Calculation of the gating charge for the Kv1.2 voltage-activated potassium Channel. Biophysical journal, 98(10):2189{2198, 2010. [54] K Kiyohara and K Asaka. Monte Carlo simulation of electrolytes in the constant voltage ensemble. Journal of chemical physics, 126(21), 2007. [55] BJ Klein and GR Pack. Calculations of the spatial-distribution of charge- density in the environment of DNA. Biopolymers, 22(11):2331{2352, 1983. [56] R Kubo. Fluctuation-dissipation theorem. Reports on progress in physics, 29(Part 1):255{&, 1966. 103 [57] S Kumar, PW Payne, and M Vasquez. Method for free-energy calculations using iterative techniques. Journal of computational chemistry, 17(10):1269{ 1275, 1996. [58] RA Kumpf and DA Dougherty. A mechanism for ion selectivity in potas- sium channels - computational studies of cation-pi interactions. SCIENCE, 261(5129):1708{1710, 1993. [59] P Lauger, W Lesslaue, E Marti, and J Richter. Electrical properties of bi- molecular phospholipid membranes. Biochimica et biophysica acta, 135(1):20{ &, 1967. [60] RJ Law, RH Henchman, and JA McCammon. A gating mechanism proposed from a simulation of a human alpha 7 nicotinic acetylcholine receptor. Pro- ceedings of the national academy of sciences of the United States of America, 102(19):6813{6818, 2005. [61] FS Lee, ZT Chu, and A Warshel. Microscopic and semimicroscopic calcula- tions of electrostatic energies in proteins by the polaris and enzymix programs. Journal of computational chemistry, 14(2):161{185, 1993. [62] XF Li, SA Hassan, and EL Mehler. Long dynamics simulations of proteins using atomistic force elds and a continuum representation of solvent eects: Calculation of structural and dynamic properties. Proteins-structure function and bioinformatics, 60(3):464{484, 2005. [63] E Lindahl, B Hess, and D van der Spoel. GROMACS 3.0: a package for molecular simulation and trajectory analysis. Journal of molecular modeling, 7(8):306{317, 2001. [64] SB Long, EB Campbell, and R MacKinnon. Crystal structure of a mammalian voltage-dependent Shaker family K+ channel. Science, 309(5736):897{903, 2005. [65] SB Long, X Tao, EB Campbell, and R MacKinnon. Atomic structure of a voltage-dependent K+ channel in a lipid membrane-like environment. Nature, 450(7168):376{U3, 2007. [66] Q Lu and J Wang. Single molecule conformational dynamics of adenylate kinase: Energy landscape, structural correlations, and transition state en- sembles. Journal of the american chemical society, 130(14):4772{4783, 2008. 104 [67] VB Luzhkov and J Aqvist. K+/Na+ selectivity of the KcsA potassium chan- nel from microscopic free energy perturbation calculations. Biochimica et biophysica acta-protein structure and molecular enzymology, 1548(2):194{202, 2001. [68] SJ Marrink, HJ Risselada, S Yemov, DP Tieleman, and AH de Vries. The MARTINI force eld: Coarse grained model for biomolecular simulations. Journal of physical chemistry B, 111(27):7812{7824, 2007. [69] A Martin-Molina, C Calero, J Faraudo, M Quesada-Perez, A Travesset, and RH Alvarez. The hydrophobic eect as a driving force for charge inversion in colloids. Soft matter, 5(7):1350{1353, 2009. [70] B Martinac. Mechanosensitive ion channels: molecules of mechanotransduc- tion. Journal of cell science, 117(12):2449{2460, 2004. [71] A Masunov and T Lazaridis. Potentials of mean force between ionizable amino acid side chains in water. Journal of the american chemical society, 125(7):1722{1730, 2003. [72] S Mclaughlin. The electrostatic properties of membranes. Annual review of biophysics and biophysical chemistry, 18:113{136, 1989. [73] BM Messer, M Roca, ZT Chu, S Vicatos, A Vardi Kilshtain, and A Warshel. Multiscale simulations of protein landscapes: Using coarse-grained models as reference potentials to full explicit models. Proteins-structure function and bioinformatics, 78(5):1212{1227, 2010. [74] S Mukherjee and A Warshel. Electrostatic origin of the mechanochemical rotary mechanism and the catalytic dwell of F1-ATPase. Proceedings of the national academy of sciences of the United States of America, 108(51):20550{ 20555, 2011. [75] JT Padding and AA Louis. Hydrodynamic interactions and Brownian forces in colloidal suspensions: Coarse-graining over time and length scales. Physical review E, 74(3), 2006. [76] S Park and K Schulten. Calculating potentials of mean force from steered molecular dynamics simulations. Journal of chemical physics, 120(13):5946{ 5961, 2004. [77] MM Pathak, V Yarov-Yarovoy, G Agarwal, B Roux, P Barth, S Kohout, F Tombola, and EY Isaco. Closing in on the resting state of the shaker K+ channel. Neuron, 56(1):124{140, 2007. 105 [78] DA Pearlman, DA Case, JW Caldwell, WS Ross, TE Cheatham, S Debolt, D Ferguson, G Seibel, and P Kollman. Amber, a package of computer- programs for applying molecular mechanics, normal-mode analysis, molecular- dynamics and free-energy calculations to simulate the structural and energetic properties of molecules. Computer physics communications, 91(1-3):1{41, 1995. [79] PS Phale, A Philippsen, C Widmer, VP Phale, JP Rosenbusch, and T Schirmer. Role of charged residues at the OmpF porin channel constric- tion probed by mutagenesis and simulation. Biochemistry, 40(21):6319{6325, 2001. [80] JC Phillips, R Braun, W Wang, J Gumbart, E Tajkhorshid, E Villa, C Chipot, RD Skeel, L Kale, and K Schulten. Scalable molecular dynamics with NAMD. Journal of computational chemistry, 26(16):1781{1802, 2005. [81] AV Pisliakov, J Cao, SCL Kamerlin, and A Warshel. Enzyme millisecond conformational dynamics do not catalyze the chemical step. Proceedings of the national academy of sciences of the United States of America, 106(41):17359{ 17364, 2009. [82] M Roca, B Messer, and A Warshel. Electrostatic contributions to protein stability and folding energy. FEBS letters, 581(10):2065{2071, 2007. [83] B Roux. In uence of the membrane potential on the free energy of an intrinsic protein. Biophysical journal, 73(6):2980{2989, 1997. [84] B Roux, T Allen, S Berneche, and W Im. Theoretical and computational models of biological ion channels. Quarterly reviews of biophysics, 37(1):15{ 103, 2004. [85] B Roux and M Karplus. Molecular-dynamics simulations of the gramicidin channel. Annual review of biophysics and biomolecular structure, 23:731{761, 1994. [86] B Roux and R MacKinnon. The cavity and pore helices the KcsA K+ channel: Electrostatic stabilization of monovalent cations. Science, 285(5424):100{102, 1999. [87] B Rudy. Diversity and ubiquity of K-channels. Neuroscience, 25(3):729{749, 1988. 106 [88] A Rychkova, S Vicatos, and A Warshel. On the energetics of translocon- assisted insertion of charged transmembrane helices into membranes. Pro- ceedings of the national academy of sciences of the United States oF America, 107(41):17598{17603, 2010. [89] Chung S, Anderson OS, and Krishnamurthy VV. Biological Membrane Ion Channels. Springer, 1 edition, 2010. [90] MSP Sansom, IH Shrivastava, JN Bright, J Tate, CE Capener, and PC Biggin. Potassium channels: structures, models, simulations. Biochimicaetbiophysica acta-biomembranes, 1565(2):294{307, 2002. [91] T Schlick, E Barth, and M Mandziuk. Biomolecular dynamics at long timesteps: Bridging the timescale gap between simulation and experimen- tation. Annual review of biophysics and biomolecular structure, 26:181{222, 1997. [92] T Schlick, RD Skeel, AT Brunger, LV Kale, JA Board, J Hermans, and K Schulten. Algorithmic challenges in computational molecular biophysics. Journal of computational physics, 151(1):9{48, 1999. [93] SA Seoh, D Sigg, DM Papazian, and F Bezanilla. Voltage-sensing residues in the S2 and S4 segments of the Shaker K+ channel. Neuron, 16(6):1159{1167, 1996. [94] YY Sham, ZT Chu, and A Warshel. Consistent calculations of pK(a)'s of ionizable residues in proteins: Semi-microscopic and microscopic approaches. Journal of physical chemistry B, 101(22):4458{4472, 1997. [95] IH Shrivastava and MSP Sansom. Simulations of ion permeation through a potassium channel: Molecular dynamics of KcsA in a phospholipid bilayer. Biophysical journal, 78(2):557{570, 2000. [96] D Sigg, F Bezanilla, and E Stefani. Fast gating in the shaker K+ channel and the energy landscape of activation. Proceedings of the national academy of sciences of the United States OF America, 100(13):7611{7615, 2003. [97] Y Sugita and Y Okamoto. Replica-exchange molecular dynamics method for protein folding. Chemical physics letters, 314(1-2):141{151, 1999. [98] C Tanford and R Roxby. Interpretation of protein titration curves - applica- tion to lysozyme. Biochemistry, 11(11):2192{2199, 1972. 107 [99] DP Tieleman and HJC Berendsen. A molecular dynamics study of the pores formed by Escherichia coli OmpF porin in a fully hydrated palmi- toyloleoylphosphatidylcholine bilayer. Biophysical journal, 74(6):2786{2801, 1998. [100] DP Tieleman, PC Biggin, GR Smith, and MSP Sansom. Simulation ap- proaches to ion channel structure-function relationships. Quarterly reviews of biophysics, 34(4):473{561, 2001. [101] GM Torrie and JP Valleau. Electrical double-layers .1. Monte-Carlo study of a uniformly charged surface. Journal of chemical physics, 73(11):5807{5816, 1980. [102] W Treptow and M Tarek. Environment of the gating charges in the Kv1.2 Shaker potassium channel. Biophysical journal, 90(9):L64{L66, 2006. [103] DJ Triggle, M Gopalakrishnan, and D Rampe. Voltage-gated ion channels as drug targets. Wiley-VCH, 1 edition, 2006. [104] D Trzesniak, AE Kunz, and WF van Gunsteren. A comparison of methods to compute the potential of mean force. ChemPhysChem, 8(1):162{169, 2007. [105] GE Uhlenbeck and LS Ornstein. On the theory of the Brownian motion. Physical review, 36(5):0823{0841, 1930. [106] WF Vangunsteren and HJC Berendsen. Computer-simulation of molecular- dynamics - methodology, applications, and perspectives in chemistry. Ange- wandte chemie-international edition in english, 29(9):992{1023, 1990. [107] S Vicatos, M Roca, and A Warshel. Eective approach for calculations of absolute stability of proteins using focused dielectric constants. Proteins- structure function and bioinformatics, 77(3):670{684, 2009. [108] J Villa and A Warshel. Energetics and dynamics of enzymatic reactions. Journal of physical chemistry B, 105(33):7887{7907, 2001. [109] A Warshel. Bicycle-pedal model for 1st step in vision process. Nature, 260(5553):679{683, 1976. [110] A Warshel. Computer Modeling of Chemical Reactions in Enzymes and Solu- tions. John Wiley & Sons, New York, 1991. [111] A Warshel. Molecular dynamics simulations of biological reactions. Accounts of chemical research, 35(6):385{395, 2002. 108 [112] A Warshel. Computer simulations of enzyme catalysis: Methods, progress, and insights. Annual review of biophysics and biomolecular structure, 32:425{ 443, 2003. [113] A Warshel and ST Russell. Calculations of electrostatic interactions in biological-systems and in solutions. Quarterly reviews of biophysics, 17(3):283{422, 1984. [114] A Warshel, PK Sharma, M Kato, and WW Parson. Modeling electrostatic eects in proteins. Biochimica et biophysica acta-proteins and proteomics, 1764(11):1647{1676, 2006. [115] JA White, JT Rubinstein, and AR Kay. Channel noise in neurons. Trends in neurosciences, 23(3):131{137, 2000. [116] YT Wong, TW Clark, J Shen, and JA McCammon. Molecular-dynamics simulation of substrate-enzyme interactions in the active-site channel of superoxide-dismutase. Molecular simulation, 10(2-6):277{288, 1993. [117] M Yi, H Nymeyer, and H Zhou. Test of the Gouy-Chapman theory for a charged lipid membrane against explicit-solvent molecular dynamics simula- tions. Physical review letters, 101(3):448{452, 2008. [118] YF Zhou, JH Morais-Cabral, A Kaufman, and R MacKinnon. Chemistry of ion coordination and hydration revealed by a K+ channel-Fab complex at 2.0 angstrom resolution. Nature, 414(6859):43{48, 2001. 109

Abstract (if available)

Abstract The relationship between the membrane potential and the gating of voltage-activated ion channels, as well as ion current through pores, is currently a problem of great interest. For example, simulations of conformational changes during the voltage-induced transition and ion current represent a major challenge that cannot be overcome at present by brute-force molecular dynamics approaches. To progress on this front, we developed and validated a renormalization approach and extended our previously developed coarse grained (CG) model so that it can describe the effect of the membrane potential and electrolyte solutions. Here we apply our model in realistic simulations of the energetics of the Kv1.2 channel. The simulations reproduce the observed experimental trend for both the gating charge and the effect of the membrane potential on the stability of the closed and open conformations. This indicates that our CG modeling offers a powerful tool for simulating voltage-activated channels.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Investigating voltage activated potassium and proton channels

PDF

Exploring the nature of the translocon-assisted protein insertion

PDF

The effects of divalent counterions on the formation and stabilization of RNA tertiary structure

PDF

Quantitative computer-aided studies of artificial and enantioselective enzymes

PDF

Advancing ab initio QM/MM free energy calculations: refining, validating and quantifying the reference potential approach

PDF

Shock-induced poration, cholesterol flip-flop and small interfering RNA transfection in a phospholipid membrane: multimillion atom, microsecond molecular dynamics simulations

PDF

Theoretical studies of lipid bilayer electroporation using molecular dynamics simulations

PDF

Simulating the helicase motor of SV40 large tumor antigen

PDF

Understanding electrostatic effects in the function of biological systems

PDF

Molecular simulations of water and monovalent ion dynamics in the electroporation of phospholipid bilayers

PDF

Photoexcitation and nonradiative relaxation in molecular systems: methodology, optical properties, and dynamics

PDF

The symmetry of membrane protein polyhedral nanoparticles

PDF

Investigation of the molecular mechanisms underlying polarized trafficking of the potassium channels Kv4.2 and Kv1.3 in neurons

PDF

Molecular dynamics and virial expansions for equilibrium thermodynamics in the solar interior

PDF

Two state free energy barrier evaluation technique for dehalogenation reaction

PDF

Advancing computational methods for free energy calculations in proteins and the applications for studies of catalysis in nucleotide-associated proteins

PDF

Modeling and simulation of complex recovery processes

PDF

Molecular dynamics simulations of lipid bilayers in megavolt per meter electric fields

PDF

Deterministic evolutionary game theory models for the nonlinear dynamics and control of chemotherapeutic resistance

PDF

Molecular dynamics simulations of nanostructures

Asset Metadata

Creator Dryga, Anatoly (author)

Core Title Capturing the complexity of ion channels: simulations of long-time dynamics, conformational changes and the effect of the membrane potential

School College of Letters, Arts and Sciences

Degree Doctor of Philosophy

Degree Program Chemistry

Publication Date 05/09/2012

Defense Date 03/20/2012

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag coarse-grain,Kv1.2,Langevin dynamics,molecular dynamics,OAI-PMH Harvest,renormalization

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Warshel, Arieh (committee chair), Nakano, Aiichiro (committee member)

Creator Email anatolydryga@gmail.com,dryga@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c3-37354

Unique identifier UC11289494

Identifier usctheses-c3-37354 (legacy record id)

Legacy Identifier etd-DrygaAnato-824.pdf

Dmrecord 37354

Document Type Dissertation

Rights Dryga, Anatoly

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA